Merge branch 'develop' of https://github.com/paboyle/Grid into develop

.gitignore

@@ -93,6 +93,7 @@ build*/*
 *.xcodeproj/*
 build.sh
 .vscode
+*.code-workspace
 
 # Eigen source #
 ################
@@ -122,4 +123,3 @@ make-bin-BUCK.sh
 #####################
 lib/qcd/spin/gamma-gen/*.h
 lib/qcd/spin/gamma-gen/*.cc
-

extras/Hadrons/Application.cc

@@ -73,12 +73,6 @@ Application::Application(const std::string parameterFileName)
     parameterFileName_ = parameterFileName;
 }
 
-// environment shortcut ////////////////////////////////////////////////////////
-Environment & Application::env(void) const
-{
-    return Environment::getInstance();
-}
-
 // access //////////////////////////////////////////////////////////////////////
 void Application::setPar(const Application::GlobalPar &par)
 {
@@ -94,14 +88,13 @@ const Application::GlobalPar & Application::getPar(void)
 // execute /////////////////////////////////////////////////////////////////////
 void Application::run(void)
 {
-    if (!parameterFileName_.empty() and (env().getNModule() == 0))
+    if (!parameterFileName_.empty() and (vm().getNModule() == 0))
     {
         parseParameterFile(parameterFileName_);
     }
-    if (!scheduled_)
+    vm().printContent();
-    {
+    env().printContent();
-        schedule();
+    schedule();
-    }
     printSchedule();
     configLoop();
 }
@@ -124,12 +117,20 @@ void Application::parseParameterFile(const std::string parameterFileName)
     LOG(Message) << "Building application from '" << parameterFileName << "'..." << std::endl;
     read(reader, "parameters", par);
     setPar(par);
-    push(reader, "modules");
+    if (!push(reader, "modules"))
-    push(reader, "module");
+    {
+        HADRON_ERROR(Parsing, "Cannot open node 'modules' in parameter file '" 
+                              + parameterFileName + "'");
+    }
+    if (!push(reader, "module"))
+    {
+        HADRON_ERROR(Parsing, "Cannot open node 'modules/module' in parameter file '" 
+                              + parameterFileName + "'");
+    }
     do
     {
         read(reader, "id", id);
-        env().createModule(id.name, id.type, reader);
+        vm().createModule(id.name, id.type, reader);
     } while (reader.nextElement("module"));
     pop(reader);
     pop(reader);
@@ -139,7 +140,7 @@ void Application::saveParameterFile(const std::string parameterFileName)
 {
     XmlWriter          writer(parameterFileName);
     ObjectId           id;
-    const unsigned int nMod = env().getNModule();
+    const unsigned int nMod = vm().getNModule();
     
     LOG(Message) << "Saving application to '" << parameterFileName << "'..." << std::endl;
     write(writer, "parameters", getPar());
@@ -147,10 +148,10 @@ void Application::saveParameterFile(const std::string parameterFileName)
     for (unsigned int i = 0; i < nMod; ++i)
     {
         push(writer, "module");
-        id.name = env().getModuleName(i);
+        id.name = vm().getModuleName(i);
-        id.type = env().getModule(i)->getRegisteredName();
+        id.type = vm().getModule(i)->getRegisteredName();
         write(writer, "id", id);
-        env().getModule(i)->saveParameters(writer, "options");
+        vm().getModule(i)->saveParameters(writer, "options");
         pop(writer);
     }
     pop(writer);
@@ -158,95 +159,13 @@ void Application::saveParameterFile(const std::string parameterFileName)
 }
 
 // schedule computation ////////////////////////////////////////////////////////
-#define MEM_MSG(size)\
-sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
-
-#define DEFINE_MEMPEAK \
-GeneticScheduler<unsigned int>::ObjFunc memPeak = \
-[this](const std::vector<unsigned int> &program)\
-{\
-    unsigned int memPeak;\
-    bool         msg;\
-    \
-    msg = HadronsLogMessage.isActive();\
-    HadronsLogMessage.Active(false);\
-    env().dryRun(true);\
-    memPeak = env().executeProgram(program);\
-    env().dryRun(false);\
-    env().freeAll();\
-    HadronsLogMessage.Active(true);\
-    \
-    return memPeak;\
-}
-
 void Application::schedule(void)
 {
-    DEFINE_MEMPEAK;
+    if (!scheduled_ and !loadedSchedule_)
-    
-    // build module dependency graph
-    LOG(Message) << "Building module graph..." << std::endl;
-    auto graph = env().makeModuleGraph();
-    auto con = graph.getConnectedComponents();
-    
-    // constrained topological sort using a genetic algorithm
-    LOG(Message) << "Scheduling computation..." << std::endl;
-    LOG(Message) << "               #module= " << graph.size() << std::endl;
-    LOG(Message) << "       population size= " << par_.genetic.popSize << std::endl;
-    LOG(Message) << "       max. generation= " << par_.genetic.maxGen << std::endl;
-    LOG(Message) << "  max. cst. generation= " << par_.genetic.maxCstGen << std::endl;
-    LOG(Message) << "         mutation rate= " << par_.genetic.mutationRate << std::endl;
-    
-    unsigned int                               k = 0, gen, prevPeak, nCstPeak = 0;
-    std::random_device                         rd;
-    GeneticScheduler<unsigned int>::Parameters par;
-    
-    par.popSize      = par_.genetic.popSize;
-    par.mutationRate = par_.genetic.mutationRate;
-    par.seed         = rd();
-    memPeak_         = 0;
-    CartesianCommunicator::BroadcastWorld(0, &(par.seed), sizeof(par.seed));
-    for (unsigned int i = 0; i < con.size(); ++i)
     {
-        GeneticScheduler<unsigned int> scheduler(con[i], memPeak, par);
+        program_   = vm().schedule(par_.genetic);
-        
+        scheduled_ = true;
-        gen = 0;
-        do
-        {
-            LOG(Debug) << "Generation " << gen << ":" << std::endl;
-            scheduler.nextGeneration();
-            if (gen != 0)
-            {
-                if (prevPeak == scheduler.getMinValue())
-                {
-                    nCstPeak++;
-                }
-                else
-                {
-                    nCstPeak = 0;
-                }
-            }
-            
-            prevPeak = scheduler.getMinValue();
-            if (gen % 10 == 0)
-            {
-                LOG(Iterative) << "Generation " << gen << ": "
-                               << MEM_MSG(scheduler.getMinValue()) << std::endl;
-            }
-            
-            gen++;
-        } while ((gen < par_.genetic.maxGen)
-                 and (nCstPeak < par_.genetic.maxCstGen));
-        auto &t = scheduler.getMinSchedule();
-        if (scheduler.getMinValue() > memPeak_)
-        {
-            memPeak_ = scheduler.getMinValue();
-        }
-        for (unsigned int j = 0; j < t.size(); ++j)
-        {
-            program_.push_back(t[j]);
-        }
     }
-    scheduled_ = true;
 }
 
 void Application::saveSchedule(const std::string filename)
@@ -256,21 +175,19 @@ void Application::saveSchedule(const std::string filename)
     
     if (!scheduled_)
     {
-        HADRON_ERROR("Computation not scheduled");
+        HADRON_ERROR(Definition, "Computation not scheduled");
     }
     LOG(Message) << "Saving current schedule to '" << filename << "'..."
                  << std::endl;
     for (auto address: program_)
     {
-        program.push_back(env().getModuleName(address));
+        program.push_back(vm().getModuleName(address));
     }
     write(writer, "schedule", program);
 }
 
 void Application::loadSchedule(const std::string filename)
 {
-    DEFINE_MEMPEAK;
-    
     TextReader               reader(filename);
     std::vector<std::string> program;
     
@@ -280,24 +197,24 @@ void Application::loadSchedule(const std::string filename)
     program_.clear();
     for (auto &name: program)
     {
-        program_.push_back(env().getModuleAddress(name));
+        program_.push_back(vm().getModuleAddress(name));
     }
-    scheduled_ = true;
+    loadedSchedule_ = true;
-    memPeak_   = memPeak(program_);
 }
 
 void Application::printSchedule(void)
 {
     if (!scheduled_)
     {
-        HADRON_ERROR("Computation not scheduled");
+        HADRON_ERROR(Definition, "Computation not scheduled");
     }
-    LOG(Message) << "Schedule (memory peak: " << MEM_MSG(memPeak_) << "):"
+    auto peak = vm().memoryNeeded(program_);
+    LOG(Message) << "Schedule (memory needed: " << sizeString(peak) << "):"
                  << std::endl;
     for (unsigned int i = 0; i < program_.size(); ++i)
     {
         LOG(Message) << std::setw(4) << i + 1 << ": "
-                     << env().getModuleName(program_[i]) << std::endl;
+                     << vm().getModuleName(program_[i]) << std::endl;
     }
 }
 
@@ -310,8 +227,8 @@ void Application::configLoop(void)
     {
         LOG(Message) << BIG_SEP << " Starting measurement for trajectory " << t
                      << " " << BIG_SEP << std::endl;
-        env().setTrajectory(t);
+        vm().setTrajectory(t);
-        env().executeProgram(program_);
+        vm().executeProgram(program_);
     }
     LOG(Message) << BIG_SEP << " End of measurement " << BIG_SEP << std::endl;
     env().freeAll();

extras/Hadrons/Application.hpp

@@ -31,8 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define Hadrons_Application_hpp_
 
 #include <Grid/Hadrons/Global.hpp>
-#include <Grid/Hadrons/Environment.hpp>
+#include <Grid/Hadrons/VirtualMachine.hpp>
-#include <Grid/Hadrons/ModuleFactory.hpp>
 #include <Grid/Hadrons/Modules.hpp>
 
 BEGIN_HADRONS_NAMESPACE
@@ -51,25 +50,13 @@ public:
                                         unsigned int, end,
                                         unsigned int, step);
     };
-    class GeneticPar: Serializable
-    {
-    public:
-        GeneticPar(void):
-            popSize{20}, maxGen{1000}, maxCstGen{100}, mutationRate{.1} {};
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(GeneticPar,
-                                        unsigned int, popSize,
-                                        unsigned int, maxGen,
-                                        unsigned int, maxCstGen,
-                                        double      , mutationRate);
-    };
     class GlobalPar: Serializable
     {
     public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(GlobalPar,
-                                        TrajRange,   trajCounter,
+                                        TrajRange,                  trajCounter,
-                                        GeneticPar,  genetic,
+                                        VirtualMachine::GeneticPar, genetic,
-                                        std::string, seed);
+                                        std::string,                seed);
     };
 public:
     // constructors
@@ -100,14 +87,15 @@ public:
     void configLoop(void);
 private:
     // environment shortcut
-    Environment & env(void) const;
+    DEFINE_ENV_ALIAS;
+    // virtual machine shortcut
+    DEFINE_VM_ALIAS;
 private:
-    long unsigned int         locVol_;
+    long unsigned int       locVol_;
-    std::string               parameterFileName_{""};
+    std::string             parameterFileName_{""};
-    GlobalPar                 par_;
+    GlobalPar               par_;
-    std::vector<unsigned int> program_;
+    VirtualMachine::Program program_;
-    Environment::Size         memPeak_;
+    bool                    scheduled_{false}, loadedSchedule_{false};
-    bool                      scheduled_{false};
 };
 
 /******************************************************************************
@@ -117,14 +105,16 @@ private:
 template <typename M>
 void Application::createModule(const std::string name)
 {
-    env().createModule<M>(name);
+    vm().createModule<M>(name);
+    scheduled_ = false;
 }
 
 template <typename M>
 void Application::createModule(const std::string name,
                                const typename M::Par &par)
 {
-    env().createModule<M>(name, par);
+    vm().createModule<M>(name, par);
+    scheduled_ = false;
 }
 
 END_HADRONS_NAMESPACE

extras/Hadrons/Environment.cc

@@ -35,6 +35,9 @@ using namespace Grid;
 using namespace QCD;
 using namespace Hadrons;
 
+#define ERROR_NO_ADDRESS(address)\
+HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
+
 /******************************************************************************
  *                       Environment implementation                           *
  ******************************************************************************/
@@ -56,28 +59,6 @@ Environment::Environment(void)
     rng4d_.reset(new GridParallelRNG(grid4d_.get()));
 }
 
-// dry run /////////////////////////////////////////////////////////////////////
-void Environment::dryRun(const bool isDry)
-{
-    dryRun_ = isDry;
-}
-
-bool Environment::isDryRun(void) const
-{
-    return dryRun_;
-}
-
-// trajectory number ///////////////////////////////////////////////////////////
-void Environment::setTrajectory(const unsigned int traj)
-{
-    traj_ = traj;
-}
-
-unsigned int Environment::getTrajectory(void) const
-{
-    return traj_;
-}
-
 // grids ///////////////////////////////////////////////////////////////////////
 void Environment::createGrid(const unsigned int Ls)
 {
@@ -105,7 +86,7 @@ GridCartesian * Environment::getGrid(const unsigned int Ls) const
     }
     catch(std::out_of_range &)
     {
-        HADRON_ERROR("no grid with Ls= " << Ls);
+        HADRON_ERROR(Definition, "no grid with Ls= " + std::to_string(Ls));
     }
 }
 
@@ -124,7 +105,7 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
     }
     catch(std::out_of_range &)
     {
-        HADRON_ERROR("no red-black 5D grid with Ls= " << Ls);
+        HADRON_ERROR(Definition, "no red-black 5D grid with Ls= " + std::to_string(Ls));
     }
 }
 
@@ -143,6 +124,11 @@ int Environment::getDim(const unsigned int mu) const
     return dim_[mu];
 }
 
+unsigned long int Environment::getLocalVolume(void) const
+{
+    return locVol_;
+}
+
 // random number generator /////////////////////////////////////////////////////
 void Environment::setSeed(const std::vector<int> &seed)
 {
@@ -154,291 +140,6 @@ GridParallelRNG * Environment::get4dRng(void) const
     return rng4d_.get();
 }
 
-// module management ///////////////////////////////////////////////////////////
-void Environment::pushModule(Environment::ModPt &pt)
-{
-    std::string name = pt->getName();
-    
-    if (!hasModule(name))
-    {
-        std::vector<unsigned int> inputAddress;
-        unsigned int              address;
-        ModuleInfo                m;
-        
-        m.data = std::move(pt);
-        m.type = typeIdPt(*m.data.get());
-        m.name = name;
-        auto input  = m.data->getInput();
-        for (auto &in: input)
-        {
-            if (!hasObject(in))
-            {
-                addObject(in , -1);
-            }
-            m.input.push_back(objectAddress_[in]);
-        }
-        auto output = m.data->getOutput();
-        module_.push_back(std::move(m));
-        address              = static_cast<unsigned int>(module_.size() - 1);
-        moduleAddress_[name] = address;
-        for (auto &out: output)
-        {
-            if (!hasObject(out))
-            {
-                addObject(out, address);
-            }
-            else
-            {
-                if (object_[objectAddress_[out]].module < 0)
-                {
-                    object_[objectAddress_[out]].module = address;
-                }
-                else
-                {
-                    HADRON_ERROR("object '" + out
-                                 + "' is already produced by module '"
-                                 + module_[object_[getObjectAddress(out)].module].name
-                                 + "' (while pushing module '" + name + "')");
-                }
-            }
-        }
-    }
-    else
-    {
-        HADRON_ERROR("module '" + name + "' already exists");
-    }
-}
-
-unsigned int Environment::getNModule(void) const
-{
-    return module_.size();
-}
-
-void Environment::createModule(const std::string name, const std::string type,
-                               XmlReader &reader)
-{
-    auto &factory = ModuleFactory::getInstance();
-    auto pt       = factory.create(type, name);
-    
-    pt->parseParameters(reader, "options");
-    pushModule(pt);
-}
-
-ModuleBase * Environment::getModule(const unsigned int address) const
-{
-    if (hasModule(address))
-    {
-        return module_[address].data.get();
-    }
-    else
-    {
-        HADRON_ERROR("no module with address " + std::to_string(address));
-    }
-}
-
-ModuleBase * Environment::getModule(const std::string name) const
-{
-    return getModule(getModuleAddress(name));
-}
-
-unsigned int Environment::getModuleAddress(const std::string name) const
-{
-    if (hasModule(name))
-    {
-        return moduleAddress_.at(name);
-    }
-    else
-    {
-        HADRON_ERROR("no module with name '" + name + "'");
-    }
-}
-
-std::string Environment::getModuleName(const unsigned int address) const
-{
-    if (hasModule(address))
-    {
-        return module_[address].name;
-    }
-    else
-    {
-        HADRON_ERROR("no module with address " + std::to_string(address));
-    }
-}
-
-std::string Environment::getModuleType(const unsigned int address) const
-{
-    if (hasModule(address))
-    {
-        return typeName(module_[address].type);
-    }
-    else
-    {
-        HADRON_ERROR("no module with address " + std::to_string(address));
-    }
-}
-
-std::string Environment::getModuleType(const std::string name) const
-{
-    return getModuleType(getModuleAddress(name));
-}
-
-std::string Environment::getModuleNamespace(const unsigned int address) const
-{
-    std::string type = getModuleType(address), ns;
-    
-    auto pos2 = type.rfind("::");
-    auto pos1 = type.rfind("::", pos2 - 2);
-    
-    return type.substr(pos1 + 2, pos2 - pos1 - 2);
-}
-
-std::string Environment::getModuleNamespace(const std::string name) const
-{
-    return getModuleNamespace(getModuleAddress(name));
-}
-
-bool Environment::hasModule(const unsigned int address) const
-{
-    return (address < module_.size());
-}
-
-bool Environment::hasModule(const std::string name) const
-{
-    return (moduleAddress_.find(name) != moduleAddress_.end());
-}
-
-Graph<unsigned int> Environment::makeModuleGraph(void) const
-{
-    Graph<unsigned int> moduleGraph;
-    
-    for (unsigned int i = 0; i < module_.size(); ++i)
-    {
-        moduleGraph.addVertex(i);
-        for (auto &j: module_[i].input)
-        {
-            moduleGraph.addEdge(object_[j].module, i);
-        }
-    }
-    
-    return moduleGraph;
-}
-
-#define BIG_SEP "==============="
-#define SEP     "---------------"
-#define MEM_MSG(size)\
-sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
-
-Environment::Size
-Environment::executeProgram(const std::vector<unsigned int> &p)
-{
-    Size                                memPeak = 0, sizeBefore, sizeAfter;
-    std::vector<std::set<unsigned int>> freeProg;
-    bool                                continueCollect, nothingFreed;
-    
-    // build garbage collection schedule
-    freeProg.resize(p.size());
-    for (unsigned int i = 0; i < object_.size(); ++i)
-    {
-        auto pred = [i, this](const unsigned int j)
-        {
-            auto &in = module_[j].input;
-            auto it  = std::find(in.begin(), in.end(), i);
-            
-            return (it != in.end()) or (j == object_[i].module);
-        };
-        auto it = std::find_if(p.rbegin(), p.rend(), pred);
-        if (it != p.rend())
-        {
-            freeProg[p.rend() - it - 1].insert(i);
-        }
-    }
-    
-    // program execution
-    for (unsigned int i = 0; i < p.size(); ++i)
-    {
-        // execute module
-        if (!isDryRun())
-        {
-            LOG(Message) << SEP << " Measurement step " << i+1 << "/"
-                         << p.size() << " (module '" << module_[p[i]].name
-                         << "') " << SEP << std::endl;
-        }
-        (*module_[p[i]].data)();
-        sizeBefore = getTotalSize();
-        // print used memory after execution
-        if (!isDryRun())
-        {
-            LOG(Message) << "Allocated objects: " << MEM_MSG(sizeBefore)
-                         << std::endl;
-        }
-        if (sizeBefore > memPeak)
-        {
-            memPeak = sizeBefore;
-        }
-        // garbage collection for step i
-        if (!isDryRun())
-        {
-            LOG(Message) << "Garbage collection..." << std::endl;
-        }
-        nothingFreed = true;
-        do
-        {
-            continueCollect = false;
-            auto toFree = freeProg[i];
-            for (auto &j: toFree)
-            {
-                // continue garbage collection while there are still
-                // objects without owners
-                continueCollect = continueCollect or !hasOwners(j);
-                if(freeObject(j))
-                {
-                    // if an object has been freed, remove it from
-                    // the garbage collection schedule
-                    freeProg[i].erase(j);
-                    nothingFreed = false;
-                }
-            }
-        } while (continueCollect);
-        // any remaining objects in step i garbage collection schedule
-        // is scheduled for step i + 1
-        if (i + 1 < p.size())
-        {
-            for (auto &j: freeProg[i])
-            {
-                freeProg[i + 1].insert(j);
-            }
-        }
-        // print used memory after garbage collection if necessary
-        if (!isDryRun())
-        {
-            sizeAfter = getTotalSize();
-            if (sizeBefore != sizeAfter)
-            {
-                LOG(Message) << "Allocated objects: " << MEM_MSG(sizeAfter)
-                             << std::endl;
-            }
-            else
-            {
-                LOG(Message) << "Nothing to free" << std::endl;
-            }
-        }
-    }
-    
-    return memPeak;
-}
-
-Environment::Size Environment::executeProgram(const std::vector<std::string> &p)
-{
-    std::vector<unsigned int> pAddress;
-    
-    for (auto &n: p)
-    {
-        pAddress.push_back(getModuleAddress(n));
-    }
-    
-    return executeProgram(pAddress);
-}
-
 // general memory management ///////////////////////////////////////////////////
 void Environment::addObject(const std::string name, const int moduleAddress)
 {
@@ -448,46 +149,25 @@ void Environment::addObject(const std::string name, const int moduleAddress)
         
         info.name   = name;
         info.module = moduleAddress;
+        info.data   = nullptr;
         object_.push_back(std::move(info));
         objectAddress_[name] = static_cast<unsigned int>(object_.size() - 1);
     }
     else
     {
-        HADRON_ERROR("object '" + name + "' already exists");
+        HADRON_ERROR(Definition, "object '" + name + "' already exists");
     }
 }
 
-void Environment::registerObject(const unsigned int address,
+void Environment::setObjectModule(const unsigned int objAddress,
-                                 const unsigned int size, const unsigned int Ls)
+                                  const int modAddress)
 {
-    if (!hasRegisteredObject(address))
+    object_[objAddress].module = modAddress;
-    {
-        if (hasObject(address))
-        {
-            object_[address].size         = size;
-            object_[address].Ls           = Ls;
-            object_[address].isRegistered = true;
-        }
-        else
-        {
-            HADRON_ERROR("no object with address " + std::to_string(address));
-        }
-    }
-    else
-    {
-        HADRON_ERROR("object with address " + std::to_string(address)
-                     + " already registered");
-    }
 }
 
-void Environment::registerObject(const std::string name,
+unsigned int Environment::getMaxAddress(void) const
-                                 const unsigned int size, const unsigned int Ls)
 {
-    if (!hasObject(name))
+    return object_.size();
-    {
-        addObject(name);
-    }
-    registerObject(getObjectAddress(name), size, Ls);
 }
 
 unsigned int Environment::getObjectAddress(const std::string name) const
@@ -498,7 +178,7 @@ unsigned int Environment::getObjectAddress(const std::string name) const
     }
     else
     {
-        HADRON_ERROR("no object with name '" + name + "'");
+        HADRON_ERROR(Definition, "no object with name '" + name + "'");
     }
 }
 
@@ -510,13 +190,13 @@ std::string Environment::getObjectName(const unsigned int address) const
     }
     else
     {
-        HADRON_ERROR("no object with address " + std::to_string(address));
+        ERROR_NO_ADDRESS(address);
     }
 }
 
 std::string Environment::getObjectType(const unsigned int address) const
 {
-    if (hasRegisteredObject(address))
+    if (hasObject(address))
     {
         if (object_[address].type)
         {
@@ -527,14 +207,9 @@ std::string Environment::getObjectType(const unsigned int address) const
             return "<no type>";
         }
     }
-    else if (hasObject(address))
-    {
-        HADRON_ERROR("object with address " + std::to_string(address)
-                     + " exists but is not registered");
-    }
     else
     {
-        HADRON_ERROR("no object with address " + std::to_string(address));
+        ERROR_NO_ADDRESS(address);
     }
 }
 
@@ -545,18 +220,13 @@ std::string Environment::getObjectType(const std::string name) const
 
 Environment::Size Environment::getObjectSize(const unsigned int address) const
 {
-    if (hasRegisteredObject(address))
+    if (hasObject(address))
     {
         return object_[address].size;
     }
-    else if (hasObject(address))
-    {
-        HADRON_ERROR("object with address " + std::to_string(address)
-                     + " exists but is not registered");
-    }
     else
     {
-        HADRON_ERROR("no object with address " + std::to_string(address));
+        ERROR_NO_ADDRESS(address);
     }
 }
 
@@ -565,7 +235,24 @@ Environment::Size Environment::getObjectSize(const std::string name) const
     return getObjectSize(getObjectAddress(name));
 }
 
-unsigned int Environment::getObjectModule(const unsigned int address) const
+Environment::Storage Environment::getObjectStorage(const unsigned int address) const
+{
+    if (hasObject(address))
+    {
+        return object_[address].storage;
+    }
+    else
+    {
+        ERROR_NO_ADDRESS(address);
+    }
+}
+
+Environment::Storage Environment::getObjectStorage(const std::string name) const
+{
+    return getObjectStorage(getObjectAddress(name));
+}
+
+int Environment::getObjectModule(const unsigned int address) const
 {
     if (hasObject(address))
     {
@@ -573,29 +260,24 @@ unsigned int Environment::getObjectModule(const unsigned int address) const
     }
     else
     {
-        HADRON_ERROR("no object with address " + std::to_string(address));
+        ERROR_NO_ADDRESS(address);
     }
 }
 
-unsigned int Environment::getObjectModule(const std::string name) const
+int Environment::getObjectModule(const std::string name) const
 {
     return getObjectModule(getObjectAddress(name));
 }
 
 unsigned int Environment::getObjectLs(const unsigned int address) const
 {
-    if (hasRegisteredObject(address))
+    if (hasObject(address))
     {
         return object_[address].Ls;
     }
-    else if (hasObject(address))
-    {
-        HADRON_ERROR("object with address " + std::to_string(address)
-                     + " exists but is not registered");
-    }
     else
     {
-        HADRON_ERROR("no object with address " + std::to_string(address));
+        ERROR_NO_ADDRESS(address);
     }
 }
 
@@ -616,30 +298,6 @@ bool Environment::hasObject(const std::string name) const
     return ((it != objectAddress_.end()) and hasObject(it->second));
 }
 
-bool Environment::hasRegisteredObject(const unsigned int address) const
-{
-    if (hasObject(address))
-    {
-        return object_[address].isRegistered;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-bool Environment::hasRegisteredObject(const std::string name) const
-{
-    if (hasObject(name))
-    {
-        return hasRegisteredObject(getObjectAddress(name));
-    }
-    else
-    {
-        return false;
-    }
-}
-
 bool Environment::hasCreatedObject(const unsigned int address) const
 {
     if (hasObject(address))
@@ -680,92 +338,27 @@ Environment::Size Environment::getTotalSize(void) const
     
     for (auto &o: object_)
     {
-        if (o.isRegistered)
+        size += o.size;
-        {
-            size += o.size;
-        }
     }
     
     return size;
 }
 
-void Environment::addOwnership(const unsigned int owner,
+void Environment::freeObject(const unsigned int address)
-                               const unsigned int property)
 {
-    if (hasObject(property))
+    if (hasCreatedObject(address))
     {
-        object_[property].owners.insert(owner);
+        LOG(Message) << "Destroying object '" << object_[address].name
-    }
+                     << "'" << std::endl;
-    else
-    {
-        HADRON_ERROR("no object with address " + std::to_string(property));
-    }
-    if (hasObject(owner))
-    {
-        object_[owner].properties.insert(property);
-    }
-    else
-    {
-        HADRON_ERROR("no object with address " + std::to_string(owner));
     }
+    object_[address].size = 0;
+    object_[address].type = nullptr;
+    object_[address].data.reset(nullptr);
 }
 
-void Environment::addOwnership(const std::string owner,
+void Environment::freeObject(const std::string name)
-                               const std::string property)
 {
-    addOwnership(getObjectAddress(owner), getObjectAddress(property));
+    freeObject(getObjectAddress(name));
-}
-
-bool Environment::hasOwners(const unsigned int address) const
-{
-    
-    if (hasObject(address))
-    {
-        return (!object_[address].owners.empty());
-    }
-    else
-    {
-        HADRON_ERROR("no object with address " + std::to_string(address));
-    }
-}
-
-bool Environment::hasOwners(const std::string name) const
-{
-    return hasOwners(getObjectAddress(name));
-}
-
-bool Environment::freeObject(const unsigned int address)
-{
-    if (!hasOwners(address))
-    {
-        if (!isDryRun() and object_[address].isRegistered)
-        {
-            LOG(Message) << "Destroying object '" << object_[address].name
-                         << "'" << std::endl;
-        }
-        for (auto &p: object_[address].properties)
-        {
-            object_[p].owners.erase(address);
-        }
-        object_[address].size         = 0;
-        object_[address].Ls           = 0;
-        object_[address].isRegistered = false;
-        object_[address].type         = nullptr;
-        object_[address].owners.clear();
-        object_[address].properties.clear();
-        object_[address].data.reset(nullptr);
-        
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-bool Environment::freeObject(const std::string name)
-{
-    return freeObject(getObjectAddress(name));
 }
 
 void Environment::freeAll(void)
@@ -776,18 +369,24 @@ void Environment::freeAll(void)
     }
 }
 
-void Environment::printContent(void)
+void Environment::protectObjects(const bool protect)
 {
-    LOG(Message) << "Modules: " << std::endl;
+    protect_ = protect;
-    for (unsigned int i = 0; i < module_.size(); ++i)
+}
-    {
+
-        LOG(Message) << std::setw(4) << i << ": "
+bool Environment::objectsProtected(void) const
-                     << getModuleName(i) << std::endl;
+{
-    }
+    return protect_;
-    LOG(Message) << "Objects: " << std::endl;
+}
+
+// print environment content ///////////////////////////////////////////////////
+void Environment::printContent(void) const
+{
+    LOG(Debug) << "Objects: " << std::endl;
     for (unsigned int i = 0; i < object_.size(); ++i)
     {
-        LOG(Message) << std::setw(4) << i << ": "
+        LOG(Debug) << std::setw(4) << i << ": "
-                     << getObjectName(i) << std::endl;
+                   << getObjectName(i) << " ("
+                   << sizeString(getObjectSize(i)) << ")" << std::endl;
     }
 }

extras/Hadrons/Environment.hpp

@@ -31,20 +31,12 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define Hadrons_Environment_hpp_
 
 #include <Grid/Hadrons/Global.hpp>
-#include <Grid/Hadrons/Graph.hpp>
-
-#ifndef SITE_SIZE_TYPE
-#define SITE_SIZE_TYPE unsigned int
-#endif
 
 BEGIN_HADRONS_NAMESPACE
 
 /******************************************************************************
  *                         Global environment                                 *
  ******************************************************************************/
-// forward declaration of Module
-class ModuleBase;
-
 class Object
 {
 public:
@@ -66,123 +58,78 @@ private:
     std::unique_ptr<T> objPt_{nullptr};
 };
 
+#define DEFINE_ENV_ALIAS \
+inline Environment & env(void) const\
+{\
+    return Environment::getInstance();\
+}
+
 class Environment
 {
     SINGLETON(Environment);
 public:
     typedef SITE_SIZE_TYPE                         Size;
-    typedef std::unique_ptr<ModuleBase>            ModPt;
     typedef std::unique_ptr<GridCartesian>         GridPt;
     typedef std::unique_ptr<GridRedBlackCartesian> GridRbPt;
     typedef std::unique_ptr<GridParallelRNG>       RngPt;
-    typedef std::unique_ptr<LatticeBase>           LatticePt;
+    enum class Storage {object, cache, temporary};
 private:
-    struct ModuleInfo
-    {
-        const std::type_info      *type{nullptr};
-        std::string               name;
-        ModPt                     data{nullptr};
-        std::vector<unsigned int> input;
-    };
     struct ObjInfo
     {
         Size                    size{0};
+        Storage                 storage{Storage::object};
         unsigned int            Ls{0};
-        bool                    isRegistered{false};
         const std::type_info    *type{nullptr};
         std::string             name;
         int                     module{-1};
-        std::set<unsigned int>  owners, properties;
         std::unique_ptr<Object> data{nullptr};
     };
 public:
-    // dry run
-    void                    dryRun(const bool isDry);
-    bool                    isDryRun(void) const;
-    // trajectory number
-    void                    setTrajectory(const unsigned int traj);
-    unsigned int            getTrajectory(void) const;
     // grids
     void                    createGrid(const unsigned int Ls);
     GridCartesian *         getGrid(const unsigned int Ls = 1) const;
     GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
     std::vector<int>        getDim(void) const;
     int                     getDim(const unsigned int mu) const;
+    unsigned long int       getLocalVolume(void) const;
     unsigned int            getNd(void) const;
     // random number generator
     void                    setSeed(const std::vector<int> &seed);
     GridParallelRNG *       get4dRng(void) const;
-    // module management
-    void                    pushModule(ModPt &pt);
-    template <typename M>
-    void                    createModule(const std::string name);
-    template <typename M>
-    void                    createModule(const std::string name,
-                                         const typename M::Par &par);
-    void                    createModule(const std::string name,
-                                         const std::string type,
-                                         XmlReader &reader);
-    unsigned int            getNModule(void) const;
-    ModuleBase *            getModule(const unsigned int address) const;
-    ModuleBase *            getModule(const std::string name) const;
-    template <typename M>
-    M *                     getModule(const unsigned int address) const;
-    template <typename M>
-    M *                     getModule(const std::string name) const;
-    unsigned int            getModuleAddress(const std::string name) const;
-    std::string             getModuleName(const unsigned int address) const;
-    std::string             getModuleType(const unsigned int address) const;
-    std::string             getModuleType(const std::string name) const;
-    std::string             getModuleNamespace(const unsigned int address) const;
-    std::string             getModuleNamespace(const std::string name) const;
-    bool                    hasModule(const unsigned int address) const;
-    bool                    hasModule(const std::string name) const;
-    Graph<unsigned int>     makeModuleGraph(void) const;
-    Size                    executeProgram(const std::vector<unsigned int> &p);
-    Size                    executeProgram(const std::vector<std::string> &p);
     // general memory management
     void                    addObject(const std::string name,
                                       const int moduleAddress = -1);
-    void                    registerObject(const unsigned int address,
+    template <typename B, typename T, typename ... Ts>
-                                           const unsigned int size,
+    void                    createDerivedObject(const std::string name,
-                                           const unsigned int Ls = 1);
+                                                const Environment::Storage storage,
-    void                    registerObject(const std::string name,
+                                                const unsigned int Ls,
-                                           const unsigned int size,
+                                                Ts && ... args);
-                                           const unsigned int Ls = 1);
+    template <typename T, typename ... Ts>
-    template <typename T>
+    void                    createObject(const std::string name,
-    unsigned int            lattice4dSize(void) const;
+                                         const Environment::Storage storage,
-    template <typename T>
+                                         const unsigned int Ls,
-    void                    registerLattice(const unsigned int address,
+                                         Ts && ... args);
-                                            const unsigned int Ls = 1);
+    void                    setObjectModule(const unsigned int objAddress,
-    template <typename T>
+                                            const int modAddress);
-    void                    registerLattice(const std::string name,
-                                            const unsigned int Ls = 1);
-    template <typename T>
-    void                    setObject(const unsigned int address, T *object);
-    template <typename T>
-    void                    setObject(const std::string name, T *object);
     template <typename T>
     T *                     getObject(const unsigned int address) const;
     template <typename T>
     T *                     getObject(const std::string name) const;
-    template <typename T>
+    unsigned int            getMaxAddress(void) const;
-    T *                     createLattice(const unsigned int address);
-    template <typename T>
-    T *                     createLattice(const std::string name);
     unsigned int            getObjectAddress(const std::string name) const;
     std::string             getObjectName(const unsigned int address) const;
     std::string             getObjectType(const unsigned int address) const;
     std::string             getObjectType(const std::string name) const;
     Size                    getObjectSize(const unsigned int address) const;
     Size                    getObjectSize(const std::string name) const;
-    unsigned int            getObjectModule(const unsigned int address) const;
+    Storage                 getObjectStorage(const unsigned int address) const;
-    unsigned int            getObjectModule(const std::string name) const;
+    Storage                 getObjectStorage(const std::string name) const;
+    int                     getObjectModule(const unsigned int address) const;
+    int                     getObjectModule(const std::string name) const;
     unsigned int            getObjectLs(const unsigned int address) const;
     unsigned int            getObjectLs(const std::string name) const;
     bool                    hasObject(const unsigned int address) const;
     bool                    hasObject(const std::string name) const;
-    bool                    hasRegisteredObject(const unsigned int address) const;
-    bool                    hasRegisteredObject(const std::string name) const;
     bool                    hasCreatedObject(const unsigned int address) const;
     bool                    hasCreatedObject(const std::string name) const;
     bool                    isObject5d(const unsigned int address) const;
@@ -192,20 +139,17 @@ public:
     template <typename T>
     bool                    isObjectOfType(const std::string name) const;
     Environment::Size       getTotalSize(void) const;
-    void                    addOwnership(const unsigned int owner,
+    void                    freeObject(const unsigned int address);
-                                         const unsigned int property);
+    void                    freeObject(const std::string name);
-    void                    addOwnership(const std::string owner,
-                                         const std::string property);
-    bool                    hasOwners(const unsigned int address) const;
-    bool                    hasOwners(const std::string name) const;
-    bool                    freeObject(const unsigned int address);
-    bool                    freeObject(const std::string name);
     void                    freeAll(void);
-    void                    printContent(void);
+    void                    protectObjects(const bool protect);
+    bool                    objectsProtected(void) const;
+    // print environment content
+    void                    printContent(void) const;
 private:
     // general
-    bool                                   dryRun_{false};
+    unsigned long int                      locVol_;
-    unsigned int                           traj_, locVol_;
+    bool                                   protect_{true};
     // grids
     std::vector<int>                       dim_;
     GridPt                                 grid4d_;
@@ -215,11 +159,6 @@ private:
     unsigned int                           nd_;
     // random number generator
     RngPt                                  rng4d_;
-    // module and related maps
-    std::vector<ModuleInfo>                module_;
-    std::map<std::string, unsigned int>    moduleAddress_;
-    // lattice store
-    std::map<unsigned int, LatticePt>      lattice_;
     // object store
     std::vector<ObjInfo>                   object_;
     std::map<std::string, unsigned int>    objectAddress_;
@@ -256,116 +195,85 @@ void Holder<T>::reset(T *pt)
 /******************************************************************************
  *                     Environment template implementation                    *
  ******************************************************************************/
-// module management ///////////////////////////////////////////////////////////
+// general memory management ///////////////////////////////////////////////////
-template <typename M>
+template <typename B, typename T, typename ... Ts>
-void Environment::createModule(const std::string name)
+void Environment::createDerivedObject(const std::string name,
+                                      const Environment::Storage storage,
+                                      const unsigned int Ls,
+                                      Ts && ... args)
 {
-    ModPt pt(new M(name));
+    if (!hasObject(name))
+    {
+        addObject(name);
+    }
     
-    pushModule(pt);
+    unsigned int address = getObjectAddress(name);
-}
-
-template <typename M>
-void Environment::createModule(const std::string name,
-                               const typename M::Par &par)
-{
-    ModPt pt(new M(name));
     
-    static_cast<M *>(pt.get())->setPar(par);
+    if (!object_[address].data or !objectsProtected())
-    pushModule(pt);
-}
-
-template <typename M>
-M * Environment::getModule(const unsigned int address) const
-{
-    if (auto *pt = dynamic_cast<M *>(getModule(address)))
     {
-        return pt;
+        MemoryStats memStats;
+    
+        if (!MemoryProfiler::stats)
+        {
+            MemoryProfiler::stats = &memStats;
+        }
+        size_t initMem           = MemoryProfiler::stats->currentlyAllocated;
+        object_[address].storage = storage;
+        object_[address].Ls      = Ls;
+        object_[address].data.reset(new Holder<B>(new T(std::forward<Ts>(args)...)));
+        object_[address].size    = MemoryProfiler::stats->maxAllocated - initMem;
+        object_[address].type    = &typeid(T);
+        if (MemoryProfiler::stats == &memStats)
+        {
+            MemoryProfiler::stats = nullptr;
+        }
     }
-    else
+    // object already exists, no error if it is a cache, error otherwise
+    else if ((object_[address].storage != Storage::cache) or 
+             (object_[address].storage != storage)        or
+             (object_[address].name    != name)           or
+             (object_[address].type    != &typeid(T)))
     {
-        HADRON_ERROR("module '" + module_[address].name
+        HADRON_ERROR(Definition, "object '" + name + "' already allocated");
-                     + "' does not have type " + typeid(M).name()
-                     + "(object type: " + getModuleType(address) + ")");
     }
 }
 
-template <typename M>
+template <typename T, typename ... Ts>
-M * Environment::getModule(const std::string name) const
+void Environment::createObject(const std::string name, 
+                               const Environment::Storage storage,
+                               const unsigned int Ls,
+                               Ts && ... args)
 {
-    return getModule<M>(getModuleAddress(name));
+    createDerivedObject<T, T>(name, storage, Ls, std::forward<Ts>(args)...);
-}
-
-template <typename T>
-unsigned int Environment::lattice4dSize(void) const
-{
-    return sizeof(typename T::vector_object)/getGrid()->Nsimd();
-}
-
-template <typename T>
-void Environment::registerLattice(const unsigned int address,
-                                  const unsigned int Ls)
-{
-    createGrid(Ls);
-    registerObject(address, Ls*lattice4dSize<T>(), Ls);
-}
-
-template <typename T>
-void Environment::registerLattice(const std::string name, const unsigned int Ls)
-{
-    createGrid(Ls);
-    registerObject(name, Ls*lattice4dSize<T>(), Ls);
-}
-
-template <typename T>
-void Environment::setObject(const unsigned int address, T *object)
-{
-    if (hasRegisteredObject(address))
-    {
-        object_[address].data.reset(new Holder<T>(object));
-        object_[address].type = &typeid(T);
-    }
-    else if (hasObject(address))
-    {
-        HADRON_ERROR("object with address " + std::to_string(address) +
-                     " exists but is not registered");
-    }
-    else
-    {
-        HADRON_ERROR("no object with address " + std::to_string(address));
-    }
-}
-
-template <typename T>
-void Environment::setObject(const std::string name, T *object)
-{
-    setObject(getObjectAddress(name), object);
 }
 
 template <typename T>
 T * Environment::getObject(const unsigned int address) const
 {
-    if (hasRegisteredObject(address))
+    if (hasObject(address))
     {
-        if (auto h = dynamic_cast<Holder<T> *>(object_[address].data.get()))
+        if (hasCreatedObject(address))
         {
-            return h->getPt();
+            if (auto h = dynamic_cast<Holder<T> *>(object_[address].data.get()))
+            {
+                return h->getPt();
+            }
+            else
+            {
+                HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
+                            " does not have type '" + typeName(&typeid(T)) +
+                            "' (has type '" + getObjectType(address) + "')");
+            }
         }
         else
         {
-            HADRON_ERROR("object with address " + std::to_string(address) +
+            HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
-                         " does not have type '" + typeName(&typeid(T)) +
+                         " is empty");
-                         "' (has type '" + getObjectType(address) + "')");
         }
     }
-    else if (hasObject(address))
-    {
-        HADRON_ERROR("object with address " + std::to_string(address) +
-                     " exists but is not registered");
-    }
     else
     {
-        HADRON_ERROR("no object with address " + std::to_string(address));
+        HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
     }
 }
 
@@ -375,26 +283,10 @@ T * Environment::getObject(const std::string name) const
     return getObject<T>(getObjectAddress(name));
 }
 
-template <typename T>
-T * Environment::createLattice(const unsigned int address)
-{
-    GridCartesian *g = getGrid(getObjectLs(address));
-    
-    setObject(address, new T(g));
-    
-    return getObject<T>(address);
-}
-
-template <typename T>
-T * Environment::createLattice(const std::string name)
-{
-    return createLattice<T>(getObjectAddress(name));
-}
-
 template <typename T>
 bool Environment::isObjectOfType(const unsigned int address) const
 {
-    if (hasRegisteredObject(address))
+    if (hasObject(address))
     {
         if (auto h = dynamic_cast<Holder<T> *>(object_[address].data.get()))
         {
@@ -405,14 +297,9 @@ bool Environment::isObjectOfType(const unsigned int address) const
             return false;
         }
     }
-    else if (hasObject(address))
-    {
-        HADRON_ERROR("object with address " + std::to_string(address) +
-                     " exists but is not registered");
-    }
     else
     {
-        HADRON_ERROR("no object with address " + std::to_string(address));
+        HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
     }
 }
 

extras/Hadrons/Exceptions.cc

@@ -0,0 +1,57 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Exceptions.cc
+
+Copyright (C) 2017
+
+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/Exceptions.hpp>
+
+#ifndef ERR_SUFF
+#define ERR_SUFF " (" + loc + ")"
+#endif
+
+#define CONST_EXC(name, init) \
+name::name(std::string msg, std::string loc)\
+:init\
+{}
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace Exceptions;
+
+// logic errors
+CONST_EXC(Logic, logic_error(msg + ERR_SUFF))
+CONST_EXC(Definition, Logic("definition error: " + msg, loc))
+CONST_EXC(Implementation, Logic("implementation error: " + msg, loc))
+CONST_EXC(Range, Logic("range error: " + msg, loc))
+CONST_EXC(Size, Logic("size error: " + msg, loc))
+// runtime errors
+CONST_EXC(Runtime, runtime_error(msg + ERR_SUFF))
+CONST_EXC(Argument, Runtime("argument error: " + msg, loc))
+CONST_EXC(Io, Runtime("IO error: " + msg, loc))
+CONST_EXC(Memory, Runtime("memory error: " + msg, loc))
+CONST_EXC(Parsing, Runtime("parsing error: " + msg, loc))
+CONST_EXC(Program, Runtime("program error: " + msg, loc))
+CONST_EXC(System, Runtime("system error: " + msg, loc))

extras/Hadrons/Exceptions.hpp

@@ -0,0 +1,72 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Exceptions.hpp
+
+Copyright (C) 2017
+
+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_Exceptions_hpp_
+#define Hadrons_Exceptions_hpp_
+
+#include <stdexcept>
+#ifndef Hadrons_Global_hpp_
+#include <Grid/Hadrons/Global.hpp>
+#endif
+
+#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) \
+class name: public base\
+{\
+public:\
+    name(std::string msg, std::string loc);\
+}
+
+BEGIN_HADRONS_NAMESPACE
+
+namespace Exceptions
+{
+    // logic errors
+    DECL_EXC(Logic, std::logic_error);
+    DECL_EXC(Definition, Logic);
+    DECL_EXC(Implementation, Logic);
+    DECL_EXC(Range, Logic);
+    DECL_EXC(Size, Logic);
+    // runtime errors
+    DECL_EXC(Runtime, std::runtime_error);
+    DECL_EXC(Argument, Runtime);
+    DECL_EXC(Io, Runtime);
+    DECL_EXC(Memory, Runtime);
+    DECL_EXC(Parsing, Runtime);
+    DECL_EXC(Program, Runtime);
+    DECL_EXC(System, Runtime);
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Exceptions_hpp_

extras/Hadrons/Factory.hpp

@@ -95,7 +95,7 @@ std::unique_ptr<T> Factory<T>::create(const std::string type,
     }
     catch (std::out_of_range &)
     {
-        HADRON_ERROR("object of type '" + type + "' unknown");
+        HADRON_ERROR(Argument, "object of type '" + type + "' unknown");
     }
     
     return func(name);

extras/Hadrons/GeneticScheduler.hpp

@@ -38,13 +38,13 @@ BEGIN_HADRONS_NAMESPACE
 /******************************************************************************
  *                   Scheduler based on a genetic algorithm                   *
  ******************************************************************************/
-template <typename T>
+template <typename V, typename T>
 class GeneticScheduler
 {
 public:
-    typedef std::vector<T>                   Gene;
+    typedef std::vector<T>                 Gene;
-    typedef std::pair<Gene *, Gene *>        GenePair;
+    typedef std::pair<Gene *, Gene *>      GenePair;
-    typedef std::function<int(const Gene &)> ObjFunc;
+    typedef std::function<V(const Gene &)> ObjFunc;
     struct Parameters
     {
         double       mutationRate;
@@ -65,7 +65,7 @@ public:
     void benchmarkCrossover(const unsigned int nIt);
     // print population
     friend std::ostream & operator<<(std::ostream &out,
-                                     const GeneticScheduler<T> &s)
+                                     const GeneticScheduler<V, T> &s)
     {
         out << "[";
         for (auto &p: s.population_)
@@ -87,19 +87,19 @@ private:
     void     mutation(Gene &m, const Gene &c);
     
 private:
-    Graph<T>                 &graph_;
+    Graph<T>               &graph_;
-    const ObjFunc            &func_;
+    const ObjFunc          &func_;
-    const Parameters         par_;
+    const Parameters       par_;
-    std::multimap<int, Gene> population_;
+    std::multimap<V, Gene> population_;
-    std::mt19937             gen_;
+    std::mt19937           gen_;
 };
 
 /******************************************************************************
  *                       template implementation                              *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
-template <typename T>
+template <typename V, typename T>
-GeneticScheduler<T>::GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
+GeneticScheduler<V, T>::GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
                                       const Parameters &par)
 : graph_(graph)
 , func_(func)
@@ -109,22 +109,22 @@ GeneticScheduler<T>::GeneticScheduler(Graph<T> &graph, const ObjFunc &func,
 }
 
 // access //////////////////////////////////////////////////////////////////////
-template <typename T>
+template <typename V, typename T>
-const typename GeneticScheduler<T>::Gene &
+const typename GeneticScheduler<V, T>::Gene &
-GeneticScheduler<T>::getMinSchedule(void)
+GeneticScheduler<V, T>::getMinSchedule(void)
 {
     return population_.begin()->second;
 }
 
-template <typename T>
+template <typename V, typename T>
-int GeneticScheduler<T>::getMinValue(void)
+int GeneticScheduler<V, T>::getMinValue(void)
 {
     return population_.begin()->first;
 }
 
 // breed a new generation //////////////////////////////////////////////////////
-template <typename T>
+template <typename V, typename T>
-void GeneticScheduler<T>::nextGeneration(void)
+void GeneticScheduler<V, T>::nextGeneration(void)
 {
     // random initialization of the population if necessary
     if (population_.size() != par_.popSize)
@@ -158,8 +158,8 @@ void GeneticScheduler<T>::nextGeneration(void)
 }
 
 // evolution steps /////////////////////////////////////////////////////////////
-template <typename T>
+template <typename V, typename T>
-void GeneticScheduler<T>::initPopulation(void)
+void GeneticScheduler<V, T>::initPopulation(void)
 {
     population_.clear();
     for (unsigned int i = 0; i < par_.popSize; ++i)
@@ -170,8 +170,8 @@ void GeneticScheduler<T>::initPopulation(void)
     }
 }
 
-template <typename T>
+template <typename V, typename T>
-void GeneticScheduler<T>::doCrossover(void)
+void GeneticScheduler<V, T>::doCrossover(void)
 {
     auto p = selectPair();
     Gene &p1 = *(p.first), &p2 = *(p.second);
@@ -185,8 +185,8 @@ void GeneticScheduler<T>::doCrossover(void)
     }
 }
 
-template <typename T>
+template <typename V, typename T>
-void GeneticScheduler<T>::doMutation(void)
+void GeneticScheduler<V, T>::doMutation(void)
 {
     std::uniform_real_distribution<double>      mdis(0., 1.);
     std::uniform_int_distribution<unsigned int> pdis(0, population_.size() - 1);
@@ -206,40 +206,35 @@ void GeneticScheduler<T>::doMutation(void)
 }
 
 // genetic operators ///////////////////////////////////////////////////////////
-template <typename T>
+template <typename V, typename T>
-typename GeneticScheduler<T>::GenePair GeneticScheduler<T>::selectPair(void)
+typename GeneticScheduler<V, T>::GenePair GeneticScheduler<V, T>::selectPair(void)
 {
     std::vector<double> prob;
     unsigned int        ind;
     Gene                *p1, *p2;
+    const double        max = population_.rbegin()->first;
     
+
     for (auto &c: population_)
     {
-        prob.push_back(1./c.first);
+        prob.push_back(std::exp((c.first-1.)/max));
-    }
+    }        
-    do
+    std::discrete_distribution<unsigned int> dis1(prob.begin(), prob.end());
-    {
+    auto rIt = population_.begin();
-        double probCpy;
+    ind = dis1(gen_);
-        
+    std::advance(rIt, ind);
-        std::discrete_distribution<unsigned int> dis1(prob.begin(), prob.end());
+    p1 = &(rIt->second);
-        auto rIt = population_.begin();
+    prob[ind] = 0.;
-        ind = dis1(gen_);
+    std::discrete_distribution<unsigned int> dis2(prob.begin(), prob.end());
-        std::advance(rIt, ind);
+    rIt = population_.begin();
-        p1 = &(rIt->second);
+    std::advance(rIt, dis2(gen_));
-        probCpy   = prob[ind];
+    p2 = &(rIt->second);
-        prob[ind] = 0.;
-        std::discrete_distribution<unsigned int> dis2(prob.begin(), prob.end());
-        rIt = population_.begin();
-        std::advance(rIt, dis2(gen_));
-        p2 = &(rIt->second);
-        prob[ind] = probCpy;
-    } while (p1 == p2);
     
     return std::make_pair(p1, p2);
 }
 
-template <typename T>
+template <typename V, typename T>
-void GeneticScheduler<T>::crossover(Gene &c1, Gene &c2, const Gene &p1,
+void GeneticScheduler<V, T>::crossover(Gene &c1, Gene &c2, const Gene &p1,
                                     const Gene &p2)
 {
     Gene                                        buf;
@@ -273,8 +268,8 @@ void GeneticScheduler<T>::crossover(Gene &c1, Gene &c2, const Gene &p1,
     }
 }
 
-template <typename T>
+template <typename V, typename T>
-void GeneticScheduler<T>::mutation(Gene &m, const Gene &c)
+void GeneticScheduler<V, T>::mutation(Gene &m, const Gene &c)
 {
     Gene                                        buf;
     std::uniform_int_distribution<unsigned int> dis(0, c.size() - 1);
@@ -303,8 +298,8 @@ void GeneticScheduler<T>::mutation(Gene &m, const Gene &c)
     }
 }
 
-template <typename T>
+template <typename V, typename T>
-void GeneticScheduler<T>::benchmarkCrossover(const unsigned int nIt)
+void GeneticScheduler<V, T>::benchmarkCrossover(const unsigned int nIt)
 {
     Gene   p1, p2, c1, c2;
     double neg = 0., eq = 0., pos = 0., total;

extras/Hadrons/Global.cc

@@ -39,33 +39,6 @@ HadronsLogger Hadrons::HadronsLogMessage(1,"Message");
 HadronsLogger Hadrons::HadronsLogIterative(1,"Iterative");
 HadronsLogger Hadrons::HadronsLogDebug(1,"Debug");
 
-// pretty size formatting //////////////////////////////////////////////////////
-std::string Hadrons::sizeString(long unsigned int bytes)
-
-{
-    constexpr unsigned int bufSize = 256;
-    const char             *suffixes[7] = {"", "K", "M", "G", "T", "P", "E"};
-    char                   buf[256];
-    long unsigned int      s     = 0;
-    double                 count = bytes;
-    
-    while (count >= 1024 && s < 7)
-    {
-        s++;
-        count /= 1024;
-    }
-    if (count - floor(count) == 0.0)
-    {
-        snprintf(buf, bufSize, "%d %sB", (int)count, suffixes[s]);
-    }
-    else
-    {
-        snprintf(buf, bufSize, "%.1f %sB", count, suffixes[s]);
-    }
-    
-    return std::string(buf);
-}
-
 // type utilities //////////////////////////////////////////////////////////////
 constexpr unsigned int maxNameSize = 1024u;
 

extras/Hadrons/Global.hpp

@@ -35,6 +35,10 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Grid.h>
 #include <cxxabi.h>
 
+#ifndef SITE_SIZE_TYPE
+#define SITE_SIZE_TYPE size_t
+#endif
+
 #define BEGIN_HADRONS_NAMESPACE \
 namespace Grid {\
 using namespace QCD;\
@@ -62,12 +66,11 @@ BEGIN_HADRONS_NAMESPACE
 
 // type aliases
 #define FERM_TYPE_ALIASES(FImpl, suffix)\
-typedef FermionOperator<FImpl>                       FMat##suffix;             \
+typedef FermionOperator<FImpl>                        FMat##suffix;            \
-typedef typename FImpl::FermionField                 FermionField##suffix;     \
+typedef typename FImpl::FermionField                  FermionField##suffix;    \
-typedef typename FImpl::PropagatorField              PropagatorField##suffix;  \
+typedef typename FImpl::PropagatorField               PropagatorField##suffix; \
-typedef typename FImpl::SitePropagator               SitePropagator##suffix;   \
+typedef typename FImpl::SitePropagator::scalar_object SitePropagator##suffix;  \
-typedef std::vector<typename FImpl::SitePropagator::scalar_object>             \
+typedef std::vector<SitePropagator##suffix>           SlicedPropagator##suffix;
-                                                     SlicedPropagator##suffix;
 
 #define GAUGE_TYPE_ALIASES(FImpl, suffix)\
 typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;
@@ -97,11 +100,6 @@ public:
 };
 
 #define LOG(channel) std::cout << HadronsLog##channel
-#define HADRON_ERROR(msg)\
-LOG(Error) << msg << " (" << __FUNCTION__ << " at " << __FILE__ << ":"\
-           << __LINE__ << ")" << std::endl;\
-abort();
-
 #define DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl;
 
 extern HadronsLogger HadronsLogError;
@@ -135,9 +133,6 @@ public:\
 private:\
     name(void) = default;
 
-// pretty size formating
-std::string sizeString(long unsigned int bytes);
-
 // type utilities
 template <typename T>
 const std::type_info * typeIdPt(const T &x)
@@ -176,4 +171,6 @@ typedef XmlWriter CorrWriter;
 
 END_HADRONS_NAMESPACE
 
+#include <Grid/Hadrons/Exceptions.hpp>
+
 #endif // Hadrons_Global_hpp_

extras/Hadrons/Graph.hpp

@@ -185,7 +185,7 @@ void Graph<T>::removeVertex(const T &value)
     }
     else
     {
-        HADRON_ERROR("vertex " << value << " does not exists");
+        HADRON_ERROR(Range, "vertex does not exists");
     }
 
     // remove all edges containing the vertex
@@ -214,7 +214,7 @@ void Graph<T>::removeEdge(const Edge &e)
     }
     else
     {
-        HADRON_ERROR("edge "  << e << " does not exists");
+        HADRON_ERROR(Range, "edge does not exists");
     }
 }
 
@@ -260,7 +260,7 @@ void Graph<T>::mark(const T &value, const bool doMark)
     }
     else
     {
-        HADRON_ERROR("vertex " << value << " does not exists");
+        HADRON_ERROR(Range, "vertex does not exists");
     }
 }
 
@@ -298,7 +298,7 @@ bool Graph<T>::isMarked(const T &value) const
     }
     else
     {
-        HADRON_ERROR("vertex " << value << " does not exists");
+        HADRON_ERROR(Range, "vertex does not exists");
         
         return false;
     }
@@ -430,7 +430,7 @@ std::vector<T> Graph<T>::getAdjacentVertices(const T &value) const
     {
         return ((e.first == value) or (e.second == value));
     };
-    auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred);
+    auto eIt = std::find_if(edgeSet_.begin(), edgeSet_.end(), pred);
     
     while (eIt != edgeSet_.end())
     {
@@ -442,7 +442,7 @@ std::vector<T> Graph<T>::getAdjacentVertices(const T &value) const
         {
             adjacentVertex.push_back((*eIt).first);
         }
-        eIt = find_if(++eIt, edgeSet_.end(), pred);
+        eIt = std::find_if(++eIt, edgeSet_.end(), pred);
     }
     
     return adjacentVertex;
@@ -458,12 +458,12 @@ std::vector<T> Graph<T>::getChildren(const T &value) const
     {
         return (e.first == value);
     };
-    auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred);
+    auto eIt = std::find_if(edgeSet_.begin(), edgeSet_.end(), pred);
     
     while (eIt != edgeSet_.end())
     {
         child.push_back((*eIt).second);
-        eIt = find_if(++eIt, edgeSet_.end(), pred);
+        eIt = std::find_if(++eIt, edgeSet_.end(), pred);
     }
     
     return child;
@@ -479,12 +479,12 @@ std::vector<T> Graph<T>::getParents(const T &value) const
     {
         return (e.second == value);
     };
-    auto eIt = find_if(edgeSet_.begin(), edgeSet_.end(), pred);
+    auto eIt = std::find_if(edgeSet_.begin(), edgeSet_.end(), pred);
     
     while (eIt != edgeSet_.end())
     {
         parent.push_back((*eIt).first);
-        eIt = find_if(++eIt, edgeSet_.end(), pred);
+        eIt = std::find_if(++eIt, edgeSet_.end(), pred);
     }
     
     return parent;
@@ -544,7 +544,7 @@ std::vector<T> Graph<T>::topoSort(void)
     {
         if (tmpMarked.at(v))
         {
-            HADRON_ERROR("cannot topologically sort a cyclic graph");
+            HADRON_ERROR(Range, "cannot topologically sort a cyclic graph");
         }
         if (!isMarked(v))
         {
@@ -603,7 +603,7 @@ std::vector<T> Graph<T>::topoSort(Gen &gen)
     {
         if (tmpMarked.at(v))
         {
-            HADRON_ERROR("cannot topologically sort a cyclic graph");
+            HADRON_ERROR(Range, "cannot topologically sort a cyclic graph");
         }
         if (!isMarked(v))
         {

extras/Hadrons/Makefile.am

@@ -7,20 +7,24 @@ libHadrons_a_SOURCES = \
     $(modules_cc)      \
     Application.cc     \
     Environment.cc     \
+	Exceptions.cc      \
     Global.cc          \
-    Module.cc
+    Module.cc		   \
+	VirtualMachine.cc
 libHadrons_adir = $(pkgincludedir)/Hadrons
 nobase_libHadrons_a_HEADERS = \
 	$(modules_hpp)            \
 	Application.hpp           \
 	Environment.hpp           \
+	Exceptions.hpp            \
 	Factory.hpp               \
 	GeneticScheduler.hpp      \
 	Global.hpp                \
 	Graph.hpp                 \
 	Module.hpp                \
 	Modules.hpp               \
-	ModuleFactory.hpp
+	ModuleFactory.hpp         \
+	VirtualMachine.hpp
 
 HadronsXmlRun_SOURCES = HadronsXmlRun.cc
 HadronsXmlRun_LDADD   = libHadrons.a -lGrid

extras/Hadrons/Module.cc

@@ -39,7 +39,6 @@ using namespace Hadrons;
 // constructor /////////////////////////////////////////////////////////////////
 ModuleBase::ModuleBase(const std::string name)
 : name_(name)
-, env_(Environment::getInstance())
 {}
 
 // access //////////////////////////////////////////////////////////////////////
@@ -48,15 +47,10 @@ std::string ModuleBase::getName(void) const
     return name_;
 }
 
-Environment & ModuleBase::env(void) const
-{
-    return env_;
-}
-
 // get factory registration name if available
 std::string ModuleBase::getRegisteredName(void)
 {
-    HADRON_ERROR("module '" + getName() + "' has a type not registered"
+    HADRON_ERROR(Definition, "module '" + getName() + "' has no registered type"
                  + " in the factory");
 }
 
@@ -64,8 +58,5 @@ std::string ModuleBase::getRegisteredName(void)
 void ModuleBase::operator()(void)
 {
     setup();
-    if (!env().isDryRun())
+    execute();
-    {
-        execute();
-    }
 }

extras/Hadrons/Module.hpp

@@ -31,7 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define Hadrons_Module_hpp_
 
 #include <Grid/Hadrons/Global.hpp>
-#include <Grid/Hadrons/Environment.hpp>
+#include <Grid/Hadrons/VirtualMachine.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -87,6 +87,56 @@ public:\
 static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance;
 
 #define ARG(...) __VA_ARGS__
+#define MACRO_REDIRECT(arg1, arg2, arg3, macro, ...) macro
+
+#define envGet(type, name)\
+*env().template getObject<type>(name)
+
+#define envGetTmp(type, var)\
+type &var = *env().template getObject<type>(getName() + "_tmp_" + #var)
+
+#define envHasType(type, name)\
+env().template isObjectOfType<type>(name)
+
+#define envCreate(type, name, Ls, ...)\
+env().template createObject<type>(name, Environment::Storage::object, Ls, __VA_ARGS__)
+
+#define envCreateDerived(base, type, name, Ls, ...)\
+env().template createDerivedObject<base, type>(name, Environment::Storage::object, Ls, __VA_ARGS__)
+
+#define envCreateLat4(type, name)\
+envCreate(type, name, 1, env().getGrid())
+
+#define envCreateLat5(type, name, Ls)\
+envCreate(type, name, Ls, env().getGrid(Ls))
+
+#define envCreateLat(...)\
+MACRO_REDIRECT(__VA_ARGS__, envCreateLat5, envCreateLat4)(__VA_ARGS__)
+
+#define envCache(type, name, Ls, ...)\
+env().template createObject<type>(name, Environment::Storage::cache, Ls, __VA_ARGS__)
+
+#define envCacheLat4(type, name)\
+envCache(type, name, 1, env().getGrid())
+
+#define envCacheLat5(type, name, Ls)\
+envCache(type, name, Ls, env().getGrid(Ls))
+
+#define envCacheLat(...)\
+MACRO_REDIRECT(__VA_ARGS__, envCacheLat5, envCacheLat4)(__VA_ARGS__)
+
+#define envTmp(type, name, Ls, ...)\
+env().template createObject<type>(getName() + "_tmp_" + name,         \
+                                  Environment::Storage::temporary, Ls, __VA_ARGS__)
+
+#define envTmpLat4(type, name)\
+envTmp(type, name, 1, env().getGrid())
+
+#define envTmpLat5(type, name, Ls)\
+envTmp(type, name, Ls, env().getGrid(Ls))
+
+#define envTmpLat(...)\
+MACRO_REDIRECT(__VA_ARGS__, envTmpLat5, envTmpLat4)(__VA_ARGS__)
 
 /******************************************************************************
  *                            Module class                                    *
@@ -101,23 +151,30 @@ public:
     virtual ~ModuleBase(void) = default;
     // access
     std::string getName(void) const;
-    Environment &env(void) const;
     // get factory registration name if available
     virtual std::string getRegisteredName(void);
     // dependencies/products
     virtual std::vector<std::string> getInput(void) = 0;
+    virtual std::vector<std::string> getReference(void)
+    {
+        return std::vector<std::string>(0);
+    };
     virtual std::vector<std::string> getOutput(void) = 0;
     // parse parameters
     virtual void parseParameters(XmlReader &reader, const std::string name) = 0;
     virtual void saveParameters(XmlWriter &writer, const std::string name) = 0;
     // setup
     virtual void setup(void) {};
+    virtual void execute(void) = 0;
     // execution
     void operator()(void);
-    virtual void execute(void) = 0;
+protected:
+    // environment shortcut
+    DEFINE_ENV_ALIAS;
+    // virtual machine shortcut
+    DEFINE_VM_ALIAS;
 private:
     std::string name_;
-    Environment &env_;
 };
 
 // derived class, templating the parameter class

extras/Hadrons/Modules.hpp

@@ -1,25 +1,60 @@
-#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
+/*************************************************************************************
-#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+Copyright (C) 2017
+
+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/Modules/MContraction/Baryon.hpp>
-#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
-#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
 #include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
-#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
 #include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
-#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
-#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
-#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
-#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
-#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
-#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
-#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
-#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
-#include <Grid/Hadrons/Modules/MSink/Point.hpp>
-#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
-#include <Grid/Hadrons/Modules/MSource/Point.hpp>
 #include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MSource/Point.hpp>
 #include <Grid/Hadrons/Modules/MSource/Wall.hpp>
 #include <Grid/Hadrons/Modules/MSource/Z2.hpp>
+#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/RBPrecCG.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Load.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
+#include <Grid/Hadrons/Modules/MGauge/Random.hpp>
+#include <Grid/Hadrons/Modules/MGauge/StochEm.hpp>
+#include <Grid/Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
+#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
+#include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
+#include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
+#include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>
+#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
+#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>

extras/Hadrons/Modules/MAction/DWF.hpp

@@ -65,6 +65,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
@@ -103,35 +104,29 @@ std::vector<std::string> TDWF<FImpl>::getOutput(void)
 template <typename FImpl>
 void TDWF<FImpl>::setup(void)
 {
-    unsigned int size;
+    LOG(Message) << "Setting up domain wall fermion matrix with m= "
-    
+                 << par().mass << ", M5= " << par().M5 << " and Ls= "
-    size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
+                 << par().Ls << " using gauge field '" << par().gauge << "'"
-    env().registerObject(getName(), size, par().Ls);
+                 << 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);
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename DomainWallFermion<FImpl>::ImplParams implParams(boundary);
+    envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
+                     grb5, g4, grb4, par().mass, par().M5, implParams);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TDWF<FImpl>::execute(void)
-{
+{}
-    LOG(Message) << "Setting up domain wall fermion matrix with m= "
-                 << par().mass << ", M5= " << par().M5 << " and Ls= "
-                 << par().Ls << " using gauge field '" << par().gauge << "'"
-                 << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary 
-                 << std::endl;
-    env().createGrid(par().Ls);
-    auto &U      = *env().template getObject<LatticeGaugeField>(par().gauge);
-    auto &g4     = *env().getGrid();
-    auto &grb4   = *env().getRbGrid();
-    auto &g5     = *env().getGrid(par().Ls);
-    auto &grb5   = *env().getRbGrid(par().Ls);
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
-    typename DomainWallFermion<FImpl>::ImplParams implParams(boundary);
-    FMat *fMatPt = new DomainWallFermion<FImpl>(U, g5, grb5, g4, grb4,
-                                                par().mass, par().M5,
-                                                implParams);
-    env().setObject(getName(), fMatPt);
-}
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MAction/Wilson.hpp

@@ -63,6 +63,7 @@ public:
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
@@ -101,29 +102,24 @@ std::vector<std::string> TWilson<FImpl>::getOutput(void)
 template <typename FImpl>
 void TWilson<FImpl>::setup(void)
 {
-    unsigned int size;
+    LOG(Message) << "Setting up TWilson fermion matrix with m= " << par().mass
-    
+                 << " using gauge field '" << par().gauge << "'" << std::endl;
-    size = 2*env().template lattice4dSize<typename FImpl::DoubledGaugeField>();
+    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-    env().registerObject(getName(), size);
+                 << std::endl;
+                 
+    auto &U      = envGet(LatticeGaugeField, par().gauge);
+    auto &grid   = *env().getGrid();
+    auto &gridRb = *env().getRbGrid();
+    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
+    typename WilsonFermion<FImpl>::ImplParams implParams(boundary);
+    envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
+                     par().mass, implParams);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TWilson<FImpl>::execute()
-{
+{}
-    LOG(Message) << "Setting up TWilson fermion matrix with m= " << par().mass
-                 << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary 
-                 << std::endl;
-    auto &U      = *env().template getObject<LatticeGaugeField>(par().gauge);
-    auto &grid   = *env().getGrid();
-    auto &gridRb = *env().getRbGrid();
-    std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
-    typename WilsonFermion<FImpl>::ImplParams implParams(boundary);
-    FMat *fMatPt = new WilsonFermion<FImpl>(U, grid, gridRb, par().mass,
-                                            implParams);
-    env().setObject(getName(), fMatPt);
-}
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MContraction/Baryon.hpp

@@ -72,6 +72,9 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
+    // setup
+    virtual void setup(void);
     // execution
     virtual void execute(void);
 };
@@ -99,11 +102,18 @@ std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getInput(void)
 template <typename FImpl1, typename FImpl2, typename FImpl3>
 std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {};
     
     return out;
 }
 
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+void TBaryon<FImpl1, FImpl2, FImpl3>::setup(void)
+{
+    envTmpLat(LatticeComplex, "c");
+}
+
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl1, typename FImpl2, typename FImpl3>
 void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
@@ -112,12 +122,12 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
                  << " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
                  << par().q3 << "'" << std::endl;
     
-    CorrWriter             writer(par().output);
+    CorrWriter writer(par().output);
-    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
+    auto       &q1 = envGet(PropagatorField1, par().q1);
-    PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
+    auto       &q2 = envGet(PropagatorField2, par().q2);
-    PropagatorField3      &q3 = *env().template getObject<PropagatorField3>(par().q2);
+    auto       &q3 = envGet(PropagatorField3, par().q2);
-    LatticeComplex        c(env().getGrid());
+    envGetTmp(LatticeComplex, c);
-    Result                result;
+    Result     result;
     
     // FIXME: do contractions
     

extras/Hadrons/Modules/MContraction/DiscLoop.hpp

@@ -68,6 +68,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
@@ -97,7 +98,7 @@ std::vector<std::string> TDiscLoop<FImpl>::getInput(void)
 template <typename FImpl>
 std::vector<std::string> TDiscLoop<FImpl>::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {};
     
     return out;
 }
@@ -106,7 +107,7 @@ std::vector<std::string> TDiscLoop<FImpl>::getOutput(void)
 template <typename FImpl>
 void TDiscLoop<FImpl>::setup(void)
 {
-    
+    envTmpLat(LatticeComplex, "c");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -118,12 +119,12 @@ void TDiscLoop<FImpl>::execute(void)
                  << " insertion." << std::endl;
 
     CorrWriter            writer(par().output);
-    PropagatorField       &q_loop = *env().template getObject<PropagatorField>(par().q_loop);
+    auto                  &q_loop = envGet(PropagatorField, par().q_loop);
-    LatticeComplex        c(env().getGrid());
     Gamma                 gamma(par().gamma);
     std::vector<TComplex> buf;
     Result                result;
 
+    envGetTmp(LatticeComplex, c);
     c = trace(gamma*q_loop);
     sliceSum(c, buf, Tp);
 

extras/Hadrons/Modules/MContraction/Gamma3pt.hpp

@@ -51,6 +51,14 @@ BEGIN_HADRONS_NAMESPACE
  *                   q1
  *
  *      trace(g5*q1*adj(q2)*g5*gamma*q3)
+ * 
+ *  options:
+ *   - q1: sink smeared propagator, source at i
+ *   - q2: propagator, source at i
+ *   - q3: propagator, source at f
+ *   - gamma: gamma matrix to insert
+ *   - tSnk: sink position for propagator q1.
+ *
  */
 
 /******************************************************************************
@@ -66,6 +74,7 @@ public:
                                     std::string,    q2,
                                     std::string,    q3,
                                     Gamma::Algebra, gamma,
+                                    unsigned int,   tSnk,
                                     std::string,    output);
 };
 
@@ -90,6 +99,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
@@ -119,7 +129,7 @@ std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getInput(void)
 template <typename FImpl1, typename FImpl2, typename FImpl3>
 std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {};
     
     return out;
 }
@@ -128,7 +138,7 @@ std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getOutput(void)
 template <typename FImpl1, typename FImpl2, typename FImpl3>
 void TGamma3pt<FImpl1, FImpl2, FImpl3>::setup(void)
 {
-    
+    envTmpLat(LatticeComplex, "c");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -140,17 +150,22 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
                  << par().q3 << "', with " << par().gamma << " insertion." 
                  << std::endl;
 
+    // Initialise variables. q2 and q3 are normal propagators, q1 may be 
+    // sink smeared.
     CorrWriter            writer(par().output);
-    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
+    auto                  &q1 = envGet(SlicedPropagator1, par().q1);
-    PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
+    auto                  &q2 = envGet(PropagatorField2, par().q2);
-    PropagatorField3      &q3 = *env().template getObject<PropagatorField3>(par().q3);
+    auto                  &q3 = envGet(PropagatorField2, par().q3);
-    LatticeComplex        c(env().getGrid());
     Gamma                 g5(Gamma::Algebra::Gamma5);
     Gamma                 gamma(par().gamma);
     std::vector<TComplex> buf;
     Result                result;
-
+    
-    c = trace(g5*q1*adj(q2)*(g5*gamma)*q3);
+    // Extract relevant timeslice of sinked propagator q1, then contract &
+    // sum over all spacial positions of gamma insertion.
+    SitePropagator1 q1Snk = q1[par().tSnk];
+    envGetTmp(LatticeComplex, c);
+    c = trace(g5*q1Snk*adj(q2)*(g5*gamma)*q3);
     sliceSum(c, buf, Tp);
 
     result.gamma = par().gamma;

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -51,8 +51,7 @@ BEGIN_HADRONS_NAMESPACE
            in a sequence (e.g. "<Gamma5 Gamma5><Gamma5 GammaT>").
 
            Special values: "all" - perform all possible contractions.
- - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0."),
+ - sink: module to compute the sink to use in contraction (string).
-        given as multiples of (2*pi) / L.
 */
 
 /******************************************************************************
@@ -98,6 +97,9 @@ public:
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
     virtual void parseGammaString(std::vector<GammaPair> &gammaList);
+protected:
+    // execution
+    virtual void setup(void);
     // execution
     virtual void execute(void);
 };
@@ -125,7 +127,7 @@ std::vector<std::string> TMeson<FImpl1, FImpl2>::getInput(void)
 template <typename FImpl1, typename FImpl2>
 std::vector<std::string> TMeson<FImpl1, FImpl2>::getOutput(void)
 {
-    std::vector<std::string> output = {getName()};
+    std::vector<std::string> output = {};
     
     return output;
 }
@@ -151,9 +153,15 @@ void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
     {
         // Parse individual contractions from input string.
         gammaList = strToVec<GammaPair>(par().gammas);
-    }
+    } 
 }
 
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2>
+void TMeson<FImpl1, FImpl2>::setup(void)
+{
+    envTmpLat(LatticeComplex, "c");
+}
 
 // execution ///////////////////////////////////////////////////////////////////
 #define mesonConnected(q1, q2, gSnk, gSrc) \
@@ -181,11 +189,11 @@ void TMeson<FImpl1, FImpl2>::execute(void)
         result[i].gamma_src = gammaList[i].second;
         result[i].corr.resize(nt);
     }
-    if (env().template isObjectOfType<SlicedPropagator1>(par().q1) and
+    if (envHasType(SlicedPropagator1, par().q1) and
-        env().template isObjectOfType<SlicedPropagator2>(par().q2))
+        envHasType(SlicedPropagator2, par().q2))
     {
-        SlicedPropagator1 &q1 = *env().template getObject<SlicedPropagator1>(par().q1);
+        auto &q1 = envGet(SlicedPropagator1, par().q1);
-        SlicedPropagator2 &q2 = *env().template getObject<SlicedPropagator2>(par().q2);
+        auto &q2 = envGet(SlicedPropagator2, par().q2);
         
         LOG(Message) << "(propagator already sinked)" << std::endl;
         for (unsigned int i = 0; i < result.size(); ++i)
@@ -201,10 +209,10 @@ void TMeson<FImpl1, FImpl2>::execute(void)
     }
     else
     {
-        PropagatorField1 &q1   = *env().template getObject<PropagatorField1>(par().q1);
+        auto &q1 = envGet(PropagatorField1, par().q1);
-        PropagatorField2 &q2   = *env().template getObject<PropagatorField2>(par().q2);
+        auto &q2 = envGet(PropagatorField2, par().q2);
-        LatticeComplex   c(env().getGrid());
         
+        envGetTmp(LatticeComplex, c);
         LOG(Message) << "(using sink '" << par().sink << "')" << std::endl;
         for (unsigned int i = 0; i < result.size(); ++i)
         {
@@ -212,18 +220,17 @@ void TMeson<FImpl1, FImpl2>::execute(void)
             Gamma       gSrc(gammaList[i].second);
             std::string ns;
                 
-            ns = env().getModuleNamespace(env().getObjectModule(par().sink));
+            ns = vm().getModuleNamespace(env().getObjectModule(par().sink));
             if (ns == "MSource")
             {
-                PropagatorField1 &sink =
+                PropagatorField1 &sink = envGet(PropagatorField1, par().sink);
-                    *env().template getObject<PropagatorField1>(par().sink);
                 
                 c = trace(mesonConnected(q1, q2, gSnk, gSrc)*sink);
                 sliceSum(c, buf, Tp);
             }
             else if (ns == "MSink")
             {
-                SinkFnScalar &sink = *env().template getObject<SinkFnScalar>(par().sink);
+                SinkFnScalar &sink = envGet(SinkFnScalar, par().sink);
                 
                 c   = trace(mesonConnected(q1, q2, gSnk, gSrc));
                 buf = sink(c);

extras/Hadrons/Modules/MContraction/WardIdentity.hpp

@@ -0,0 +1,223 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WardIdentity.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_MContraction_WardIdentity_hpp_
+#define Hadrons_MContraction_WardIdentity_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+  Ward Identity contractions
+ -----------------------------
+ 
+ * options:
+ - q:          propagator, 5D if available (string)
+ - action:     action module used for propagator solution (string)
+ - mass:       mass of quark (double)
+ - test_axial: whether or not to test PCAC relation.
+*/
+
+/******************************************************************************
+ *                              WardIdentity                                  *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class WardIdentityPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WardIdentityPar,
+                                    std::string, q,
+                                    std::string, action,
+                                    double,      mass,
+                                    bool,        test_axial);
+};
+
+template <typename FImpl>
+class TWardIdentity: public Module<WardIdentityPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TWardIdentity(const std::string name);
+    // destructor
+    virtual ~TWardIdentity(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+protected:
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+private:
+    unsigned int Ls_;
+};
+
+MODULE_REGISTER_NS(WardIdentity, TWardIdentity<FIMPL>, MContraction);
+
+/******************************************************************************
+ *                     TWardIdentity implementation                           *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWardIdentity<FImpl>::TWardIdentity(const std::string name)
+: Module<WardIdentityPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWardIdentity<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q, par().action};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWardIdentity<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWardIdentity<FImpl>::setup(void)
+{
+    Ls_ = env().getObjectLs(par().q);
+    if (Ls_ != env().getObjectLs(par().action))
+    {
+        HADRON_ERROR(Size, "Ls mismatch between quark action and propagator");
+    }
+    envTmpLat(PropagatorField, "tmp");
+    envTmpLat(PropagatorField, "vector_WI");
+    if (par().test_axial)
+    {
+        envTmpLat(PropagatorField, "psi");
+        envTmpLat(LatticeComplex,  "PP");
+        envTmpLat(LatticeComplex,  "axial_defect");
+        envTmpLat(LatticeComplex,  "PJ5q");
+    }
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWardIdentity<FImpl>::execute(void)
+{
+    LOG(Message) << "Performing Ward Identity checks for quark '" << par().q
+                 << "'." << std::endl;
+
+    auto  &q   = envGet(PropagatorField, par().q);
+    auto  &act = envGet(FMat, par().action);
+    Gamma g5(Gamma::Algebra::Gamma5);
+
+    // Compute D_mu V_mu, D here is backward derivative.
+    envGetTmp(PropagatorField, tmp);
+    envGetTmp(PropagatorField, vector_WI);
+    vector_WI    = zero;
+    for (unsigned int mu = 0; mu < Nd; ++mu)
+    {
+        act.ContractConservedCurrent(q, q, tmp, Current::Vector, mu);
+        tmp -= Cshift(tmp, mu, -1);
+        vector_WI += tmp;
+    }
+
+    // Test ward identity D_mu V_mu = 0;
+    LOG(Message) << "Vector Ward Identity check Delta_mu V_mu = " 
+                 << norm2(vector_WI) << std::endl;
+
+    if (par().test_axial)
+    {
+        envGetTmp(PropagatorField, psi);
+        envGetTmp(LatticeComplex, PP);
+        envGetTmp(LatticeComplex, axial_defect);
+        envGetTmp(LatticeComplex, PJ5q);
+        std::vector<TComplex> axial_buf;
+
+        // Compute <P|D_mu A_mu>, D is backwards derivative.
+        axial_defect = zero;
+        for (unsigned int mu = 0; mu < Nd; ++mu)
+        {
+            act.ContractConservedCurrent(q, q, tmp, Current::Axial, mu);
+            tmp -= Cshift(tmp, mu, -1);
+            axial_defect += trace(g5*tmp);
+        }
+
+        // Get <P|J5q> for 5D (zero for 4D) and <P|P>.
+        PJ5q = zero;
+        if (Ls_ > 1)
+        {
+            // <P|P>
+            ExtractSlice(tmp, q, 0, 0);
+            psi  = 0.5 * (tmp - g5*tmp);
+            ExtractSlice(tmp, q, Ls_ - 1, 0);
+            psi += 0.5 * (tmp + g5*tmp);
+            PP = trace(adj(psi)*psi);
+
+            // <P|5Jq>
+            ExtractSlice(tmp, q, Ls_/2 - 1, 0);
+            psi  = 0.5 * (tmp + g5*tmp);
+            ExtractSlice(tmp, q, Ls_/2, 0);
+            psi += 0.5 * (tmp - g5*tmp);
+            PJ5q = trace(adj(psi)*psi);
+        }
+        else
+        {
+            PP = trace(adj(q)*q);
+        }
+
+        // Test ward identity <P|D_mu A_mu> = 2m<P|P> + 2<P|J5q>
+        LOG(Message) << "|D_mu A_mu|^2 = " << norm2(axial_defect) << std::endl;
+        LOG(Message) << "|PP|^2        = " << norm2(PP) << std::endl;
+        LOG(Message) << "|PJ5q|^2      = " << norm2(PJ5q) << std::endl;
+        LOG(Message) << "Axial Ward Identity defect Delta_mu A_mu = "
+                     << norm2(axial_defect) << std::endl;
+    
+        // Axial defect by timeslice.
+        axial_defect -= 2.*(par().mass*PP + PJ5q);
+        LOG(Message) << "Check Axial defect by timeslice" << std::endl;
+        sliceSum(axial_defect, axial_buf, Tp);
+        for (int t = 0; t < axial_buf.size(); ++t)
+        {
+            LOG(Message) << "t = " << t << ": " 
+                         << TensorRemove(axial_buf[t]) << std::endl;
+        }
+    }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WardIdentity_hpp_

extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp

@@ -76,6 +76,7 @@ public:
                                     std::string, q2,
                                     std::string, q3,
                                     std::string, q4,
+                                    unsigned int, tSnk,
                                     std::string, output);
 };
 
@@ -99,11 +100,13 @@ public:\
     /* dependency relation */ \
     virtual std::vector<std::string> getInput(void);\
     virtual std::vector<std::string> getOutput(void);\
+public:\
+    std::vector<std::string> VA_label = {"V", "A"};\
+protected:\
     /* setup */ \
     virtual void setup(void);\
     /* execution */ \
     virtual void execute(void);\
-    std::vector<std::string> VA_label = {"V", "A"};\
 };\
 MODULE_REGISTER_NS(modname, T##modname, MContraction);
 

extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc

@@ -54,6 +54,8 @@ using namespace MContraction;
  * 
  * S: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1]*q4*gL[mu][p_2])
  * E: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1])*trace(q4*gL[mu][p_2])
+ * 
+ * Note q1 must be sink smeared.
  */
 
 /******************************************************************************
@@ -74,7 +76,7 @@ std::vector<std::string> TWeakHamiltonianEye::getInput(void)
 
 std::vector<std::string> TWeakHamiltonianEye::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {};
     
     return out;
 }
@@ -82,7 +84,15 @@ std::vector<std::string> TWeakHamiltonianEye::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TWeakHamiltonianEye::setup(void)
 {
+    unsigned int ndim = env().getNd();
 
+    envTmpLat(LatticeComplex,  "expbuf");
+    envTmpLat(PropagatorField, "tmp1");
+    envTmpLat(LatticeComplex,  "tmp2");
+    envTmp(std::vector<PropagatorField>, "S_body", 1, ndim, PropagatorField(env().getGrid()));
+    envTmp(std::vector<PropagatorField>, "S_loop", 1, ndim, PropagatorField(env().getGrid()));
+    envTmp(std::vector<LatticeComplex>,  "E_body", 1, ndim, LatticeComplex(env().getGrid()));
+    envTmp(std::vector<LatticeComplex>,  "E_loop", 1, ndim, LatticeComplex(env().getGrid()));
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -94,27 +104,30 @@ void TWeakHamiltonianEye::execute(void)
                  << "'." << std::endl;
 
     CorrWriter             writer(par().output);
-    PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
+    auto                   &q1 = envGet(SlicedPropagator, par().q1);
-    PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
+    auto                   &q2 = envGet(PropagatorField, par().q2);
-    PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
+    auto                   &q3 = envGet(PropagatorField, par().q3);
-    PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
+    auto                   &q4 = envGet(PropagatorField, par().q4);
-    Gamma g5            = Gamma(Gamma::Algebra::Gamma5);
+    Gamma                  g5  = Gamma(Gamma::Algebra::Gamma5);
-    LatticeComplex        expbuf(env().getGrid());
+    std::vector<TComplex>  corrbuf;
-    std::vector<TComplex> corrbuf;
+    std::vector<Result>    result(n_eye_diag);
-    std::vector<Result>   result(n_eye_diag);
+    unsigned int ndim    = env().getNd();
-    unsigned int ndim   = env().getNd();
 
-    PropagatorField              tmp1(env().getGrid());
+    envGetTmp(LatticeComplex,               expbuf); 
-    LatticeComplex               tmp2(env().getGrid());
+    envGetTmp(PropagatorField,              tmp1);
-    std::vector<PropagatorField> S_body(ndim, tmp1);
+    envGetTmp(LatticeComplex,               tmp2);
-    std::vector<PropagatorField> S_loop(ndim, tmp1);
+    envGetTmp(std::vector<PropagatorField>, S_body);
-    std::vector<LatticeComplex>  E_body(ndim, tmp2);
+    envGetTmp(std::vector<PropagatorField>, S_loop);
-    std::vector<LatticeComplex>  E_loop(ndim, tmp2);
+    envGetTmp(std::vector<LatticeComplex>,  E_body);
+    envGetTmp(std::vector<LatticeComplex>,  E_loop);
+
+    // Get sink timeslice of q1.
+    SitePropagator q1Snk = q1[par().tSnk];
 
     // Setup for S-type contractions.
     for (int mu = 0; mu < ndim; ++mu)
     {
-        S_body[mu] = MAKE_SE_BODY(q1, q2, q3, GammaL(Gamma::gmu[mu]));
+        S_body[mu] = MAKE_SE_BODY(q1Snk, q2, q3, GammaL(Gamma::gmu[mu]));
         S_loop[mu] = MAKE_SE_LOOP(q4, GammaL(Gamma::gmu[mu]));
     }
 

extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc

@@ -76,7 +76,7 @@ std::vector<std::string> TWeakHamiltonianNonEye::getInput(void)
 
 std::vector<std::string> TWeakHamiltonianNonEye::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {};
     
     return out;
 }
@@ -84,7 +84,15 @@ std::vector<std::string> TWeakHamiltonianNonEye::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TWeakHamiltonianNonEye::setup(void)
 {
+    unsigned int ndim = env().getNd();
 
+    envTmpLat(LatticeComplex,  "expbuf");
+    envTmpLat(PropagatorField, "tmp1");
+    envTmpLat(LatticeComplex,  "tmp2");
+    envTmp(std::vector<PropagatorField>, "C_i_side_loop", 1, ndim, PropagatorField(env().getGrid()));
+    envTmp(std::vector<PropagatorField>, "C_f_side_loop", 1, ndim, PropagatorField(env().getGrid()));
+    envTmp(std::vector<LatticeComplex>,  "W_i_side_loop", 1, ndim, LatticeComplex(env().getGrid()));
+    envTmp(std::vector<LatticeComplex>,  "W_f_side_loop", 1, ndim, LatticeComplex(env().getGrid()));
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -95,23 +103,23 @@ void TWeakHamiltonianNonEye::execute(void)
                  << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
                  << "'." << std::endl;
     
-    CorrWriter             writer(par().output);
+    CorrWriter            writer(par().output);
-    PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
+    auto                  &q1 = envGet(PropagatorField, par().q1);
-    PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
+    auto                  &q2 = envGet(PropagatorField, par().q2);
-    PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
+    auto                  &q3 = envGet(PropagatorField, par().q3);
-    PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
+    auto                  &q4 = envGet(PropagatorField, par().q4);
-    Gamma g5            = Gamma(Gamma::Algebra::Gamma5);
+    Gamma                 g5  = Gamma(Gamma::Algebra::Gamma5);
-    LatticeComplex        expbuf(env().getGrid());
     std::vector<TComplex> corrbuf;
     std::vector<Result>   result(n_noneye_diag); 
-    unsigned int ndim   = env().getNd();
+    unsigned int          ndim = env().getNd();
 
-    PropagatorField              tmp1(env().getGrid());
+    envGetTmp(LatticeComplex,               expbuf); 
-    LatticeComplex               tmp2(env().getGrid());
+    envGetTmp(PropagatorField,              tmp1);
-    std::vector<PropagatorField> C_i_side_loop(ndim, tmp1);
+    envGetTmp(LatticeComplex,               tmp2);
-    std::vector<PropagatorField> C_f_side_loop(ndim, tmp1);
+    envGetTmp(std::vector<PropagatorField>, C_i_side_loop);
-    std::vector<LatticeComplex>  W_i_side_loop(ndim, tmp2);
+    envGetTmp(std::vector<PropagatorField>, C_f_side_loop);
-    std::vector<LatticeComplex>  W_f_side_loop(ndim, tmp2);
+    envGetTmp(std::vector<LatticeComplex>,  W_i_side_loop);
+    envGetTmp(std::vector<LatticeComplex>,  W_f_side_loop);
 
     // Setup for C-type contractions.
     for (int mu = 0; mu < ndim; ++mu)

extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc

@@ -78,7 +78,7 @@ std::vector<std::string> TWeakNeutral4ptDisc::getInput(void)
 
 std::vector<std::string> TWeakNeutral4ptDisc::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {};
     
     return out;
 }
@@ -86,7 +86,13 @@ std::vector<std::string> TWeakNeutral4ptDisc::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TWeakNeutral4ptDisc::setup(void)
 {
+    unsigned int ndim = env().getNd();
 
+    envTmpLat(LatticeComplex,  "expbuf");
+    envTmpLat(PropagatorField, "tmp");
+    envTmpLat(LatticeComplex,  "curr");
+    envTmp(std::vector<PropagatorField>, "meson", 1, ndim, PropagatorField(env().getGrid()));
+    envTmp(std::vector<PropagatorField>, "loop", 1, ndim,  PropagatorField(env().getGrid()));
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -97,21 +103,21 @@ void TWeakNeutral4ptDisc::execute(void)
                  << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
                  << "'." << std::endl;
 
-    CorrWriter             writer(par().output);
+    CorrWriter            writer(par().output);
-    PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
+    auto                  &q1 = envGet(PropagatorField, par().q1);
-    PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
+    auto                  &q2 = envGet(PropagatorField, par().q2);
-    PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
+    auto                  &q3 = envGet(PropagatorField, par().q3);
-    PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
+    auto                  &q4 = envGet(PropagatorField, par().q4);
-    Gamma g5            = Gamma(Gamma::Algebra::Gamma5);
+    Gamma                 g5  = Gamma(Gamma::Algebra::Gamma5);
-    LatticeComplex        expbuf(env().getGrid());
     std::vector<TComplex> corrbuf;
     std::vector<Result>   result(n_neut_disc_diag);
-    unsigned int ndim   = env().getNd();
+    unsigned int          ndim = env().getNd();
 
-    PropagatorField              tmp(env().getGrid());
+    envGetTmp(LatticeComplex,               expbuf); 
-    std::vector<PropagatorField> meson(ndim, tmp);
+    envGetTmp(PropagatorField,              tmp);
-    std::vector<PropagatorField> loop(ndim, tmp);
+    envGetTmp(LatticeComplex,               curr);
-    LatticeComplex               curr(env().getGrid());
+    envGetTmp(std::vector<PropagatorField>, meson);
+    envGetTmp(std::vector<PropagatorField>, loop);
 
     // Setup for type 1 contractions.
     for (int mu = 0; mu < ndim; ++mu)

extras/Hadrons/Modules/MFermion/GaugeProp.hpp

@@ -1,3 +1,34 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MFermion/GaugeProp.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+Copyright (C) 2017
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+        Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
 #ifndef Hadrons_MFermion_GaugeProp_hpp_
 #define Hadrons_MFermion_GaugeProp_hpp_
 
@@ -7,6 +38,27 @@
 
 BEGIN_HADRONS_NAMESPACE
 
+/******************************************************************************
+ * 5D -> 4D and 4D -> 5D conversions.                                         *
+ ******************************************************************************/
+template<class vobj> // Note that 5D object is modified.
+inline void make_4D(Lattice<vobj> &in_5d, Lattice<vobj> &out_4d, int Ls)
+{
+    axpby_ssp_pminus(in_5d, 0., in_5d, 1., in_5d, 0, 0);
+    axpby_ssp_pplus(in_5d, 1., in_5d, 1., in_5d, 0, Ls-1);
+    ExtractSlice(out_4d, in_5d, 0, 0);
+}
+
+template<class vobj>
+inline void make_5D(Lattice<vobj> &in_4d, Lattice<vobj> &out_5d, int Ls)
+{
+    out_5d = zero;
+    InsertSlice(in_4d, out_5d, 0, 0);
+    InsertSlice(in_4d, out_5d, Ls-1, 0);
+    axpby_ssp_pplus(out_5d, 0., out_5d, 1., out_5d, 0, 0);
+    axpby_ssp_pminus(out_5d, 0., out_5d, 1., out_5d, Ls-1, Ls-1);
+}
+
 /******************************************************************************
  *                                GaugeProp                                   *
  ******************************************************************************/
@@ -33,6 +85,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
@@ -75,10 +128,13 @@ template <typename FImpl>
 void TGaugeProp<FImpl>::setup(void)
 {
     Ls_ = env().getObjectLs(par().solver);
-    env().template registerLattice<PropagatorField>(getName());
+    envCreateLat(PropagatorField, getName());
+    envTmpLat(FermionField, "source", Ls_);
+    envTmpLat(FermionField, "sol", Ls_);
+    envTmpLat(FermionField, "tmp");
     if (Ls_ > 1)
     {
-        env().template registerLattice<PropagatorField>(getName() + "_5d", Ls_);
+        envCreateLat(PropagatorField, getName() + "_5d", Ls_);
     }
 }
 
@@ -87,26 +143,23 @@ template <typename FImpl>
 void TGaugeProp<FImpl>::execute(void)
 {
     LOG(Message) << "Computing quark propagator '" << getName() << "'"
-    << std::endl;
+                 << std::endl;
     
-    FermionField    source(env().getGrid(Ls_)), sol(env().getGrid(Ls_)),
+    std::string propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
-    tmp(env().getGrid());
+    auto        &prop    = envGet(PropagatorField, propName);
-    std::string     propName = (Ls_ == 1) ? getName() : (getName() + "_5d");
+    auto        &fullSrc = envGet(PropagatorField, par().source);
-    PropagatorField &prop    = *env().template createLattice<PropagatorField>(propName);
+    auto        &solver  = envGet(SolverFn, par().solver);
-    PropagatorField &fullSrc = *env().template getObject<PropagatorField>(par().source);
-    SolverFn        &solver  = *env().template getObject<SolverFn>(par().solver);
-    if (Ls_ > 1)
-    {
-        env().template createLattice<PropagatorField>(getName());
-    }
     
+    envGetTmp(FermionField, source);
+    envGetTmp(FermionField, sol);
+    envGetTmp(FermionField, tmp);
     LOG(Message) << "Inverting using solver '" << par().solver
-    << "' on source '" << par().source << "'" << std::endl;
+                 << "' on source '" << par().source << "'" << std::endl;
     for (unsigned int s = 0; s < Ns; ++s)
     for (unsigned int c = 0; c < Nc; ++c)
     {
         LOG(Message) << "Inversion for spin= " << s << ", color= " << c
-        << std::endl;
+                     << std::endl;
         // source conversion for 4D sources
         if (!env().isObject5d(par().source))
         {
@@ -116,12 +169,8 @@ void TGaugeProp<FImpl>::execute(void)
             }
             else
             {
-                source = zero;
                 PropToFerm(tmp, fullSrc, s, c);
-                InsertSlice(tmp, source, 0, 0);
+                make_5D(tmp, source, Ls_);
-                InsertSlice(tmp, source, Ls_-1, 0);
-                axpby_ssp_pplus(source, 0., source, 1., source, 0, 0);
-                axpby_ssp_pminus(source, 0., source, 1., source, Ls_-1, Ls_-1);
             }
         }
         // source conversion for 5D sources
@@ -129,7 +178,7 @@ void TGaugeProp<FImpl>::execute(void)
         {
             if (Ls_ != env().getObjectLs(par().source))
             {
-                HADRON_ERROR("Ls mismatch between quark action and source");
+                HADRON_ERROR(Size, "Ls mismatch between quark action and source");
             }
             else
             {
@@ -142,12 +191,8 @@ void TGaugeProp<FImpl>::execute(void)
         // create 4D propagators from 5D one if necessary
         if (Ls_ > 1)
         {
-            PropagatorField &p4d =
+            PropagatorField &p4d = envGet(PropagatorField, getName());
-            *env().template getObject<PropagatorField>(getName());
+            make_4D(sol, tmp, Ls_);
-            
-            axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0);
-            axpby_ssp_pplus(sol, 1., sol, 1., sol, 0, Ls_-1);
-            ExtractSlice(tmp, sol, 0, 0);
             FermToProp(p4d, tmp, s, c);
         }
     }

extras/Hadrons/Modules/MGauge/Load.cc

@@ -59,19 +59,19 @@ std::vector<std::string> TLoad::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TLoad::setup(void)
 {
-    env().registerLattice<LatticeGaugeField>(getName());
+    envCreateLat(LatticeGaugeField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void TLoad::execute(void)
 {
-    FieldMetaData  header;
+    FieldMetaData header;
-    std::string fileName = par().file + "."
+    std::string   fileName = par().file + "."
-                           + std::to_string(env().getTrajectory());
+                             + std::to_string(vm().getTrajectory());
-    
     LOG(Message) << "Loading NERSC configuration from file '" << fileName
                  << "'" << std::endl;
-    LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
+
+    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/MGauge/Load.hpp

@@ -58,6 +58,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution

extras/Hadrons/Modules/MGauge/Random.cc

@@ -44,7 +44,9 @@ TRandom::TRandom(const std::string name)
 // dependencies/products ///////////////////////////////////////////////////////
 std::vector<std::string> TRandom::getInput(void)
 {
-    return std::vector<std::string>();
+    std::vector<std::string> in;
+    
+    return in;
 }
 
 std::vector<std::string> TRandom::getOutput(void)
@@ -57,13 +59,14 @@ std::vector<std::string> TRandom::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TRandom::setup(void)
 {
-    env().registerLattice<LatticeGaugeField>(getName());
+    envCreateLat(LatticeGaugeField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void TRandom::execute(void)
 {
     LOG(Message) << "Generating random gauge configuration" << std::endl;
-    LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
+    
+    auto &U = envGet(LatticeGaugeField, getName());
     SU3::HotConfiguration(*env().get4dRng(), U);
 }

extras/Hadrons/Modules/MGauge/Random.hpp

@@ -51,6 +51,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution

extras/Hadrons/Modules/MGauge/StochEm.cc

@@ -57,32 +57,28 @@ std::vector<std::string> TStochEm::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TStochEm::setup(void)
 {
-    if (!env().hasRegisteredObject("_" + getName() + "_weight"))
+    if (!env().hasCreatedObject("_" + getName() + "_weight"))
     {
-        env().registerLattice<EmComp>("_" + getName() + "_weight");
+        envCacheLat(EmComp, "_" + getName() + "_weight");
     }
-    env().registerLattice<EmField>(getName());
+    envCreateLat(EmField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void TStochEm::execute(void)
 {
+    LOG(Message) << "Generating stochatic EM potential..." << std::endl;
+
     PhotonR photon(par().gauge, par().zmScheme);
-    EmField &a = *env().createLattice<EmField>(getName());
+    auto    &a = envGet(EmField, getName());
-    EmComp  *w;
+    auto    &w = envGet(EmComp, "_" + getName() + "_weight");
     
     if (!env().hasCreatedObject("_" + getName() + "_weight"))
     {
         LOG(Message) << "Caching stochatic EM potential weight (gauge: "
                      << par().gauge << ", zero-mode scheme: "
                      << par().zmScheme << ")..." << std::endl;
-        w = env().createLattice<EmComp>("_" + getName() + "_weight");
+        photon.StochasticWeight(w);
-        photon.StochasticWeight(*w);
     }
-    else
+    photon.StochasticField(a, *env().get4dRng(), w);
-    {
-        w = env().getObject<EmComp>("_" + getName() + "_weight");
-    }
-    LOG(Message) << "Generating stochatic EM potential..." << std::endl;
-    photon.StochasticField(a, *env().get4dRng(), *w);
 }

extras/Hadrons/Modules/MGauge/StochEm.hpp

@@ -60,6 +60,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution

extras/Hadrons/Modules/MGauge/Unit.cc

@@ -57,13 +57,14 @@ std::vector<std::string> TUnit::getOutput(void)
 // setup ///////////////////////////////////////////////////////////////////////
 void TUnit::setup(void)
 {
-    env().registerLattice<LatticeGaugeField>(getName());
+    envCreateLat(LatticeGaugeField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void TUnit::execute(void)
 {
     LOG(Message) << "Creating unit gauge configuration" << std::endl;
-    LatticeGaugeField &U = *env().createLattice<LatticeGaugeField>(getName());
+    
+    auto &U = envGet(LatticeGaugeField, getName());
     SU3::ColdConfiguration(*env().get4dRng(), U);
 }

extras/Hadrons/Modules/MGauge/Unit.hpp

@@ -51,6 +51,7 @@ public:
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution

extras/Hadrons/Modules/MLoop/NoiseLoop.hpp

@@ -74,6 +74,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
@@ -112,16 +113,16 @@ std::vector<std::string> TNoiseLoop<FImpl>::getOutput(void)
 template <typename FImpl>
 void TNoiseLoop<FImpl>::setup(void)
 {
-    env().template registerLattice<PropagatorField>(getName());
+    envCreateLat(PropagatorField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TNoiseLoop<FImpl>::execute(void)
 {
-    PropagatorField &loop = *env().template createLattice<PropagatorField>(getName());
+    auto &loop = envGet(PropagatorField, getName());
-    PropagatorField &q    = *env().template getObject<PropagatorField>(par().q);
+    auto &q    = envGet(PropagatorField, par().q);
-    PropagatorField &eta  = *env().template getObject<PropagatorField>(par().eta);
+    auto &eta  = envGet(PropagatorField, par().eta);
     loop = q*adj(eta);
 }
 

extras/Hadrons/Modules/MScalar/ChargedProp.cc

@@ -37,90 +37,44 @@ void TChargedProp::setup(void)
     {
         phaseName_.push_back("_shiftphase_" + std::to_string(mu));
     }
-    GFSrcName_ = "_" + getName() + "_DinvSrc";
+    GFSrcName_ = getName() + "_DinvSrc";
-    if (!env().hasRegisteredObject(freeMomPropName_))
+    fftName_   = getName() + "_fft";
+
+    freeMomPropDone_ = env().hasCreatedObject(freeMomPropName_);
+    GFSrcDone_       = env().hasCreatedObject(GFSrcName_);
+    phasesDone_      = env().hasCreatedObject(phaseName_[0]);
+    envCacheLat(ScalarField, freeMomPropName_);
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
     {
-        env().registerLattice<ScalarField>(freeMomPropName_);
+        envCacheLat(ScalarField, phaseName_[mu]);
     }
-    if (!env().hasRegisteredObject(phaseName_[0]))
+    envCacheLat(ScalarField, GFSrcName_);
-    {
+    envCreateLat(ScalarField, getName());
-        for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    envTmpLat(ScalarField, "buf");
-        {
+    envTmpLat(ScalarField, "result");
-            env().registerLattice<ScalarField>(phaseName_[mu]);
+    envTmpLat(ScalarField, "Amu");
-        }
+    envCache(FFT, fftName_, 1, env().getGrid());
-    }
-    if (!env().hasRegisteredObject(GFSrcName_))
-    {
-        env().registerLattice<ScalarField>(GFSrcName_);
-    }
-    env().registerLattice<ScalarField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void TChargedProp::execute(void)
 {
     // CACHING ANALYTIC EXPRESSIONS
-    ScalarField &source = *env().getObject<ScalarField>(par().source);
+    makeCaches();
-    Complex     ci(0.0,1.0);
-    FFT         fft(env().getGrid());
-    
-    // cache free scalar propagator
-    if (!env().hasCreatedObject(freeMomPropName_))
-    {
-        LOG(Message) << "Caching momentum space free scalar propagator"
-                     << " (mass= " << par().mass << ")..." << std::endl;
-        freeMomProp_ = env().createLattice<ScalarField>(freeMomPropName_);
-        SIMPL::MomentumSpacePropagator(*freeMomProp_, par().mass);
-    }
-    else
-    {
-        freeMomProp_ = env().getObject<ScalarField>(freeMomPropName_);
-    }
-    // cache G*F*src
-    if (!env().hasCreatedObject(GFSrcName_))
-        
-    {
-        GFSrc_ = env().createLattice<ScalarField>(GFSrcName_);
-        fft.FFT_all_dim(*GFSrc_, source, FFT::forward);
-        *GFSrc_ = (*freeMomProp_)*(*GFSrc_);
-    }
-    else
-    {
-        GFSrc_ = env().getObject<ScalarField>(GFSrcName_);
-    }
-    // cache phases
-    if (!env().hasCreatedObject(phaseName_[0]))
-    {
-        std::vector<int> &l = env().getGrid()->_fdimensions;
-        
-        LOG(Message) << "Caching shift phases..." << std::endl;
-        for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-        {
-            Real    twoPiL = M_PI*2./l[mu];
-            
-            phase_.push_back(env().createLattice<ScalarField>(phaseName_[mu]));
-            LatticeCoordinate(*(phase_[mu]), mu);
-            *(phase_[mu]) = exp(ci*twoPiL*(*(phase_[mu])));
-        }
-    }
-    else
-    {
-        for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-        {
-            phase_.push_back(env().getObject<ScalarField>(phaseName_[mu]));
-        }
-    }
 
     // PROPAGATOR CALCULATION
     LOG(Message) << "Computing charged scalar propagator"
                  << " (mass= " << par().mass
                  << ", charge= " << par().charge << ")..." << std::endl;
     
-    ScalarField &prop   = *env().createLattice<ScalarField>(getName());
+    auto   &prop  = envGet(ScalarField, getName());
-    ScalarField buf(env().getGrid());
+    auto   &GFSrc = envGet(ScalarField, GFSrcName_);
-    ScalarField &GFSrc = *GFSrc_, &G = *freeMomProp_;
+    auto   &G     = envGet(ScalarField, freeMomPropName_);
-    double      q = par().charge;
+    auto   &fft   = envGet(FFT, fftName_);
-    
+    double q      = par().charge;
+    envGetTmp(ScalarField, result); 
+    envGetTmp(ScalarField, buf); 
+
     // G*F*Src
     prop = GFSrc;
 
@@ -146,7 +100,7 @@ void TChargedProp::execute(void)
     if (!par().output.empty())
     {
         std::string           filename = par().output + "." +
-                                         std::to_string(env().getTrajectory());
+                                         std::to_string(vm().getTrajectory());
         
         LOG(Message) << "Saving zero-momentum projection to '"
                      << filename << "'..." << std::endl;
@@ -166,15 +120,55 @@ void TChargedProp::execute(void)
     }
 }
 
+void TChargedProp::makeCaches(void)
+{
+    auto &freeMomProp = envGet(ScalarField, freeMomPropName_);
+    auto &GFSrc       = envGet(ScalarField, GFSrcName_);
+    auto &fft         = envGet(FFT, fftName_);
+
+    if (!freeMomPropDone_)
+    {
+        LOG(Message) << "Caching momentum space free scalar propagator"
+                     << " (mass= " << par().mass << ")..." << std::endl;
+        SIMPL::MomentumSpacePropagator(freeMomProp, par().mass);
+    }
+    if (!GFSrcDone_)
+    {   
+        FFT  fft(env().getGrid());
+        auto &source = envGet(ScalarField, par().source);
+
+        LOG(Message) << "Caching G*F*src..." << std::endl;
+        fft.FFT_all_dim(GFSrc, source, FFT::forward);
+        GFSrc = freeMomProp*GFSrc;
+    }
+    if (!phasesDone_)
+    {
+        std::vector<int> &l = env().getGrid()->_fdimensions;
+        Complex          ci(0.0,1.0);
+        
+        LOG(Message) << "Caching shift phases..." << std::endl;
+        for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+        {
+            Real twoPiL = M_PI*2./l[mu];
+            auto &phmu  = envGet(ScalarField, phaseName_[mu]);
+            
+            LatticeCoordinate(phmu, mu);
+            phmu = exp(ci*twoPiL*phmu);
+            phase_.push_back(&phmu);
+        }
+    }
+}
+
 void TChargedProp::momD1(ScalarField &s, FFT &fft)
 {
-    EmField     &A = *env().getObject<EmField>(par().emField);
+    auto        &A = envGet(EmField, par().emField);
-    ScalarField buf(env().getGrid()), result(env().getGrid()),
-                Amu(env().getGrid());
     Complex     ci(0.0,1.0);
 
-    result = zero;
+    envGetTmp(ScalarField, buf);
+    envGetTmp(ScalarField, result);
+    envGetTmp(ScalarField, Amu);
 
+    result = zero;
     for (unsigned int mu = 0; mu < env().getNd(); ++mu)
     {
         Amu = peekLorentz(A, mu);
@@ -198,12 +192,13 @@ void TChargedProp::momD1(ScalarField &s, FFT &fft)
 
 void TChargedProp::momD2(ScalarField &s, FFT &fft)
 {
-    EmField     &A = *env().getObject<EmField>(par().emField);
+    auto &A = envGet(EmField, par().emField);
-    ScalarField buf(env().getGrid()), result(env().getGrid()),
+
-                Amu(env().getGrid());
+    envGetTmp(ScalarField, buf);
+    envGetTmp(ScalarField, result);
+    envGetTmp(ScalarField, Amu);
 
     result = zero;
-    
     for (unsigned int mu = 0; mu < env().getNd(); ++mu)
     {
         Amu = peekLorentz(A, mu);

extras/Hadrons/Modules/MScalar/ChargedProp.hpp

@@ -37,19 +37,20 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
     virtual void execute(void);
 private:
+    void makeCaches(void);
     void momD1(ScalarField &s, FFT &fft);
     void momD2(ScalarField &s, FFT &fft);
 private:
-    std::string                freeMomPropName_, GFSrcName_;
+    bool                       freeMomPropDone_, GFSrcDone_, phasesDone_;
+    std::string                freeMomPropName_, GFSrcName_, fftName_;
     std::vector<std::string>   phaseName_;
-    ScalarField                *freeMomProp_, *GFSrc_;
     std::vector<ScalarField *> phase_;
-    EmField                    *A;
 };
 
 MODULE_REGISTER_NS(ChargedProp, TChargedProp, MScalar);

extras/Hadrons/Modules/MScalar/FreeProp.cc

@@ -33,38 +33,31 @@ void TFreeProp::setup(void)
 {
     freeMomPropName_ = FREEMOMPROP(par().mass);
     
-    if (!env().hasRegisteredObject(freeMomPropName_))
+    freePropDone_ = env().hasCreatedObject(freeMomPropName_);
-    {
+    envCacheLat(ScalarField, freeMomPropName_);
-        env().registerLattice<ScalarField>(freeMomPropName_);
+    envCreateLat(ScalarField, getName());
-    }
-    env().registerLattice<ScalarField>(getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 void TFreeProp::execute(void)
 {
-    ScalarField &prop   = *env().createLattice<ScalarField>(getName());
+    auto &freeMomProp = envGet(ScalarField, freeMomPropName_);
-    ScalarField &source = *env().getObject<ScalarField>(par().source);
+    auto &prop        = envGet(ScalarField, getName());
-    ScalarField *freeMomProp;
+    auto &source      = envGet(ScalarField, par().source);
 
-    if (!env().hasCreatedObject(freeMomPropName_))
+    if (!freePropDone_)
     {
         LOG(Message) << "Caching momentum space free scalar propagator"
                      << " (mass= " << par().mass << ")..." << std::endl;
-        freeMomProp = env().createLattice<ScalarField>(freeMomPropName_);
+        SIMPL::MomentumSpacePropagator(freeMomProp, par().mass);
-        SIMPL::MomentumSpacePropagator(*freeMomProp, par().mass);
-    }
-    else
-    {
-        freeMomProp = env().getObject<ScalarField>(freeMomPropName_);
     }
     LOG(Message) << "Computing free scalar propagator..." << std::endl;
-    SIMPL::FreePropagator(source, prop, *freeMomProp);
+    SIMPL::FreePropagator(source, prop, freeMomProp);
     
     if (!par().output.empty())
     {
         TextWriter            writer(par().output + "." +
-                                     std::to_string(env().getTrajectory()));
+                                     std::to_string(vm().getTrajectory()));
         std::vector<TComplex> buf;
         std::vector<Complex>  result;
         

extras/Hadrons/Modules/MScalar/FreeProp.hpp

@@ -33,12 +33,14 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
     virtual void execute(void);
 private:
     std::string freeMomPropName_;
+    bool        freePropDone_;
 };
 
 MODULE_REGISTER_NS(FreeProp, TFreeProp, MScalar);

extras/Hadrons/Modules/MSink/Point.hpp

@@ -1,3 +1,31 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSink/Point.hpp
+
+Copyright (C) 2017
+
+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_MSink_Point_hpp_
 #define Hadrons_MSink_Point_hpp_
 
@@ -33,10 +61,14 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
     virtual void execute(void);
+private:
+    bool        hasPhase_{false}; 
+    std::string momphName_;
 };
 
 MODULE_REGISTER_NS(Point,       TPoint<FIMPL>,        MSink);
@@ -49,6 +81,7 @@ MODULE_REGISTER_NS(ScalarPoint, TPoint<ScalarImplCR>, MSink);
 template <typename FImpl>
 TPoint<FImpl>::TPoint(const std::string name)
 : Module<PointPar>(name)
+, momphName_ (name + "_momph")
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
@@ -72,30 +105,37 @@ std::vector<std::string> TPoint<FImpl>::getOutput(void)
 template <typename FImpl>
 void TPoint<FImpl>::setup(void)
 {
-    unsigned int size;
+    envTmpLat(LatticeComplex, "coor");
-    
+    envCacheLat(LatticeComplex, momphName_);
-    size = env().template lattice4dSize<LatticeComplex>();
+    envCreate(SinkFn, getName(), 1, nullptr);
-    env().registerObject(getName(), size);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TPoint<FImpl>::execute(void)
-{
+{   
-    std::vector<Real> p = strToVec<Real>(par().mom);
-    LatticeComplex    ph(env().getGrid()), coor(env().getGrid());
-    Complex           i(0.0,1.0);
-    
     LOG(Message) << "Setting up point sink function for momentum ["
                  << par().mom << "]" << std::endl;
-    ph = zero;
+
-    for(unsigned int mu = 0; mu < env().getNd(); mu++)
+    auto &ph = envGet(LatticeComplex, momphName_);
+    
+    if (!hasPhase_)
     {
-        LatticeCoordinate(coor, mu);
+        Complex           i(0.0,1.0);
-        ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
+        std::vector<Real> p;
+
+        envGetTmp(LatticeComplex, coor);
+        p  = strToVec<Real>(par().mom);
+        ph = zero;
+        for(unsigned int mu = 0; mu < env().getNd(); mu++)
+        {
+            LatticeCoordinate(coor, mu);
+            ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
+        }
+        ph = exp((Real)(2*M_PI)*i*ph);
+        hasPhase_ = true;
     }
-    ph = exp((Real)(2*M_PI)*i*ph);
+    auto sink = [&ph](const PropagatorField &field)
-    auto sink = [ph](const PropagatorField &field)
     {
         SlicedPropagator res;
         PropagatorField  tmp = ph*field;
@@ -104,7 +144,7 @@ void TPoint<FImpl>::execute(void)
         
         return res;
     };
-    env().setObject(getName(), new SinkFn(sink));
+    envGet(SinkFn, getName()) = sink;
 }
 
 END_MODULE_NAMESPACE

extras/Hadrons/Modules/MSink/Smear.hpp

@@ -0,0 +1,126 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSink/Smear.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_MSink_Smear_hpp_
+#define Hadrons_MSink_Smear_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                 Smear                                      *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSink)
+
+class SmearPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SmearPar,
+                                    std::string, q,
+                                    std::string, sink);
+};
+
+template <typename FImpl>
+class TSmear: public Module<SmearPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    SINK_TYPE_ALIASES();
+public:
+    // constructor
+    TSmear(const std::string name);
+    // destructor
+    virtual ~TSmear(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+protected:
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Smear, TSmear<FIMPL>, MSink);
+
+/******************************************************************************
+ *                          TSmear implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TSmear<FImpl>::TSmear(const std::string name)
+: Module<SmearPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TSmear<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q, par().sink};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TSmear<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSmear<FImpl>::setup(void)
+{
+    envCreate(SlicedPropagator, getName(), 1, env().getDim(Tp));
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSmear<FImpl>::execute(void)
+{
+    LOG(Message) << "Sink smearing propagator '" << par().q
+                 << "' using sink function '" << par().sink << "'."
+                 << std::endl;
+
+    auto &sink = envGet(SinkFn, par().sink);
+    auto &q    = envGet(PropagatorField, par().q);
+    auto &out  = envGet(SlicedPropagator, getName());
+    
+    out = sink(q);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MSink_Smear_hpp_

extras/Hadrons/Modules/MSolver/RBPrecCG.hpp

@@ -61,7 +61,9 @@ public:
     virtual ~TRBPrecCG(void) = default;
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getReference(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
@@ -83,11 +85,19 @@ TRBPrecCG<FImpl>::TRBPrecCG(const std::string name)
 template <typename FImpl>
 std::vector<std::string> TRBPrecCG<FImpl>::getInput(void)
 {
-    std::vector<std::string> in = {par().action};
+    std::vector<std::string> in = {};
     
     return in;
 }
 
+template <typename FImpl>
+std::vector<std::string> TRBPrecCG<FImpl>::getReference(void)
+{
+    std::vector<std::string> ref = {par().action};
+    
+    return ref;
+}
+
 template <typename FImpl>
 std::vector<std::string> TRBPrecCG<FImpl>::getOutput(void)
 {
@@ -100,17 +110,12 @@ std::vector<std::string> TRBPrecCG<FImpl>::getOutput(void)
 template <typename FImpl>
 void TRBPrecCG<FImpl>::setup(void)
 {
-    auto Ls = env().getObjectLs(par().action);
+    LOG(Message) << "setting up Schur red-black preconditioned CG for"
-    
+                 << " action '" << par().action << "' with residual "
-    env().registerObject(getName(), 0, Ls);
+                 << par().residual << std::endl;
-    env().addOwnership(getName(), par().action);
-}
 
-// execution ///////////////////////////////////////////////////////////////////
+    auto Ls     = env().getObjectLs(par().action);
-template <typename FImpl>
+    auto &mat   = envGet(FMat, par().action);
-void TRBPrecCG<FImpl>::execute(void)
-{
-    auto &mat   = *(env().template getObject<FMat>(par().action));
     auto solver = [&mat, this](FermionField &sol, const FermionField &source)
     {
         ConjugateGradient<FermionField>           cg(par().residual, 10000);
@@ -118,13 +123,14 @@ void TRBPrecCG<FImpl>::execute(void)
         
         schurSolver(mat, source, sol);
     };
-    
+    envCreate(SolverFn, getName(), Ls, solver);
-    LOG(Message) << "setting up Schur red-black preconditioned CG for"
-                 << " action '" << par().action << "' with residual "
-                 << par().residual << std::endl;
-    env().setObject(getName(), new SolverFn(solver));
 }
 
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TRBPrecCG<FImpl>::execute(void)
+{}
+
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE

extras/Hadrons/Modules/MSource/Point.hpp

@@ -72,6 +72,7 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
@@ -111,19 +112,20 @@ std::vector<std::string> TPoint<FImpl>::getOutput(void)
 template <typename FImpl>
 void TPoint<FImpl>::setup(void)
 {
-    env().template registerLattice<PropagatorField>(getName());
+    envCreateLat(PropagatorField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TPoint<FImpl>::execute(void)
 {
-    std::vector<int> position = strToVec<int>(par().position);
-    typename SitePropagator::scalar_object id;
-    
     LOG(Message) << "Creating point source at position [" << par().position
-                 << "]" << std::endl;
+                << "]" << std::endl;
-    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+
+    std::vector<int> position = strToVec<int>(par().position);
+    auto             &src     = envGet(PropagatorField, getName());
+    SitePropagator   id;
+    
     id  = 1.;
     src = zero;
     pokeSite(id, src, position);

extras/Hadrons/Modules/MSource/SeqConserved.hpp

@@ -0,0 +1,159 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/SeqConserved.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_MSource_SeqConserved_hpp_
+#define Hadrons_MSource_SeqConserved_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Sequential source
+ -----------------------------
+ * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * J_mu * exp(i x.mom)
+ 
+ * options:
+ - q: input propagator (string)
+ - action: fermion action used for propagator q (string)
+ - 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)
+ - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
+ 
+ */
+
+/******************************************************************************
+ *                              SeqConserved                                  *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class SeqConservedPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SeqConservedPar,
+                                    std::string,  q,
+                                    std::string,  action,
+                                    unsigned int, tA,
+                                    unsigned int, tB,
+                                    Current,      curr_type,
+                                    unsigned int, mu,
+                                    std::string,  mom);
+};
+
+template <typename FImpl>
+class TSeqConserved: public Module<SeqConservedPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TSeqConserved(const std::string name);
+    // destructor
+    virtual ~TSeqConserved(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+protected:
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(SeqConserved, TSeqConserved<FIMPL>, MSource);
+
+/******************************************************************************
+ *                      TSeqConserved implementation                          *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TSeqConserved<FImpl>::TSeqConserved(const std::string name)
+: Module<SeqConservedPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TSeqConserved<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q, par().action};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TSeqConserved<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSeqConserved<FImpl>::setup(void)
+{
+    auto Ls_ = env().getObjectLs(par().action);
+    envCreateLat(PropagatorField, getName(), Ls_);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSeqConserved<FImpl>::execute(void)
+{
+    if (par().tA == par().tB)
+    {
+        LOG(Message) << "Generating sequential source with conserved "
+                     << par().curr_type << " current insertion (mu = " 
+                     << par().mu << ") at " << "t = " << par().tA << std::endl;
+    }
+    else
+    {
+        LOG(Message) << "Generating sequential source with conserved "
+                     << par().curr_type << " current insertion (mu = " 
+                     << par().mu << ") for " << par().tA << " <= t <= " 
+                     << par().tB << std::endl;
+    }
+    auto &src = envGet(PropagatorField, getName());
+    auto &q   = envGet(PropagatorField, par().q);
+    auto &mat = envGet(FMat, par().action);
+
+    std::vector<Real> mom = strToVec<Real>(par().mom);
+    mat.SeqConservedCurrent(q, src, par().curr_type, par().mu, 
+                            mom, par().tA, par().tB);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_SeqConserved_hpp_

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -81,10 +81,14 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
     virtual void execute(void);
+private:
+    bool        hasPhase_{false};
+    std::string momphName_, tName_;
 };
 
 MODULE_REGISTER_NS(SeqGamma, TSeqGamma<FIMPL>, MSource);
@@ -96,6 +100,8 @@ MODULE_REGISTER_NS(SeqGamma, TSeqGamma<FIMPL>, MSource);
 template <typename FImpl>
 TSeqGamma<FImpl>::TSeqGamma(const std::string name)
 : Module<SeqGammaPar>(name)
+, momphName_ (name + "_momph")
+, tName_ (name + "_t")
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
@@ -119,7 +125,10 @@ std::vector<std::string> TSeqGamma<FImpl>::getOutput(void)
 template <typename FImpl>
 void TSeqGamma<FImpl>::setup(void)
 {
-    env().template registerLattice<PropagatorField>(getName());
+    envCreateLat(PropagatorField, getName());
+    envCacheLat(Lattice<iScalar<vInteger>>, tName_);
+    envCacheLat(LatticeComplex, momphName_);
+    envTmpLat(LatticeComplex, "coor");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -137,23 +146,29 @@ void TSeqGamma<FImpl>::execute(void)
                      << " sequential source for "
                      << par().tA << " <= t <= " << par().tB << std::endl;
     }
-    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    auto  &src = envGet(PropagatorField, getName());
-    PropagatorField &q   = *env().template getObject<PropagatorField>(par().q);
+    auto  &q   = envGet(PropagatorField, par().q);
-    Lattice<iScalar<vInteger>> t(env().getGrid());
+    auto  &ph  = envGet(LatticeComplex, momphName_);
-    LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
+    auto  &t   = envGet(Lattice<iScalar<vInteger>>, tName_);
-    Gamma                      g(par().gamma);
+    Gamma g(par().gamma);
-    std::vector<Real>          p;
-    Complex                    i(0.0,1.0);
     
-    p  = strToVec<Real>(par().mom);
+    if (!hasPhase_)
-    ph = zero;
-    for(unsigned int mu = 0; mu < env().getNd(); mu++)
     {
-        LatticeCoordinate(coor, mu);
+        Complex           i(0.0,1.0);
-        ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
+        std::vector<Real> p;
+
+        envGetTmp(LatticeComplex, coor);
+        p  = strToVec<Real>(par().mom);
+        ph = zero;
+        for(unsigned int mu = 0; mu < env().getNd(); mu++)
+        {
+            LatticeCoordinate(coor, mu);
+            ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
+        }
+        ph = exp((Real)(2*M_PI)*i*ph);
+        LatticeCoordinate(t, Tp);
+        hasPhase_ = true;
     }
-    ph = exp((Real)(2*M_PI)*i*ph);
-    LatticeCoordinate(t, Tp);
     src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q);
 }
 

extras/Hadrons/Modules/MSource/Wall.hpp

@@ -73,10 +73,14 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
     virtual void execute(void);
+private:
+    bool        hasPhase_{false};
+    std::string momphName_, tName_;
 };
 
 MODULE_REGISTER_NS(Wall, TWall<FIMPL>, MSource);
@@ -88,13 +92,15 @@ MODULE_REGISTER_NS(Wall, TWall<FIMPL>, MSource);
 template <typename FImpl>
 TWall<FImpl>::TWall(const std::string name)
 : Module<WallPar>(name)
+, momphName_ (name + "_momph")
+, tName_ (name + "_t")
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
 template <typename FImpl>
 std::vector<std::string> TWall<FImpl>::getInput(void)
 {
-    std::vector<std::string> in;
+    std::vector<std::string> in = {};
     
     return in;
 }
@@ -111,7 +117,7 @@ std::vector<std::string> TWall<FImpl>::getOutput(void)
 template <typename FImpl>
 void TWall<FImpl>::setup(void)
 {
-    env().template registerLattice<PropagatorField>(getName());
+    envCreateLat(PropagatorField, getName());
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -121,21 +127,28 @@ void TWall<FImpl>::execute(void)
     LOG(Message) << "Generating wall source at t = " << par().tW 
                  << " with momentum " << par().mom << std::endl;
     
-    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    auto  &src = envGet(PropagatorField, getName());
-    Lattice<iScalar<vInteger>> t(env().getGrid());
+    auto  &ph  = envGet(LatticeComplex, momphName_);
-    LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
+    auto  &t   = envGet(Lattice<iScalar<vInteger>>, tName_);
-    std::vector<Real>          p;
-    Complex                    i(0.0,1.0);
     
-    p  = strToVec<Real>(par().mom);
+    if (!hasPhase_)
-    ph = zero;
-    for(unsigned int mu = 0; mu < Nd; mu++)
     {
-        LatticeCoordinate(coor, mu);
+        Complex           i(0.0,1.0);
-        ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
+        std::vector<Real> p;
+
+        envGetTmp(LatticeComplex, coor);
+        p  = strToVec<Real>(par().mom);
+        ph = zero;
+        for(unsigned int mu = 0; mu < env().getNd(); mu++)
+        {
+            LatticeCoordinate(coor, mu);
+            ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
+        }
+        ph = exp((Real)(2*M_PI)*i*ph);
+        LatticeCoordinate(t, Tp);
+        hasPhase_ = true;
     }
-    ph = exp((Real)(2*M_PI)*i*ph);
+
-    LatticeCoordinate(t, Tp);
     src = 1.;
     src = where((t == par().tW), src*ph, 0.*src);
 }

extras/Hadrons/Modules/MSource/Z2.hpp

@@ -76,10 +76,14 @@ public:
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+protected:
     // setup
     virtual void setup(void);
     // execution
     virtual void execute(void);
+private:
+    bool        hasT_{false};
+    std::string tName_;
 };
 
 MODULE_REGISTER_NS(Z2,       TZ2<FIMPL>,        MSource);
@@ -92,6 +96,7 @@ MODULE_REGISTER_NS(ScalarZ2, TZ2<ScalarImplCR>, MSource);
 template <typename FImpl>
 TZ2<FImpl>::TZ2(const std::string name)
 : Module<Z2Par>(name)
+, tName_ (name + "_t")
 {}
 
 // dependencies/products ///////////////////////////////////////////////////////
@@ -115,29 +120,36 @@ std::vector<std::string> TZ2<FImpl>::getOutput(void)
 template <typename FImpl>
 void TZ2<FImpl>::setup(void)
 {
-    env().template registerLattice<PropagatorField>(getName());
+    envCreateLat(PropagatorField, getName());
+    envCacheLat(Lattice<iScalar<vInteger>>, tName_);
+    envTmpLat(LatticeComplex, "eta");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
 void TZ2<FImpl>::execute(void)
 {
-    Lattice<iScalar<vInteger>> t(env().getGrid());
-    LatticeComplex             eta(env().getGrid());
-    Complex                    shift(1., 1.);
-    
     if (par().tA == par().tB)
     {
         LOG(Message) << "Generating Z_2 wall source at t= " << par().tA
-        << std::endl;
+                     << std::endl;
     }
     else
     {
         LOG(Message) << "Generating Z_2 band for " << par().tA << " <= t <= "
-        << par().tB << std::endl;
+                     << par().tB << std::endl;
     }
-    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    
-    LatticeCoordinate(t, Tp);
+    auto    &src = envGet(PropagatorField, getName());
+    auto    &t   = envGet(Lattice<iScalar<vInteger>>, tName_);
+    Complex shift(1., 1.);
+
+    if (!hasT_)
+    {
+        LatticeCoordinate(t, Tp);
+        hasT_ = true;
+    }
+    envGetTmp(LatticeComplex, eta);
     bernoulli(*env().get4dRng(), eta);
     eta = (2.*eta - shift)*(1./::sqrt(2.));
     eta = where((t >= par().tA) and (t <= par().tB), eta, 0.*eta);

extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp

@@ -0,0 +1,185 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MUtilities/TestSeqConserved.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_MUtilities_TestSeqConserved_hpp_
+#define Hadrons_MUtilities_TestSeqConserved_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+  Ward Identity contractions using sequential propagators.
+ -----------------------------
+ 
+ * options:
+ - q:      point source propagator, 5D if available (string)
+ - qSeq:   result of sequential insertion of conserved current using q (string)
+ - action: action used for computation of q (string)
+ - origin: string giving point source origin of q (string)
+ - t_J:    time at which sequential current is inserted (int)
+ - mu:     Lorentz index of current inserted (int)
+ - curr:   current type, e.g. vector/axial (Current)
+*/
+
+/******************************************************************************
+ *                            TestSeqConserved                                *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MUtilities)
+
+class TestSeqConservedPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TestSeqConservedPar,
+                                    std::string,  q,
+                                    std::string,  qSeq,
+                                    std::string,  action,
+                                    std::string,  origin,
+                                    unsigned int, t_J,
+                                    unsigned int, mu,
+                                    Current,      curr);
+};
+
+template <typename FImpl>
+class TTestSeqConserved: public Module<TestSeqConservedPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TTestSeqConserved(const std::string name);
+    // destructor
+    virtual ~TTestSeqConserved(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+protected:
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(TestSeqConserved, TTestSeqConserved<FIMPL>, MUtilities);
+
+/******************************************************************************
+ *                     TTestSeqConserved implementation                       *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TTestSeqConserved<FImpl>::TTestSeqConserved(const std::string name)
+: Module<TestSeqConservedPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TTestSeqConserved<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q, par().qSeq, par().action};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TTestSeqConserved<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TTestSeqConserved<FImpl>::setup(void)
+{
+    auto Ls = env().getObjectLs(par().q);
+    if (Ls != env().getObjectLs(par().action))
+    {
+        HADRON_ERROR(Size, "Ls mismatch between quark action and propagator");
+    }
+    envTmpLat(PropagatorField, "tmp");
+    envTmpLat(LatticeComplex, "c");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TTestSeqConserved<FImpl>::execute(void)
+{
+    // Check sequential insertion of current gives same result as conserved 
+    // current sink upon contraction. Assume q uses a point source.
+
+    auto                  &q    = envGet(PropagatorField, par().q);
+    auto                  &qSeq = envGet(PropagatorField, par().qSeq);
+    auto                  &act  = envGet(FMat, par().action);
+    Gamma                 g5(Gamma::Algebra::Gamma5);
+    Gamma::Algebra        gA = (par().curr == Current::Axial) ?
+                                  Gamma::Algebra::Gamma5 :
+                                  Gamma::Algebra::Identity;
+    Gamma                 g(gA);
+    SitePropagator        qSite;
+    Complex               test_S, test_V, check_S, check_V;
+    std::vector<TComplex> check_buf;
+    std::vector<int>      siteCoord;
+
+    envGetTmp(PropagatorField, tmp);
+    envGetTmp(LatticeComplex, c);
+    siteCoord = strToVec<int>(par().origin);
+    peekSite(qSite, qSeq, siteCoord);
+    test_S = trace(qSite*g);
+    test_V = trace(qSite*g*Gamma::gmu[par().mu]);
+    act.ContractConservedCurrent(q, q, tmp, par().curr, par().mu);
+    c = trace(tmp*g);
+    sliceSum(c, check_buf, Tp);
+    check_S = TensorRemove(check_buf[par().t_J]);
+
+    c = trace(tmp*g*Gamma::gmu[par().mu]);
+    sliceSum(c, check_buf, Tp);
+    check_V = TensorRemove(check_buf[par().t_J]);
+
+    LOG(Message) << "Test S  = " << abs(test_S)   << std::endl;
+    LOG(Message) << "Test V  = " << abs(test_V) << std::endl;
+    LOG(Message) << "Check S = " << abs(check_S) << std::endl;
+    LOG(Message) << "Check V = " << abs(check_V) << std::endl;
+
+    // Check difference = 0
+    check_S -= test_S;
+    check_V -= test_V;
+
+    LOG(Message) << "Consistency check for sequential conserved " 
+                 << par().curr << " current insertion: " << std::endl; 
+    LOG(Message) << "Diff S  = " << abs(check_S) << std::endl;
+    LOG(Message) << "Diff V  = " << abs(check_V) << std::endl;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_TestSeqConserved_hpp_

extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp

@@ -0,0 +1,149 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MUtilities/TestSeqGamma.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_MUtilities_TestSeqGamma_hpp_
+#define Hadrons_MUtilities_TestSeqGamma_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                              TestSeqGamma                                  *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MUtilities)
+
+class TestSeqGammaPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TestSeqGammaPar,
+                                    std::string,    q,
+                                    std::string,    qSeq,
+                                    std::string,    origin,
+                                    Gamma::Algebra, gamma,
+                                    unsigned int,   t_g);
+};
+
+template <typename FImpl>
+class TTestSeqGamma: public Module<TestSeqGammaPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TTestSeqGamma(const std::string name);
+    // destructor
+    virtual ~TTestSeqGamma(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+protected:
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(TestSeqGamma, TTestSeqGamma<FIMPL>, MUtilities);
+
+/******************************************************************************
+ *                      TTestSeqGamma implementation                          *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TTestSeqGamma<FImpl>::TTestSeqGamma(const std::string name)
+: Module<TestSeqGammaPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TTestSeqGamma<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q, par().qSeq};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TTestSeqGamma<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TTestSeqGamma<FImpl>::setup(void)
+{
+    envTmpLat(LatticeComplex, "c");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TTestSeqGamma<FImpl>::execute(void)
+{
+    auto                  &q    = envGet(PropagatorField, par().q);
+    auto                  &qSeq = envGet(PropagatorField, par().qSeq);
+    Gamma                 g5(Gamma::Algebra::Gamma5);
+    Gamma                 g(par().gamma);
+    SitePropagator        qSite;
+    Complex               test, check;
+    std::vector<TComplex> check_buf;
+    std::vector<int>      siteCoord;
+
+    // Check sequential insertion of gamma matrix gives same result as 
+    // insertion of gamma at sink upon contraction. Assume q uses a point 
+    // source.
+    
+    envGetTmp(LatticeComplex, c);
+    siteCoord = strToVec<int>(par().origin);
+    peekSite(qSite, qSeq, siteCoord);
+    test = trace(g*qSite);
+
+    c = trace(adj(g)*g5*adj(q)*g5*g*q);
+    sliceSum(c, check_buf, Tp);
+    check = TensorRemove(check_buf[par().t_g]);
+
+    LOG(Message) << "Seq Result = " << abs(test)  << std::endl;
+    LOG(Message) << "Reference  = " << abs(check) << std::endl;
+
+    // Check difference = 0
+    check -= test;
+
+    LOG(Message) << "Consistency check for sequential " << par().gamma  
+                 << " insertion = " << abs(check) << std::endl;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_TestSeqGamma_hpp_

extras/Hadrons/VirtualMachine.cc

@@ -0,0 +1,622 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/VirtualMachine.cc
+
+Copyright (C) 2017
+
+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/VirtualMachine.hpp>
+#include <Grid/Hadrons/GeneticScheduler.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+/******************************************************************************
+ *                      VirtualMachine implementation                         *
+ ******************************************************************************/
+// trajectory counter //////////////////////////////////////////////////////////
+void VirtualMachine::setTrajectory(const unsigned int traj)
+{
+    traj_ = traj;
+}
+
+unsigned int VirtualMachine::getTrajectory(void) const
+{
+    return traj_;
+}
+
+// module management ///////////////////////////////////////////////////////////
+void VirtualMachine::pushModule(VirtualMachine::ModPt &pt)
+{
+    std::string name = pt->getName();
+    
+    if (!hasModule(name))
+    {
+        std::vector<unsigned int> inputAddress;
+        unsigned int              address;
+        ModuleInfo                m;
+        
+        // module registration -------------------------------------------------
+        m.data = std::move(pt);
+        m.type = typeIdPt(*m.data.get());
+        m.name = name;
+        // input dependencies
+        for (auto &in: m.data->getInput())
+        {
+            if (!env().hasObject(in))
+            {
+                // if object does not exist, add it with no creator module
+                env().addObject(in , -1);
+            }
+            m.input.push_back(env().getObjectAddress(in));
+        }
+        // reference dependencies
+        for (auto &ref: m.data->getReference())
+        {
+            if (!env().hasObject(ref))
+            {
+                // if object does not exist, add it with no creator module
+                env().addObject(ref , -1);
+            }
+            m.input.push_back(env().getObjectAddress(ref));
+        }
+        auto inCopy = m.input;
+        // if module has inputs with references, they need to be added as
+        // an input
+        for (auto &in: inCopy)
+        {
+            int inm = env().getObjectModule(in);
+
+            if (inm > 0)
+            {
+                if (getModule(inm)->getReference().size() > 0)
+                {
+                    for (auto &rin: getModule(inm)->getReference())
+                    {
+                        m.input.push_back(env().getObjectAddress(rin));
+                    }
+                }
+            }
+        }
+        module_.push_back(std::move(m));
+        address              = static_cast<unsigned int>(module_.size() - 1);
+        moduleAddress_[name] = address;
+        // connecting outputs to potential inputs ------------------------------
+        for (auto &out: getModule(address)->getOutput())
+        {
+            if (!env().hasObject(out))
+            {
+                // output does not exists, add it
+                env().addObject(out, address);
+            }
+            else
+            {
+                if (env().getObjectModule(env().getObjectAddress(out)) < 0)
+                {
+                    // output exists but without creator, correct it
+                    env().setObjectModule(env().getObjectAddress(out), address);
+                }
+                else
+                {
+                    // output already fully registered, error
+                    HADRON_ERROR(Definition, "object '" + out
+                                 + "' is already produced by module '"
+                                 + module_[env().getObjectModule(out)].name
+                                 + "' (while pushing module '" + name + "')");
+                }
+                if (getModule(address)->getReference().size() > 0)
+                {
+                    // module has references, dependency should be propagated
+                    // to children modules; find module with `out` as an input
+                    // and add references to their input
+                    auto pred = [this, out](const ModuleInfo &n)
+                    {
+                        auto &in = n.input;
+                        auto it  = std::find(in.begin(), in.end(), 
+                                             env().getObjectAddress(out));
+                        
+                        return (it != in.end());
+                    };
+                    auto it = std::find_if(module_.begin(), module_.end(), pred);
+                    while (it != module_.end())
+                    {
+                        for (auto &ref: getModule(address)->getReference())
+                        {
+                            it->input.push_back(env().getObjectAddress(ref));
+                        }
+                        it = std::find_if(++it, module_.end(), pred);
+                    }   
+                }
+            }
+        }
+        graphOutdated_         = true;
+        memoryProfileOutdated_ = true;
+    }
+    else
+    {
+        HADRON_ERROR(Definition, "module '" + name + "' already exists");
+    }
+}
+
+unsigned int VirtualMachine::getNModule(void) const
+{
+    return module_.size();
+}
+
+void VirtualMachine::createModule(const std::string name, const std::string type,
+                                  XmlReader &reader)
+{
+    auto &factory = ModuleFactory::getInstance();
+    auto pt       = factory.create(type, name);
+    
+    pt->parseParameters(reader, "options");
+    pushModule(pt);
+}
+
+ModuleBase * VirtualMachine::getModule(const unsigned int address) const
+{
+    if (hasModule(address))
+    {
+        return module_[address].data.get();
+    }
+    else
+    {
+        HADRON_ERROR(Definition, "no module with address " + std::to_string(address));
+    }
+}
+
+ModuleBase * VirtualMachine::getModule(const std::string name) const
+{
+    return getModule(getModuleAddress(name));
+}
+
+unsigned int VirtualMachine::getModuleAddress(const std::string name) const
+{
+    if (hasModule(name))
+    {
+        return moduleAddress_.at(name);
+    }
+    else
+    {
+        HADRON_ERROR(Definition, "no module with name '" + name + "'");
+    }
+}
+
+std::string VirtualMachine::getModuleName(const unsigned int address) const
+{
+    if (hasModule(address))
+    {
+        return module_[address].name;
+    }
+    else
+    {
+        HADRON_ERROR(Definition, "no module with address " + std::to_string(address));
+    }
+}
+
+std::string VirtualMachine::getModuleType(const unsigned int address) const
+{
+    if (hasModule(address))
+    {
+        return typeName(module_[address].type);
+    }
+    else
+    {
+        HADRON_ERROR(Definition, "no module with address " + std::to_string(address));
+    }
+}
+
+std::string VirtualMachine::getModuleType(const std::string name) const
+{
+    return getModuleType(getModuleAddress(name));
+}
+
+std::string VirtualMachine::getModuleNamespace(const unsigned int address) const
+{
+    std::string type = getModuleType(address), ns;
+    
+    auto pos2 = type.rfind("::");
+    auto pos1 = type.rfind("::", pos2 - 2);
+    
+    return type.substr(pos1 + 2, pos2 - pos1 - 2);
+}
+
+std::string VirtualMachine::getModuleNamespace(const std::string name) const
+{
+    return getModuleNamespace(getModuleAddress(name));
+}
+
+bool VirtualMachine::hasModule(const unsigned int address) const
+{
+    return (address < module_.size());
+}
+
+bool VirtualMachine::hasModule(const std::string name) const
+{
+    return (moduleAddress_.find(name) != moduleAddress_.end());
+}
+
+// print VM content ////////////////////////////////////////////////////////////
+void VirtualMachine::printContent(void) const
+{
+    LOG(Debug) << "Modules: " << std::endl;
+    for (unsigned int i = 0; i < module_.size(); ++i)
+    {
+        LOG(Debug) << std::setw(4) << i << ": "
+                   << getModuleName(i) << std::endl;
+    }
+}
+
+// module graph ////////////////////////////////////////////////////////////////
+Graph<unsigned int> VirtualMachine::getModuleGraph(void)
+{
+    if (graphOutdated_)
+    {
+        makeModuleGraph();
+        graphOutdated_ = false;
+    }
+
+    return graph_;
+}
+
+void VirtualMachine::makeModuleGraph(void)
+{
+    Graph<unsigned int> graph;
+    
+    // create vertices
+    for (unsigned int m = 0; m < module_.size(); ++m)
+    {
+        graph.addVertex(m);
+    }
+    // create edges
+    for (unsigned int m = 0; m < module_.size(); ++m)
+    {
+        for (auto &in: module_[m].input)
+        {
+            graph.addEdge(env().getObjectModule(in), m);
+        }
+    }
+    graph_ = graph;
+}
+
+// memory profile //////////////////////////////////////////////////////////////
+const VirtualMachine::MemoryProfile & VirtualMachine::getMemoryProfile(void)
+{
+    if (memoryProfileOutdated_)
+    {
+        makeMemoryProfile();
+        memoryProfileOutdated_ = false;
+    }
+
+    return profile_;
+}
+
+void VirtualMachine::makeMemoryProfile(void)
+{
+    bool protect = env().objectsProtected();
+    bool hmsg    = HadronsLogMessage.isActive();
+    bool gmsg    = GridLogMessage.isActive();
+    bool err     = HadronsLogError.isActive();
+    auto program = getModuleGraph().topoSort();
+
+    resetProfile();
+    profile_.module.resize(getNModule());
+    env().protectObjects(false);
+    GridLogMessage.Active(false);
+    HadronsLogMessage.Active(false);
+    HadronsLogError.Active(false);
+    for (auto it = program.rbegin(); it != program.rend(); ++it) 
+    {
+        auto a = *it;
+
+        if (profile_.module[a].empty())
+        {
+            LOG(Debug) << "Profiling memory for module '" << module_[a].name
+                       << "' (" << a << ")..." << std::endl;
+            memoryProfile(a);
+            env().freeAll();
+        }
+    }
+    env().protectObjects(protect);
+    GridLogMessage.Active(gmsg);
+    HadronsLogMessage.Active(hmsg);
+    HadronsLogError.Active(err);
+    LOG(Debug) << "Memory profile:" << std::endl;
+    LOG(Debug) << "----------------" << std::endl;
+    for (unsigned int a = 0; a < profile_.module.size(); ++a)
+    {
+        LOG(Debug) << getModuleName(a) << " (" << a << ")" << std::endl;
+        for (auto &o: profile_.module[a])
+        {
+            LOG(Debug) << "|__ " << env().getObjectName(o.first) << " ("
+                       << sizeString(o.second) << ")" << std::endl;
+        }
+        LOG(Debug) << std::endl;
+    }
+    LOG(Debug) << "----------------" << std::endl;
+}
+
+void VirtualMachine::resetProfile(void)
+{
+    profile_.module.clear();
+    profile_.object.clear();
+}
+
+void VirtualMachine::resizeProfile(void)
+{
+    if (env().getMaxAddress() > profile_.object.size())
+    {
+        MemoryPrint empty;
+
+        empty.size   = 0;
+        empty.module = -1;
+        profile_.object.resize(env().getMaxAddress(), empty);
+    }
+}
+
+void VirtualMachine::updateProfile(const unsigned int address)
+{
+    resizeProfile();
+    for (unsigned int a = 0; a < env().getMaxAddress(); ++a)
+    {
+        if (env().hasCreatedObject(a) and (profile_.object[a].module == -1))
+        {
+            profile_.object[a].size     = env().getObjectSize(a);
+            profile_.object[a].storage  = env().getObjectStorage(a);
+            profile_.object[a].module   = address;
+            profile_.module[address][a] = profile_.object[a].size;
+            if (env().getObjectModule(a) < 0)
+            {
+                env().setObjectModule(a, address);
+            }
+        }
+    }
+}
+
+void VirtualMachine::cleanEnvironment(void)
+{
+    resizeProfile();
+    for (unsigned int a = 0; a < env().getMaxAddress(); ++a)
+    {
+        if (env().hasCreatedObject(a) and (profile_.object[a].module == -1))
+        {
+            env().freeObject(a);
+        }
+    }
+}
+
+void VirtualMachine::memoryProfile(const unsigned int address)
+{
+    auto m = getModule(address);
+
+    LOG(Debug) << "Setting up module '" << m->getName() 
+               << "' (" << address << ")..." << std::endl;
+    try
+    {
+        m->setup();
+        updateProfile(address);
+    }
+    catch (Exceptions::Definition &)
+    {
+        cleanEnvironment();
+        for (auto &in: m->getInput())
+        {
+            memoryProfile(env().getObjectModule(in));
+        }
+        for (auto &ref: m->getReference())
+        {
+            memoryProfile(env().getObjectModule(ref));
+        }
+        m->setup();
+        updateProfile(address);
+    }
+}
+
+void VirtualMachine::memoryProfile(const std::string name)
+{
+    memoryProfile(getModuleAddress(name));
+}
+
+// garbage collector ///////////////////////////////////////////////////////////
+VirtualMachine::GarbageSchedule 
+VirtualMachine::makeGarbageSchedule(const Program &p) const
+{
+    GarbageSchedule freeProg;
+    
+    freeProg.resize(p.size());
+    for (unsigned int a = 0; a < env().getMaxAddress(); ++a)
+    {
+        if (env().getObjectStorage(a) == Environment::Storage::temporary)
+        {
+            auto it = std::find(p.begin(), p.end(), env().getObjectModule(a));
+
+            if (it != p.end())
+            {
+                freeProg[std::distance(p.begin(), it)].insert(a);
+            }
+        }
+        else if (env().getObjectStorage(a) == Environment::Storage::object)
+        {
+            auto pred = [a, this](const unsigned int b)
+            {
+                auto &in = module_[b].input;
+                auto it  = std::find(in.begin(), in.end(), a);
+                
+                return (it != in.end()) or (b == env().getObjectModule(a));
+            };
+            auto it = std::find_if(p.rbegin(), p.rend(), pred);
+            if (it != p.rend())
+            {
+                freeProg[std::distance(it, p.rend()) - 1].insert(a);
+            }
+        }
+    }
+
+    return freeProg;
+}
+
+// high-water memory function //////////////////////////////////////////////////
+VirtualMachine::Size VirtualMachine::memoryNeeded(const Program &p)
+{
+    const MemoryProfile &profile = getMemoryProfile();
+    GarbageSchedule     freep    = makeGarbageSchedule(p);
+    Size                current = 0, max = 0;
+
+    for (unsigned int i = 0; i < p.size(); ++i)
+    {
+        for (auto &o: profile.module[p[i]])
+        {
+            current += o.second;
+        }
+        max = std::max(current, max);
+        for (auto &o: freep[i])
+        {
+            current -= profile.object[o].size;
+        }
+    }
+
+    return max;
+}
+
+// genetic scheduler ///////////////////////////////////////////////////////////
+VirtualMachine::Program VirtualMachine::schedule(const GeneticPar &par)
+{
+    typedef GeneticScheduler<Size, unsigned int> Scheduler;
+
+    auto graph = getModuleGraph();
+
+    //constrained topological sort using a genetic algorithm
+    LOG(Message) << "Scheduling computation..." << std::endl;
+    LOG(Message) << "               #module= " << graph.size() << std::endl;
+    LOG(Message) << "       population size= " << par.popSize << std::endl;
+    LOG(Message) << "       max. generation= " << par.maxGen << std::endl;
+    LOG(Message) << "  max. cst. generation= " << par.maxCstGen << std::endl;
+    LOG(Message) << "         mutation rate= " << par.mutationRate << std::endl;
+    
+    unsigned int          k = 0, gen, prevPeak, nCstPeak = 0;
+    std::random_device    rd;
+    Scheduler::Parameters gpar;
+    
+    gpar.popSize      = par.popSize;
+    gpar.mutationRate = par.mutationRate;
+    gpar.seed         = rd();
+    CartesianCommunicator::BroadcastWorld(0, &(gpar.seed), sizeof(gpar.seed));
+    Scheduler::ObjFunc memPeak = [this](const Program &p)->Size
+    {
+        return memoryNeeded(p);
+    };
+    Scheduler scheduler(graph, memPeak, gpar);
+    gen = 0;
+    do
+    {
+        LOG(Debug) << "Generation " << gen << ":" << std::endl;
+        scheduler.nextGeneration();
+        if (gen != 0)
+        {
+            if (prevPeak == scheduler.getMinValue())
+            {
+                nCstPeak++;
+            }
+            else
+            {
+                nCstPeak = 0;
+            }
+        }
+        
+        prevPeak = scheduler.getMinValue();
+        if (gen % 10 == 0)
+        {
+            LOG(Iterative) << "Generation " << gen << ": "
+                           << sizeString(scheduler.getMinValue()) << std::endl;
+        }
+        
+        gen++;
+    } while ((gen < par.maxGen) and (nCstPeak < par.maxCstGen));
+    
+    return scheduler.getMinSchedule();
+}
+
+// general execution ///////////////////////////////////////////////////////////
+#define BIG_SEP "==============="
+#define SEP     "---------------"
+#define MEM_MSG(size) sizeString(size)
+
+void VirtualMachine::executeProgram(const Program &p) const
+{
+    Size            memPeak = 0, sizeBefore, sizeAfter;
+    GarbageSchedule freeProg;
+    
+    // build garbage collection schedule
+    LOG(Debug) << "Building garbage collection schedule..." << std::endl;
+    freeProg = makeGarbageSchedule(p);
+
+    // program execution
+    LOG(Debug) << "Executing program..." << std::endl;
+    for (unsigned int i = 0; i < p.size(); ++i)
+    {
+        // execute module
+        LOG(Message) << SEP << " Measurement step " << i + 1 << "/"
+                     << p.size() << " (module '" << module_[p[i]].name
+                     << "') " << SEP << std::endl;
+        (*module_[p[i]].data)();
+        sizeBefore = env().getTotalSize();
+        // print used memory after execution
+        LOG(Message) << "Allocated objects: " << MEM_MSG(sizeBefore)
+                     << std::endl;
+        if (sizeBefore > memPeak)
+        {
+            memPeak = sizeBefore;
+        }
+        // garbage collection for step i
+        LOG(Message) << "Garbage collection..." << std::endl;
+        for (auto &j: freeProg[i])
+        {
+            env().freeObject(j);
+        }
+        // print used memory after garbage collection if necessary
+        sizeAfter = env().getTotalSize();
+        if (sizeBefore != sizeAfter)
+        {
+            LOG(Message) << "Allocated objects: " << MEM_MSG(sizeAfter)
+                            << std::endl;
+        }
+        else
+        {
+            LOG(Message) << "Nothing to free" << std::endl;
+        }
+    }
+}
+
+void VirtualMachine::executeProgram(const std::vector<std::string> &p) const
+{
+    Program pAddress;
+    
+    for (auto &n: p)
+    {
+        pAddress.push_back(getModuleAddress(n));
+    }
+    executeProgram(pAddress);
+}

extras/Hadrons/VirtualMachine.hpp

@@ -0,0 +1,207 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/VirtualMachine.hpp
+
+Copyright (C) 2017
+
+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_VirtualMachine_hpp_
+#define Hadrons_VirtualMachine_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Graph.hpp>
+#include <Grid/Hadrons/Environment.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+#define DEFINE_VM_ALIAS \
+inline VirtualMachine & vm(void) const\
+{\
+    return VirtualMachine::getInstance();\
+}
+
+/******************************************************************************
+ *                   Virtual machine for module execution                     *
+ ******************************************************************************/
+// forward declaration of Module
+class ModuleBase;
+
+class VirtualMachine
+{
+    SINGLETON_DEFCTOR(VirtualMachine);
+public:
+    typedef SITE_SIZE_TYPE                      Size;
+    typedef std::unique_ptr<ModuleBase>         ModPt;
+    typedef std::vector<std::set<unsigned int>> GarbageSchedule;
+    typedef std::vector<unsigned int>           Program;
+    struct MemoryPrint
+    {
+        Size                 size;
+        Environment::Storage storage;
+        int                  module;
+    };
+    struct MemoryProfile
+    {
+        std::vector<std::map<unsigned int, Size>> module;
+        std::vector<MemoryPrint>                  object;
+    };
+    class GeneticPar: Serializable
+    {
+    public:
+        GeneticPar(void):
+            popSize{20}, maxGen{1000}, maxCstGen{100}, mutationRate{.1} {};
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(GeneticPar,
+                                        unsigned int, popSize,
+                                        unsigned int, maxGen,
+                                        unsigned int, maxCstGen,
+                                        double      , mutationRate);
+    };
+private:
+    struct ModuleInfo
+    {
+        const std::type_info      *type{nullptr};
+        std::string               name;
+        ModPt                     data{nullptr};
+        std::vector<unsigned int> input;
+        size_t                    maxAllocated;
+    };
+public:
+    // trajectory counter
+    void                setTrajectory(const unsigned int traj);
+    unsigned int        getTrajectory(void) const;
+    // module management
+    void                pushModule(ModPt &pt);
+    template <typename M>
+    void                createModule(const std::string name);
+    template <typename M>
+    void                createModule(const std::string name,
+                                         const typename M::Par &par);
+    void                createModule(const std::string name,
+                                         const std::string type,
+                                         XmlReader &reader);
+    unsigned int        getNModule(void) const;
+    ModuleBase *        getModule(const unsigned int address) const;
+    ModuleBase *        getModule(const std::string name) const;
+    template <typename M>
+    M *                 getModule(const unsigned int address) const;
+    template <typename M>
+    M *                 getModule(const std::string name) const;
+    unsigned int        getModuleAddress(const std::string name) const;
+    std::string         getModuleName(const unsigned int address) const;
+    std::string         getModuleType(const unsigned int address) const;
+    std::string         getModuleType(const std::string name) const;
+    std::string         getModuleNamespace(const unsigned int address) const;
+    std::string         getModuleNamespace(const std::string name) const;
+    bool                hasModule(const unsigned int address) const;
+    bool                hasModule(const std::string name) const;
+    // print VM content
+    void                printContent(void) const;
+    // module graph (could be a const reference if topoSort was const)
+    Graph<unsigned int> getModuleGraph(void);
+    // memory profile
+    const MemoryProfile &getMemoryProfile(void);
+    // garbage collector
+    GarbageSchedule     makeGarbageSchedule(const Program &p) const;
+    // high-water memory function
+    Size                memoryNeeded(const Program &p);
+    // genetic scheduler
+    Program             schedule(const GeneticPar &par);
+    // general execution
+    void                executeProgram(const Program &p) const;
+    void                executeProgram(const std::vector<std::string> &p) const;
+private:
+    // environment shortcut
+    DEFINE_ENV_ALIAS;
+    // module graph
+    void makeModuleGraph(void);
+    // memory profile
+    void makeMemoryProfile(void);
+    void resetProfile(void);
+    void resizeProfile(void);
+    void updateProfile(const unsigned int address);
+    void cleanEnvironment(void);
+    void memoryProfile(const std::string name);
+    void memoryProfile(const unsigned int address);
+private:
+    // general
+    unsigned int                        traj_;
+    // module and related maps
+    std::vector<ModuleInfo>             module_;
+    std::map<std::string, unsigned int> moduleAddress_;
+    std::string                         currentModule_{""};
+    // module graph
+    bool                                graphOutdated_{true};
+    Graph<unsigned int>                 graph_;
+    // memory profile
+    bool                                memoryProfileOutdated_{true};
+    MemoryProfile                       profile_;
+};
+
+/******************************************************************************
+ *                   VirtualMachine template implementation                   *
+ ******************************************************************************/
+// module management ///////////////////////////////////////////////////////////
+template <typename M>
+void VirtualMachine::createModule(const std::string name)
+{
+    ModPt pt(new M(name));
+    
+    pushModule(pt);
+}
+
+template <typename M>
+void VirtualMachine::createModule(const std::string name,
+                               const typename M::Par &par)
+{
+    ModPt pt(new M(name));
+    
+    static_cast<M *>(pt.get())->setPar(par);
+    pushModule(pt);
+}
+
+template <typename M>
+M * VirtualMachine::getModule(const unsigned int address) const
+{
+    if (auto *pt = dynamic_cast<M *>(getModule(address)))
+    {
+        return pt;
+    }
+    else
+    {
+        HADRON_ERROR(Definition, "module '" + module_[address].name
+                     + "' does not have type " + typeid(M).name()
+                     + "(has type: " + getModuleType(address) + ")");
+    }
+}
+
+template <typename M>
+M * VirtualMachine::getModule(const std::string name) const
+{
+    return getModule<M>(getModuleAddress(name));
+}
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_VirtualMachine_hpp_

extras/Hadrons/make_module_list.sh

@@ -7,6 +7,36 @@ echo 'modules_hpp =\' >> modules.inc
 find Modules -name '*.hpp' -type f -print | sed 's/^/  /;$q;s/$/ \\/' >> modules.inc
 echo '' >> modules.inc
 rm -f Modules.hpp
+echo "/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules.hpp
+
+Copyright (C) 2015
+Copyright (C) 2016
+Copyright (C) 2017
+
+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 */
+" > Modules.hpp
 for f in `find Modules -name '*.hpp'`; do
 	echo "#include <Grid/Hadrons/${f}>" >> Modules.hpp
 done

extras/Hadrons/modules.inc

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

lib/allocator/AlignedAllocator.cc

@@ -3,9 +3,12 @@
 
 namespace Grid {
 
+MemoryStats *MemoryProfiler::stats = nullptr;
+bool         MemoryProfiler::debug = false;
+
 int PointerCache::victim;
 
-  PointerCache::PointerCacheEntry PointerCache::Entries[PointerCache::Ncache];
+PointerCache::PointerCacheEntry PointerCache::Entries[PointerCache::Ncache];
 
 void *PointerCache::Insert(void *ptr,size_t bytes) {
 
@@ -94,4 +97,29 @@ void check_huge_pages(void *Buf,uint64_t BYTES)
 #endif
 }
 
+std::string sizeString(const size_t bytes)
+{
+  constexpr unsigned int bufSize = 256;
+  const char             *suffixes[7] = {"", "K", "M", "G", "T", "P", "E"};
+  char                   buf[256];
+  size_t                 s     = 0;
+  double                 count = bytes;
+  
+  while (count >= 1024 && s < 7)
+  {
+      s++;
+      count /= 1024;
+  }
+  if (count - floor(count) == 0.0)
+  {
+      snprintf(buf, bufSize, "%d %sB", (int)count, suffixes[s]);
+  }
+  else
+  {
+      snprintf(buf, bufSize, "%.1f %sB", count, suffixes[s]);
+  }
+  
+  return std::string(buf);
+}
+
 }

lib/allocator/AlignedAllocator.h

@@ -63,6 +63,64 @@ namespace Grid {
     static void *Lookup(size_t bytes) ;
 
   };
+  
+  std::string sizeString(size_t bytes);
+
+  struct MemoryStats
+  {
+    size_t totalAllocated{0}, maxAllocated{0}, 
+           currentlyAllocated{0}, totalFreed{0};
+  };
+    
+  class MemoryProfiler
+  {
+  public:
+    static MemoryStats *stats;
+    static bool        debug;
+  };
+
+  #define memString(bytes) std::to_string(bytes) + " (" + sizeString(bytes) + ")"
+  #define profilerDebugPrint \
+  if (MemoryProfiler::stats)\
+  {\
+    auto s = MemoryProfiler::stats;\
+    std::cout << GridLogDebug << "[Memory debug] Stats " << MemoryProfiler::stats << std::endl;\
+    std::cout << GridLogDebug << "[Memory debug] total  : " << memString(s->totalAllocated) \
+              << std::endl;\
+    std::cout << GridLogDebug << "[Memory debug] max    : " << memString(s->maxAllocated) \
+              << std::endl;\
+    std::cout << GridLogDebug << "[Memory debug] current: " << memString(s->currentlyAllocated) \
+              << std::endl;\
+    std::cout << GridLogDebug << "[Memory debug] freed  : " << memString(s->totalFreed) \
+              << std::endl;\
+  }
+
+  #define profilerAllocate(bytes)\
+  if (MemoryProfiler::stats)\
+  {\
+    auto s = MemoryProfiler::stats;\
+    s->totalAllocated     += (bytes);\
+    s->currentlyAllocated += (bytes);\
+    s->maxAllocated        = std::max(s->maxAllocated, s->currentlyAllocated);\
+  }\
+  if (MemoryProfiler::debug)\
+  {\
+    std::cout << GridLogDebug << "[Memory debug] allocating " << memString(bytes) << std::endl;\
+    profilerDebugPrint;\
+  }
+
+  #define profilerFree(bytes)\
+  if (MemoryProfiler::stats)\
+  {\
+    auto s = MemoryProfiler::stats;\
+    s->totalFreed         += (bytes);\
+    s->currentlyAllocated -= (bytes);\
+  }\
+  if (MemoryProfiler::debug)\
+  {\
+    std::cout << GridLogDebug << "[Memory debug] freeing " << memString(bytes) << std::endl;\
+    profilerDebugPrint;\
+  }
 
   void check_huge_pages(void *Buf,uint64_t BYTES);
 
@@ -92,6 +150,7 @@ public:
   pointer allocate(size_type __n, const void* _p= 0)
   { 
     size_type bytes = __n*sizeof(_Tp);
+    profilerAllocate(bytes);
 
     _Tp *ptr = (_Tp *) PointerCache::Lookup(bytes);
     //    if ( ptr != NULL ) 
@@ -122,6 +181,8 @@ public:
   void deallocate(pointer __p, size_type __n) { 
     size_type bytes = __n * sizeof(_Tp);
 
+    profilerFree(bytes);
+
     pointer __freeme = (pointer)PointerCache::Insert((void *)__p,bytes);
 
 #ifdef HAVE_MM_MALLOC_H
@@ -172,10 +233,13 @@ public:
 #ifdef GRID_COMMS_SHMEM
   pointer allocate(size_type __n, const void* _p= 0)
   {
+    size_type bytes = __n*sizeof(_Tp);
+
+    profilerAllocate(bytes);
 #ifdef CRAY
-    _Tp *ptr = (_Tp *) shmem_align(__n*sizeof(_Tp),64);
+    _Tp *ptr = (_Tp *) shmem_align(bytes,64);
 #else
-    _Tp *ptr = (_Tp *) shmem_align(64,__n*sizeof(_Tp));
+    _Tp *ptr = (_Tp *) shmem_align(64,bytes);
 #endif
 #ifdef PARANOID_SYMMETRIC_HEAP
     static void * bcast;
@@ -193,18 +257,23 @@ public:
 #endif 
     return ptr;
   }
-  void deallocate(pointer __p, size_type) { 
+  void deallocate(pointer __p, size_type __n) { 
+    size_type bytes = __n*sizeof(_Tp);
+
+    profilerFree(bytes);
     shmem_free((void *)__p);
   }
 #else
   pointer allocate(size_type __n, const void* _p= 0) 
   {
-#ifdef HAVE_MM_MALLOC_H
-    _Tp * ptr = (_Tp *) _mm_malloc(__n*sizeof(_Tp),GRID_ALLOC_ALIGN);
-#else
-    _Tp * ptr = (_Tp *) memalign(GRID_ALLOC_ALIGN,__n*sizeof(_Tp));
-#endif
     size_type bytes = __n*sizeof(_Tp);
+    
+    profilerAllocate(bytes);
+#ifdef HAVE_MM_MALLOC_H
+    _Tp * ptr = (_Tp *) _mm_malloc(bytes, GRID_ALLOC_ALIGN);
+#else
+    _Tp * ptr = (_Tp *) memalign(GRID_ALLOC_ALIGN, bytes);
+#endif
     uint8_t *cp = (uint8_t *)ptr;
     if ( ptr ) { 
     // One touch per 4k page, static OMP loop to catch same loop order
@@ -215,7 +284,10 @@ public:
     }
     return ptr;
   }
-  void deallocate(pointer __p, size_type) { 
+  void deallocate(pointer __p, size_type __n) {
+    size_type bytes = __n*sizeof(_Tp);
+
+    profilerFree(bytes);
 #ifdef HAVE_MM_MALLOC_H
     _mm_free((void *)__p); 
 #else

lib/qcd/QCD.h

@@ -492,6 +492,14 @@ namespace QCD {
       return traceIndex<ColourIndex>(lhs);
     }
 
+    //////////////////////////////////////////
+    // Current types
+    //////////////////////////////////////////
+    GRID_SERIALIZABLE_ENUM(Current, undef,
+                           Vector,  0,
+                           Axial,   1,
+                           Tadpole, 2);
+
 }   //namespace QCD
 } // Grid
 

lib/qcd/action/fermion/FermionOperator.h

@@ -47,6 +47,7 @@ namespace Grid {
       INHERIT_IMPL_TYPES(Impl);
 
       FermionOperator(const ImplParams &p= ImplParams()) : Impl(p) {};
+      virtual ~FermionOperator(void) = default;
 
       virtual FermionField &tmp(void) = 0;
 
@@ -112,6 +113,21 @@ namespace Grid {
       ///////////////////////////////////////////////
       virtual void ImportGauge(const GaugeField & _U)=0;
 
+      //////////////////////////////////////////////////////////////////////
+      // Conserved currents, either contract at sink or insert sequentially.
+      //////////////////////////////////////////////////////////////////////
+      virtual void ContractConservedCurrent(PropagatorField &q_in_1,
+                                            PropagatorField &q_in_2,
+                                            PropagatorField &q_out,
+                                            Current curr_type,
+                                            unsigned int mu)=0;
+      virtual void SeqConservedCurrent(PropagatorField &q_in, 
+                                       PropagatorField &q_out,
+                                       Current curr_type,
+                                       unsigned int mu,
+                                       std::vector<Real> mom,
+                                       unsigned int tmin, 
+                                       unsigned int tmax)=0;
     };
 
   }

lib/qcd/action/fermion/FermionOperatorImpl.h

@@ -212,6 +212,13 @@ namespace QCD {
                          StencilImpl &St) {
       mult(&phi(), &U(mu), &chi());
     }
+    
+    inline void multLinkProp(SitePropagator &phi,
+                             const SiteDoubledGaugeField &U,
+                             const SitePropagator &chi,
+                             int mu) {
+       mult(&phi(), &U(mu), &chi());
+    }
       
     template <class ref>
     inline void loadLinkElement(Simd &reg, ref &memory) {
@@ -340,7 +347,20 @@ class DomainWallVec5dImpl :  public PeriodicGaugeImpl< GaugeImplTypes< S,Nrepres
     }
     mult(&phi(), &UU(), &chi());
   }
-      
+
+  inline void multLinkProp(SitePropagator &phi,
+                           const SiteDoubledGaugeField &U,
+                           const SitePropagator &chi,
+                           int mu) {
+    SiteGaugeLink UU;
+    for (int i = 0; i < Nrepresentation; i++) {
+      for (int j = 0; j < Nrepresentation; j++) {
+        vsplat(UU()()(i, j), U(mu)()(i, j));
+      }
+    }
+    mult(&phi(), &UU(), &chi());
+  }
+
   inline void DoubleStore(GridBase *GaugeGrid, DoubledGaugeField &Uds,const GaugeField &Umu) 
   {
     SiteScalarGaugeField  ScalarUmu;
@@ -537,7 +557,12 @@ class GparityWilsonImpl : public ConjugateGaugeImpl<GaugeImplTypes<S, Nrepresent
    }
    
  }
-
+    // Fixme: Gparity prop * link
+    inline void multLinkProp(SitePropagator &phi, const SiteDoubledGaugeField &U,
+                             const SitePropagator &chi, int mu)
+    {
+        assert(0);
+    }
 
  template <class ref>
  inline void loadLinkElement(Simd &reg, ref &memory) {

lib/qcd/action/fermion/ImprovedStaggeredFermion.cc

@@ -393,6 +393,31 @@ void ImprovedStaggeredFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder
   }
 };
 
+//////////////////////////////////////////////////////// 
+// Conserved current - not yet implemented.
+////////////////////////////////////////////////////////
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
+                                                        PropagatorField &q_in_2,
+                                                        PropagatorField &q_out,
+                                                        Current curr_type,
+                                                        unsigned int mu)
+{
+    assert(0);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
+                                                         PropagatorField &q_out,
+                                                         Current curr_type,
+                                                         unsigned int mu, 
+                                                         std::vector<Real> mom,
+                                                         unsigned int tmin,
+                                                         unsigned int tmax)
+{
+    assert(0);
+}
+
 FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion);
 
   //AdjointFermOpTemplateInstantiate(ImprovedStaggeredFermion);

lib/qcd/action/fermion/ImprovedStaggeredFermion.h

@@ -157,6 +157,22 @@ class ImprovedStaggeredFermion : public StaggeredKernels<Impl>, public ImprovedS
 
   LebesgueOrder Lebesgue;
   LebesgueOrder LebesgueEvenOdd;
+  
+  ///////////////////////////////////////////////////////////////
+  // Conserved current utilities
+  ///////////////////////////////////////////////////////////////
+  void ContractConservedCurrent(PropagatorField &q_in_1,
+                                PropagatorField &q_in_2,
+                                PropagatorField &q_out,
+                                Current curr_type,
+                                unsigned int mu);
+  void SeqConservedCurrent(PropagatorField &q_in,
+                           PropagatorField &q_out,
+                           Current curr_type,
+                           unsigned int mu, 
+                           std::vector<Real> mom,
+                           unsigned int tmin,
+                           unsigned int tmax);
 };
 
 typedef ImprovedStaggeredFermion<StaggeredImplF> ImprovedStaggeredFermionF;

lib/qcd/action/fermion/ImprovedStaggeredFermion5D.cc

@@ -405,6 +405,30 @@ void ImprovedStaggeredFermion5D<Impl>::MooeeInvDag(const FermionField &in,
   MooeeInv(in, out);
 }
 
+//////////////////////////////////////////////////////// 
+// Conserved current - not yet implemented.
+////////////////////////////////////////////////////////
+template <class Impl>
+void ImprovedStaggeredFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
+                                                         PropagatorField &q_in_2,
+                                                         PropagatorField &q_out,
+                                                         Current curr_type,
+                                                         unsigned int mu)
+{
+    assert(0);
+}
+
+template <class Impl>
+void ImprovedStaggeredFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
+                                                          PropagatorField &q_out,
+                                                          Current curr_type,
+                                                          unsigned int mu, 
+                                                          std::vector<Real> mom,
+                                                          unsigned int tmin,
+                                                          unsigned int tmax)
+{
+    assert(0);
+}
 
 FermOpStaggeredTemplateInstantiate(ImprovedStaggeredFermion5D);
 FermOpStaggeredVec5dTemplateInstantiate(ImprovedStaggeredFermion5D);

lib/qcd/action/fermion/ImprovedStaggeredFermion5D.h

@@ -170,6 +170,21 @@ namespace QCD {
     // Comms buffer
     std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  comm_buf;
     
+    ///////////////////////////////////////////////////////////////
+    // Conserved current utilities
+    ///////////////////////////////////////////////////////////////
+    void ContractConservedCurrent(PropagatorField &q_in_1,
+                                  PropagatorField &q_in_2,
+                                  PropagatorField &q_out,
+                                  Current curr_type,
+                                  unsigned int mu);
+    void SeqConservedCurrent(PropagatorField &q_in,
+                             PropagatorField &q_out,
+                             Current curr_type,
+                             unsigned int mu, 
+                             std::vector<Real> mom,
+                             unsigned int tmin,
+                             unsigned int tmax);
   };
 
 }}

lib/qcd/action/fermion/WilsonFermion.cc

@@ -345,6 +345,112 @@ void WilsonFermion<Impl>::DhopInternal(StencilImpl &st, LebesgueOrder &lo,
   }
 };
 
+/*******************************************************************************
+ * Conserved current utilities for Wilson fermions, for contracting propagators
+ * to make a conserved current sink or inserting the conserved current 
+ * sequentially.
+ ******************************************************************************/
+template <class Impl>
+void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
+                                                   PropagatorField &q_in_2,
+                                                   PropagatorField &q_out,
+                                                   Current curr_type,
+                                                   unsigned int mu)
+{
+    Gamma g5(Gamma::Algebra::Gamma5);
+    conformable(_grid, q_in_1._grid);
+    conformable(_grid, q_in_2._grid);
+    conformable(_grid, q_out._grid);
+    PropagatorField tmp1(_grid), tmp2(_grid);
+    q_out = zero;
+
+    // Forward, need q1(x + mu), q2(x). Backward, need q1(x), q2(x + mu).
+    // Inefficient comms method but not performance critical.
+    tmp1 = Cshift(q_in_1, mu, 1);
+    tmp2 = Cshift(q_in_2, mu, 1);
+    parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
+    {
+        Kernels::ContractConservedCurrentSiteFwd(tmp1._odata[sU],
+                                                 q_in_2._odata[sU],
+                                                 q_out._odata[sU],
+                                                 Umu, sU, mu);
+        Kernels::ContractConservedCurrentSiteBwd(q_in_1._odata[sU],
+                                                 tmp2._odata[sU],
+                                                 q_out._odata[sU],
+                                                 Umu, sU, mu);
+    }
+}
+
+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)
+{
+    conformable(_grid, q_in._grid);
+    conformable(_grid, q_out._grid);
+    Lattice<iSinglet<Simd>> ph(_grid), coor(_grid);
+    Complex 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((Real)(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;
+    tmp = ph*q_in;
+    tmpBwd = Cshift(tmp, mu, -1);
+
+    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)
+        {
+            Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sU], 
+                                                q_out._odata[sU], 
+                                                Umu, sU, mu, t_mask);
+        }
+
+        // Repeat for backward direction.
+        t_mask     = ((coords._odata[sU] >= (tmin + tshift)) && 
+                      (coords._odata[sU] <= (tmax + tshift)));
+
+	//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 ));
+
+        timeSlices = Reduce(t_mask);
+
+        if (timeSlices > 0)
+        {
+            Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sU], 
+                                                q_out._odata[sU], 
+                                                Umu, sU, mu, t_mask);
+        }
+    }
+}
+
 FermOpTemplateInstantiate(WilsonFermion);
 AdjointFermOpTemplateInstantiate(WilsonFermion);
 TwoIndexFermOpTemplateInstantiate(WilsonFermion);

lib/qcd/action/fermion/WilsonFermion.h

@@ -146,6 +146,22 @@ class WilsonFermion : public WilsonKernels<Impl>, public WilsonFermionStatic {
 
   LebesgueOrder Lebesgue;
   LebesgueOrder LebesgueEvenOdd;
+  
+  ///////////////////////////////////////////////////////////////
+  // Conserved current utilities
+  ///////////////////////////////////////////////////////////////
+  void ContractConservedCurrent(PropagatorField &q_in_1,
+                                PropagatorField &q_in_2,
+                                PropagatorField &q_out,
+                                Current curr_type,
+                                unsigned int mu);
+  void SeqConservedCurrent(PropagatorField &q_in,
+                           PropagatorField &q_out,
+                           Current curr_type,
+                           unsigned int mu, 
+                           std::vector<Real> mom,
+                           unsigned int tmin,
+                           unsigned int tmax);
 };
 
 typedef WilsonFermion<WilsonImplF> WilsonFermionF;

lib/qcd/action/fermion/WilsonFermion5D.cc

@@ -12,6 +12,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
 Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
 Author: paboyle <paboyle@ph.ed.ac.uk>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: Andrew Lawson <andrew.lawson1991@gmail.com>
 
     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
@@ -702,6 +703,168 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const Fe
 
 }
 
+/*******************************************************************************
+ * Conserved current utilities for Wilson fermions, for contracting propagators
+ * to make a conserved current sink or inserting the conserved current 
+ * sequentially.
+ ******************************************************************************/
+
+// Helper macro to reverse Simd vector. Fixme: slow, generic implementation.
+#define REVERSE_LS(qSite, qSiteRev, Nsimd) \
+{ \
+    std::vector<typename SitePropagator::scalar_object> qSiteVec(Nsimd); \
+    extract(qSite, qSiteVec); \
+    for (int i = 0; i < Nsimd / 2; ++i) \
+    { \
+        typename SitePropagator::scalar_object tmp = qSiteVec[i]; \
+        qSiteVec[i] = qSiteVec[Nsimd - i - 1]; \
+        qSiteVec[Nsimd - i - 1] = tmp; \
+    } \
+    merge(qSiteRev, qSiteVec); \
+}
+
+template <class Impl>
+void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
+                                                     PropagatorField &q_in_2,
+                                                     PropagatorField &q_out,
+                                                     Current curr_type,
+                                                     unsigned int mu)
+{
+    conformable(q_in_1._grid, FermionGrid());
+    conformable(q_in_1._grid, q_in_2._grid);
+    conformable(_FourDimGrid, q_out._grid);
+    PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid());
+    unsigned int LLs = q_in_1._grid->_rdimensions[0];
+    q_out = zero;
+
+    // Forward, need q1(x + mu, s), q2(x, Ls - 1 - s). Backward, need q1(x, s), 
+    // q2(x + mu, Ls - 1 - s). 5D lattice so shift 4D coordinate mu by one.
+    tmp1 = Cshift(q_in_1, mu + 1, 1);
+    tmp2 = Cshift(q_in_2, mu + 1, 1);
+    parallel_for (unsigned int sU = 0; sU < Umu._grid->oSites(); ++sU)
+    {
+        unsigned int sF1 = sU * LLs;
+        unsigned int sF2 = (sU + 1) * LLs - 1;
+
+        for (unsigned int s = 0; s < LLs; ++s)
+        {
+            bool axial_sign = ((curr_type == Current::Axial) && \
+                               (s < (LLs / 2)));
+            SitePropagator qSite2, qmuSite2;
+
+            // If vectorised in 5th dimension, reverse q2 vector to match up
+            // sites correctly.
+            if (Impl::LsVectorised)
+            {
+                REVERSE_LS(q_in_2._odata[sF2], qSite2, Ls / LLs);
+                REVERSE_LS(tmp2._odata[sF2], qmuSite2, Ls / LLs);
+            }
+            else
+            {
+                qSite2   = q_in_2._odata[sF2];
+                qmuSite2 = tmp2._odata[sF2];
+            }
+            Kernels::ContractConservedCurrentSiteFwd(tmp1._odata[sF1], 
+                                                     qSite2, 
+                                                     q_out._odata[sU],
+                                                     Umu, sU, mu, axial_sign);
+            Kernels::ContractConservedCurrentSiteBwd(q_in_1._odata[sF1],
+                                                     qmuSite2,
+                                                     q_out._odata[sU],
+                                                     Umu, sU, mu, axial_sign);
+            sF1++;
+            sF2--;
+        }
+    }
+}
+
+
+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)
+{
+    conformable(q_in._grid, FermionGrid());
+    conformable(q_in._grid, q_out._grid);
+    Lattice<iSinglet<Simd>> ph(FermionGrid()), coor(FermionGrid());
+    PropagatorField tmpFwd(FermionGrid()), tmpBwd(FermionGrid()),
+                    tmp(FermionGrid());
+    Complex 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((Real)(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)
+    {
+        // 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)
+        {
+            unsigned int sF = sU * LLs;
+            for (unsigned int s = 0; s < LLs; ++s)
+            {
+                bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
+                Kernels::SeqConservedCurrentSiteFwd(tmpFwd._odata[sF], 
+                                                    q_out._odata[sF], Umu, sU,
+                                                    mu, t_mask, axial_sign);
+                ++sF;
+            }
+        }
+
+        // Repeat for backward direction.
+        t_mask     = ((coords._odata[sU] >= (tmin + tshift)) && 
+                      (coords._odata[sU] <= (tmax + tshift)));
+
+	//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 ));
+
+        timeSlices = Reduce(t_mask);
+
+        if (timeSlices > 0)
+        {
+            unsigned int sF = sU * LLs;
+            for (unsigned int s = 0; s < LLs; ++s)
+            {
+                bool axial_sign = ((curr_type == Current::Axial) && (s < (LLs / 2)));
+                Kernels::SeqConservedCurrentSiteBwd(tmpBwd._odata[sF], 
+                                                    q_out._odata[sF], Umu, sU,
+                                                    mu, t_mask, axial_sign);
+                ++sF;
+            }
+        }
+    }
+}
+
 FermOpTemplateInstantiate(WilsonFermion5D);
 GparityFermOpTemplateInstantiate(WilsonFermion5D);
   

lib/qcd/action/fermion/WilsonFermion5D.h

@@ -214,6 +214,21 @@ namespace QCD {
     // Comms buffer
     std::vector<SiteHalfSpinor,alignedAllocator<SiteHalfSpinor> >  comm_buf;
     
+    ///////////////////////////////////////////////////////////////
+    // Conserved current utilities
+    ///////////////////////////////////////////////////////////////
+    void ContractConservedCurrent(PropagatorField &q_in_1,
+                                  PropagatorField &q_in_2,
+                                  PropagatorField &q_out,
+                                  Current curr_type, 
+                                  unsigned int mu);
+    void SeqConservedCurrent(PropagatorField &q_in,
+                             PropagatorField &q_out,
+                             Current curr_type,
+                             unsigned int mu,
+                             std::vector<Real> mom,
+                             unsigned int tmin,
+                             unsigned int tmax);
   };
 
 }}

lib/qcd/action/fermion/WilsonKernels.cc

@@ -281,6 +281,172 @@ void WilsonKernels<Impl>::DhopDir( StencilImpl &st, DoubledGaugeField &U,SiteHal
   vstream(out._odata[sF], result);
 }
 
+/*******************************************************************************
+ * Conserved current utilities for Wilson fermions, for contracting propagators
+ * to make a conserved current sink or inserting the conserved current 
+ * sequentially. Common to both 4D and 5D.
+ ******************************************************************************/
+// N.B. Functions below assume a -1/2 factor within U.
+#define WilsonCurrentFwd(expr, mu) ((expr - Gamma::gmu[mu]*expr))
+#define WilsonCurrentBwd(expr, mu) ((expr + Gamma::gmu[mu]*expr))
+
+/*******************************************************************************
+ * Name: ContractConservedCurrentSiteFwd
+ * Operation: (1/2) * q2[x] * U(x) * (g[mu] - 1) * q1[x + mu]
+ * Notes: - DoubledGaugeField U assumed to contain -1/2 factor.
+ *        - Pass in q_in_1 shifted in +ve mu direction.
+ ******************************************************************************/
+template<class Impl>
+void WilsonKernels<Impl>::ContractConservedCurrentSiteFwd(
+                                                  const SitePropagator &q_in_1,
+                                                  const SitePropagator &q_in_2,
+                                                  SitePropagator &q_out,
+                                                  DoubledGaugeField &U,
+                                                  unsigned int sU,
+                                                  unsigned int mu,
+                                                  bool switch_sign)
+{
+    SitePropagator result, tmp;
+    Gamma g5(Gamma::Algebra::Gamma5);
+    Impl::multLinkProp(tmp, U._odata[sU], q_in_1, mu);
+    result = g5 * adj(q_in_2) * g5 * WilsonCurrentFwd(tmp, mu);
+    if (switch_sign)
+    {
+        q_out -= result;
+    }
+    else
+    {
+        q_out += result;
+    }
+}
+
+/*******************************************************************************
+ * Name: ContractConservedCurrentSiteBwd
+ * Operation: (1/2) * q2[x + mu] * U^dag(x) * (g[mu] + 1) * q1[x]
+ * Notes: - DoubledGaugeField U assumed to contain -1/2 factor.
+ *        - Pass in q_in_2 shifted in +ve mu direction.
+ ******************************************************************************/
+template<class Impl>
+void WilsonKernels<Impl>::ContractConservedCurrentSiteBwd(
+                                                  const SitePropagator &q_in_1,
+                                                  const SitePropagator &q_in_2,
+                                                  SitePropagator &q_out,
+                                                  DoubledGaugeField &U,
+                                                  unsigned int sU,
+                                                  unsigned int mu,
+                                                  bool switch_sign)
+{
+    SitePropagator result, tmp;
+    Gamma g5(Gamma::Algebra::Gamma5);
+    Impl::multLinkProp(tmp, U._odata[sU], q_in_1, mu + Nd);
+    result = g5 * adj(q_in_2) * g5 * WilsonCurrentBwd(tmp, mu);
+    if (switch_sign)
+    {
+        q_out += result;
+    }
+    else
+    {
+        q_out -= result;
+    }
+}
+
+// G-parity requires more specialised implementation.
+#define NO_CURR_SITE(Impl) \
+template <> \
+void WilsonKernels<Impl>::ContractConservedCurrentSiteFwd( \
+                                                  const SitePropagator &q_in_1, \
+                                                  const SitePropagator &q_in_2, \
+                                                  SitePropagator &q_out,        \
+                                                  DoubledGaugeField &U,         \
+                                                  unsigned int sU,              \
+                                                  unsigned int mu,              \
+                                                  bool switch_sign)             \
+{ \
+    assert(0); \
+} \
+template <> \
+void WilsonKernels<Impl>::ContractConservedCurrentSiteBwd( \
+                                                  const SitePropagator &q_in_1, \
+                                                  const SitePropagator &q_in_2, \
+                                                  SitePropagator &q_out,        \
+                                                  DoubledGaugeField &U,         \
+                                                  unsigned int mu,              \
+                                                  unsigned int sU,              \
+                                                  bool switch_sign)             \
+{ \
+    assert(0); \
+}
+
+NO_CURR_SITE(GparityWilsonImplF);
+NO_CURR_SITE(GparityWilsonImplD);
+NO_CURR_SITE(GparityWilsonImplFH);
+NO_CURR_SITE(GparityWilsonImplDF);
+
+
+/*******************************************************************************
+ * Name: SeqConservedCurrentSiteFwd
+ * Operation: (1/2) * U(x) * (g[mu] - 1) * q[x + mu]
+ * Notes: - DoubledGaugeField U assumed to contain -1/2 factor.
+ *        - Pass in q_in shifted in +ve mu direction.
+ ******************************************************************************/
+template<class Impl>
+void WilsonKernels<Impl>::SeqConservedCurrentSiteFwd(const SitePropagator &q_in,
+                                                     SitePropagator &q_out,
+                                                     DoubledGaugeField &U,
+                                                     unsigned int sU,
+                                                     unsigned int mu,
+                                                     vInteger t_mask,
+                                                     bool switch_sign)
+{
+    SitePropagator result;
+    Impl::multLinkProp(result, U._odata[sU], q_in, mu);
+    result = WilsonCurrentFwd(result, mu);
+
+    // Zero any unwanted timeslice entries.
+    result = predicatedWhere(t_mask, result, 0.*result);
+
+    if (switch_sign)
+    {
+        q_out -= result;
+    }
+    else
+    {
+        q_out += result;
+    }
+}
+
+/*******************************************************************************
+ * Name: SeqConservedCurrentSiteFwd
+ * Operation: (1/2) * U^dag(x) * (g[mu] + 1) * q[x - mu]
+ * Notes: - DoubledGaugeField U assumed to contain -1/2 factor.
+ *        - Pass in q_in shifted in -ve mu direction.
+ ******************************************************************************/
+template<class Impl>
+void WilsonKernels<Impl>::SeqConservedCurrentSiteBwd(const SitePropagator &q_in, 
+                                                     SitePropagator &q_out,
+                                                     DoubledGaugeField &U,
+                                                     unsigned int sU,
+                                                     unsigned int mu,
+                                                     vInteger t_mask,
+                                                     bool switch_sign)
+{
+    SitePropagator result;
+    Impl::multLinkProp(result, U._odata[sU], q_in, mu + Nd);
+    result = WilsonCurrentBwd(result, mu);
+
+    // Zero any unwanted timeslice entries.
+    result = predicatedWhere(t_mask, result, 0.*result);
+
+    if (switch_sign)
+    {
+        q_out += result;
+    }
+    else
+    {
+        q_out -= result;
+    }
+}
+
 FermOpTemplateInstantiate(WilsonKernels);
 AdjointFermOpTemplateInstantiate(WilsonKernels);
 TwoIndexFermOpTemplateInstantiate(WilsonKernels);

lib/qcd/action/fermion/WilsonKernels.h

@@ -180,6 +180,38 @@ public:
   void DhopDir(StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor * buf,
 		       int sF, int sU, const FermionField &in, FermionField &out, int dirdisp, int gamma);
       
+  //////////////////////////////////////////////////////////////////////////////
+  // Utilities for inserting Wilson conserved current.
+  //////////////////////////////////////////////////////////////////////////////
+  void ContractConservedCurrentSiteFwd(const SitePropagator &q_in_1,
+                                       const SitePropagator &q_in_2,
+                                       SitePropagator &q_out,
+                                       DoubledGaugeField &U,
+                                       unsigned int sU,
+                                       unsigned int mu,
+                                       bool switch_sign = false);
+  void ContractConservedCurrentSiteBwd(const SitePropagator &q_in_1,
+                                       const SitePropagator &q_in_2,
+                                       SitePropagator &q_out,
+                                       DoubledGaugeField &U,
+                                       unsigned int sU,
+                                       unsigned int mu,
+                                       bool switch_sign = false);
+  void SeqConservedCurrentSiteFwd(const SitePropagator &q_in, 
+                                  SitePropagator &q_out,
+                                  DoubledGaugeField &U,
+                                  unsigned int sU,
+                                  unsigned int mu,
+                                  vInteger t_mask,
+                                  bool switch_sign = false);
+  void SeqConservedCurrentSiteBwd(const SitePropagator &q_in,
+                                  SitePropagator &q_out,
+                                  DoubledGaugeField &U,
+                                  unsigned int sU,
+                                  unsigned int mu,
+                                  vInteger t_mask,
+                                  bool switch_sign = false);
+
 private:
      // Specialised variants
   void GenericDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, SiteHalfSpinor * buf,

lib/util/Init.cc

@@ -204,7 +204,7 @@ std::string GridCmdVectorIntToString(const std::vector<int> & vec){
 // Reinit guard
 /////////////////////////////////////////////////////////
 static int Grid_is_initialised = 0;
-
+static MemoryStats dbgMemStats;
 
 void Grid_init(int *argc,char ***argv)
 {
@@ -251,6 +251,11 @@ void Grid_init(int *argc,char ***argv)
     assert(fp!=(FILE *)NULL);
   }
 
+  if( GridCmdOptionExists(*argv,*argv+*argc,"--debug-mem") ){
+    MemoryProfiler::debug = true;
+    MemoryProfiler::stats = &dbgMemStats;
+  }
+
   ////////////////////////////////////
   // Banner
   ////////////////////////////////////
@@ -324,6 +329,7 @@ void Grid_init(int *argc,char ***argv)
     std::cout<<GridLogMessage<<"  --decomposition : report on default omp,mpi and simd decomposition"<<std::endl;    
     std::cout<<GridLogMessage<<"  --debug-signals : catch sigsegv and print a blame report"<<std::endl;
     std::cout<<GridLogMessage<<"  --debug-stdout  : print stdout from EVERY node"<<std::endl;
+    std::cout<<GridLogMessage<<"  --debug-mem     : print Grid allocator activity"<<std::endl;
     std::cout<<GridLogMessage<<"  --notimestamp   : suppress millisecond resolution stamps"<<std::endl;
     std::cout<<GridLogMessage<<std::endl;
     std::cout<<GridLogMessage<<"Performance:"<<std::endl;

tests/hadrons/Test_hadrons.hpp

@@ -0,0 +1,659 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+
+ Source file: tests/hadrons/Test_hadrons.hpp
+
+ Copyright (C) 2017
+
+ Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/*******************************************************************************
+ * Macros to reduce code duplication.
+ ******************************************************************************/
+// Common initialisation
+#define HADRONS_DEFAULT_INIT \
+    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;
+
+#define HADRONS_DEFAULT_GLOBALS(application) \
+{ \
+    Application::GlobalPar globalPar;           \
+    globalPar.trajCounter.start    = 1500;      \
+    globalPar.trajCounter.end      = 1520;      \
+    globalPar.trajCounter.step     = 20;        \
+    globalPar.seed                 = "1 2 3 4"; \
+    globalPar.genetic.maxGen       = 1000;      \
+    globalPar.genetic.maxCstGen    = 200;       \
+    globalPar.genetic.popSize      = 20;        \
+    globalPar.genetic.mutationRate = .1;        \
+    application.setPar(globalPar);              \
+}
+
+// Useful definitions
+#define ZERO_MOM "0. 0. 0. 0."
+#define INIT_INDEX(s, n) (std::string(s) + "_" + std::to_string(n))
+#define ADD_INDEX(s, n) (s + "_" + std::to_string(n))
+#define LABEL_3PT(s, t1, t2) ADD_INDEX(INIT_INDEX(s, t1), t2)
+#define LABEL_4PT(s, t1, t2, t3) ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3)
+#define LABEL_4PT_NOISE(s, t1, t2, t3, nn) ADD_INDEX(ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3), nn)
+#define LABEL_5D(s) s + "_5d";
+
+// Wall source/sink macros
+#define NAME_3MOM_WALL_SOURCE(t, mom) ("wall_" + std::to_string(t) + "_" + mom)
+#define NAME_WALL_SOURCE(t) NAME_3MOM_WALL_SOURCE(t, ZERO_MOM)
+#define NAME_POINT_SOURCE(pos) ("point_" + pos)
+
+// Meson module "gammas" special values
+#define ALL_GAMMAS "all"
+
+#define MAKE_3MOM_WALL_PROP(tW, mom, propName, solver)\
+{\
+    std::string srcName = NAME_3MOM_WALL_SOURCE(tW, mom);\
+    makeWallSource(application, srcName, tW, mom);\
+    makePropagator(application, propName, srcName, solver);\
+}
+
+#define MAKE_WALL_PROP(tW, propName, solver)\
+        MAKE_3MOM_WALL_PROP(tW, ZERO_MOM, propName, solver)
+
+// Sequential source macros
+#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, seqPropName, solver, gamma)\
+{\
+    std::string srcName = seqPropName + "_src";\
+    makeSequentialSource(application, srcName, qSrc, tS, gamma, mom);\
+    makePropagator(application, seqPropName, srcName, solver);\
+}
+
+// Point source macros
+#define MAKE_POINT_PROP(pos, propName, solver)\
+{\
+    std::string srcName = NAME_POINT_SOURCE(pos);\
+    makePointSource(application, srcName, pos);\
+    makePropagator(application, propName, srcName, solver);\
+}
+
+/*******************************************************************************
+ * Action setups.
+ ******************************************************************************/
+
+/*******************************************************************************
+ * Name: makeWilsonAction
+ * Parameters: application - main application that stores modules.
+ *             actionName  - name of action module to create.
+ *             gaugeField  - gauge field module.     
+ *             mass        - quark mass.
+ *             boundary    - fermion boundary conditions (default to periodic
+ *                           space, antiperiodic time).
+ * Returns: None.
+ ******************************************************************************/
+inline void makeWilsonAction(Application &application, std::string actionName,
+                             std::string &gaugeField, double mass,
+                             std::string boundary = "1 1 1 -1")
+{
+    if (!(VirtualMachine::getInstance().hasModule(actionName)))
+    {
+        MAction::Wilson::Par actionPar;
+        actionPar.gauge = gaugeField;
+        actionPar.mass  = mass;
+        actionPar.boundary = boundary;
+        application.createModule<MAction::Wilson>(actionName, actionPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeDWFAction
+ * Parameters: application - main application that stores modules.
+ *             actionName  - name of action module to create.
+ *             gaugeField  - gauge field module.     
+ *             mass        - quark mass.
+ *             M5          - domain wall height.
+ *             Ls          - fifth dimension extent.
+ *             boundary    - fermion boundary conditions (default to periodic
+ *                           space, antiperiodic time).
+ * Returns: None.
+ ******************************************************************************/
+inline void makeDWFAction(Application &application, std::string actionName,
+                          std::string &gaugeField, double mass, double M5,
+                          unsigned int Ls, std::string boundary = "1 1 1 -1")
+{
+    if (!(VirtualMachine::getInstance().hasModule(actionName)))
+    {
+        MAction::DWF::Par actionPar;
+        actionPar.gauge = gaugeField;
+        actionPar.Ls    = Ls;
+        actionPar.M5    = M5;
+        actionPar.mass  = mass;
+        actionPar.boundary = boundary;
+        application.createModule<MAction::DWF>(actionName, actionPar);
+    }
+}
+
+/*******************************************************************************
+ * Functions for propagator construction.
+ ******************************************************************************/
+ 
+/*******************************************************************************
+ * Name: makeRBPrecCGSolver
+ * Purpose: Make RBPrecCG solver module for specified action.
+ * Parameters: application - main application that stores modules.
+ *             solverName  - name of solver module to create.
+ *             actionName  - action module corresponding to propagators to be
+ *                           computed.
+ *             residual    - CG target residual.
+ * Returns: None.
+ ******************************************************************************/
+inline void makeRBPrecCGSolver(Application &application, std::string &solverName,
+                               std::string &actionName, double residual = 1e-8)
+{
+    if (!(VirtualMachine::getInstance().hasModule(solverName)))
+    {
+        MSolver::RBPrecCG::Par solverPar;
+        solverPar.action   = actionName;
+        solverPar.residual = residual;
+        application.createModule<MSolver::RBPrecCG>(solverName,
+                                                    solverPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makePointSource
+ * Purpose: Construct point source and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             pos         - Position of point source.
+ * Returns: None.
+ ******************************************************************************/
+inline void makePointSource(Application &application, std::string srcName,
+                            std::string pos)
+{
+    // If the source already exists, don't make the module again.
+    if (!(VirtualMachine::getInstance().hasModule(srcName)))
+    {
+        MSource::Point::Par pointPar;
+        pointPar.position = pos;
+        application.createModule<MSource::Point>(srcName, pointPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeSequentialSource
+ * Purpose: Construct sequential source and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             qSrc        - Input quark for sequential inversion.
+ *             tS          - sequential source timeslice.
+ *             mom         - momentum insertion (default is zero).
+ * Returns: None.
+ ******************************************************************************/
+inline void makeSequentialSource(Application &application, std::string srcName,
+                                 std::string qSrc, unsigned int tS,
+                                 Gamma::Algebra gamma = Gamma::Algebra::GammaT,
+                                 std::string mom = ZERO_MOM)
+{
+    // If the source already exists, don't make the module again.
+    if (!(VirtualMachine::getInstance().hasModule(srcName)))
+    {
+        MSource::SeqGamma::Par seqPar;
+        seqPar.q   = qSrc;
+        seqPar.tA  = tS;
+        seqPar.tB  = tS;
+        seqPar.mom = mom;
+        seqPar.gamma = gamma;
+        application.createModule<MSource::SeqGamma>(srcName, seqPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeConservedSequentialSource
+ * Purpose: Construct sequential source with conserved current insertion and 
+ *          add to application module.
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             qSrc        - Input quark for sequential inversion.
+ *             actionName  - action corresponding to quark.
+ *             tS          - sequential source timeslice.
+ *             curr        - conserved current type to insert.
+ *             mu          - Lorentz index of current to insert.
+ *             mom         - momentum insertion (default is zero).
+ * Returns: None.
+ ******************************************************************************/
+inline void makeConservedSequentialSource(Application &application,
+                                          std::string &srcName,
+                                          std::string &qSrc, 
+                                          std::string &actionName,
+                                          unsigned int tS,
+                                          Current curr,
+                                          unsigned int mu,
+                                          std::string mom = ZERO_MOM)
+{
+    // If the source already exists, don't make the module again.
+    if (!(VirtualMachine::getInstance().hasModule(srcName)))
+    {
+        MSource::SeqConserved::Par seqPar;
+        seqPar.q         = qSrc;
+        seqPar.action    = actionName;
+        seqPar.tA        = tS;
+        seqPar.tB        = tS;
+        seqPar.curr_type = curr;
+        seqPar.mu        = mu;
+        seqPar.mom       = mom;
+        application.createModule<MSource::SeqConserved>(srcName, seqPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeNoiseSource
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             tA          - lower source timeslice limit.
+ *             tB          - upper source timeslice limit.
+ * Returns: None.
+ ******************************************************************************/
+inline void makeNoiseSource(Application &application, std::string &srcName,
+                            unsigned int tA, unsigned int tB)
+{
+    if (!(VirtualMachine::getInstance().hasModule(srcName)))
+    {
+        MSource::Z2::Par noisePar;
+        noisePar.tA = tA;
+        noisePar.tB = tB;
+        application.createModule<MSource::Z2>(srcName, noisePar);
+    }
+ }
+
+/*******************************************************************************
+ * Name: makeWallSource
+ * Purpose: Construct wall source and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             tW          - wall source timeslice.
+ *             mom         - momentum insertion (default is zero).
+ * Returns: None.
+ ******************************************************************************/
+inline void makeWallSource(Application &application, std::string &srcName,
+                           unsigned int tW, std::string mom = ZERO_MOM)
+{
+    // If the source already exists, don't make the module again.
+    if (!(VirtualMachine::getInstance().hasModule(srcName)))
+    {
+        MSource::Wall::Par wallPar;
+        wallPar.tW  = tW;
+        wallPar.mom = mom;
+        application.createModule<MSource::Wall>(srcName, wallPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makePointSink
+ * Purpose: Create function for point sink smearing of a propagator.
+ * Parameters: application - main application that stores modules.
+ *             propName    - name of input propagator.
+ *             sinkFnct    - name of output sink smearing module.
+ *             mom         - momentum insertion (default is zero).
+ * Returns: None.
+ ******************************************************************************/
+inline void makePointSink(Application &application, std::string &sinkFnct,
+                          std::string mom = ZERO_MOM)
+{
+    // If the sink function already exists, don't make it again.
+    if (!(VirtualMachine::getInstance().hasModule(sinkFnct)))
+    {
+        MSink::Point::Par pointPar;
+        pointPar.mom = mom;
+        application.createModule<MSink::Point>(sinkFnct, pointPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: sinkSmear
+ * Purpose: Perform sink smearing of a propagator.
+ * Parameters: application - main application that stores modules.
+ *             sinkFnct    - sink smearing module.
+ *             propName    - propagator to smear.
+ *             smearedProp - name of output smeared propagator.
+ * Returns: None.
+ ******************************************************************************/
+inline void sinkSmear(Application &application, std::string &sinkFnct,
+                      std::string &propName, std::string &smearedProp)
+{
+    // If the propagator has already been smeared, don't smear it again.
+    if (!(VirtualMachine::getInstance().hasModule(smearedProp)))
+    {
+        MSink::Smear::Par smearPar;
+        smearPar.q    = propName;
+        smearPar.sink = sinkFnct;
+        application.createModule<MSink::Smear>(smearedProp, smearPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makePropagator
+ * Purpose: Construct source and propagator then add to application module.
+ * Parameters: application - main application that stores modules.
+ *             propName    - name of propagator module to create.
+ *             srcName     - name of source module to use.
+ *             solver      - solver to use (default is CG).
+ * Returns: None.
+ ******************************************************************************/
+inline void makePropagator(Application &application, std::string &propName,
+                           std::string &srcName, std::string &solver)
+{
+    // If the propagator already exists, don't make the module again.
+    if (!(VirtualMachine::getInstance().hasModule(propName)))
+    {
+        MFermion::GaugeProp::Par quarkPar;
+        quarkPar.source = srcName;
+        quarkPar.solver = solver;
+        application.createModule<MFermion::GaugeProp>(propName, quarkPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeLoop
+ * Purpose: Use noise source and inversion result to make loop propagator, then 
+ *          add to application module.
+ * Parameters: application - main application that stores modules.
+ *             propName    - name of propagator module to create.
+ *             srcName     - name of noise source module to use.
+ *             resName     - name of inversion result on given noise source.
+ * Returns: None.
+ ******************************************************************************/
+inline void makeLoop(Application &application, std::string &propName,
+                     std::string &srcName, std::string &resName)
+{
+    // If the loop propagator already exists, don't make the module again.
+    if (!(VirtualMachine::getInstance().hasModule(propName)))
+    {
+        MLoop::NoiseLoop::Par loopPar;
+        loopPar.q   = resName;
+        loopPar.eta = srcName;
+        application.createModule<MLoop::NoiseLoop>(propName, loopPar);
+    }
+}
+
+/*******************************************************************************
+ * Contraction module creation.
+ ******************************************************************************/
+
+/*******************************************************************************
+ * Name: mesonContraction
+ * Purpose: Create meson contraction module and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             modName     - unique module name.
+ *             output      - name of output files.
+ *             q1          - quark propagator 1.
+ *             q2          - quark propagator 2.
+ *             sink        - sink smearing module.
+ *             gammas      - gamma insertions at source and sink.
+ * Returns: None.
+ ******************************************************************************/
+inline void mesonContraction(Application &application, 
+                             std::string &modName, std::string &output, 
+                             std::string &q1, std::string &q2,
+                             std::string &sink,
+                             std::string gammas = "<Gamma5 Gamma5>")
+{
+    if (!(VirtualMachine::getInstance().hasModule(modName)))
+    {
+        MContraction::Meson::Par mesPar;
+        mesPar.output = output;
+        mesPar.q1 = q1;
+        mesPar.q2 = q2;
+        mesPar.sink = sink;
+        mesPar.gammas = gammas;
+        application.createModule<MContraction::Meson>(modName, mesPar);
+    }
+ }
+
+/*******************************************************************************
+ * Name: gamma3ptContraction
+ * Purpose: Create gamma3pt contraction module and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             npt         - specify n-point correlator (for labelling).
+ *             q1          - quark propagator 1, sink smeared.
+ *             q2          - quark propagator 2.
+ *             q3          - quark propagator 3.
+ *             label       - unique label to construct module name.
+ *             tSnk        - sink position of sink for q1.
+ *             gamma       - gamma insertions between q2 and q3.
+ * Returns: None.
+ ******************************************************************************/
+inline void gamma3ptContraction(Application &application, unsigned int npt, 
+                                std::string &q1, std::string &q2,
+                                std::string &q3, std::string &label,
+                                unsigned int tSnk = 0,
+                                Gamma::Algebra gamma = Gamma::Algebra::Identity)
+{
+    std::string modName = std::to_string(npt) + "pt_" + label;
+    if (!(VirtualMachine::getInstance().hasModule(modName)))
+    {
+        MContraction::Gamma3pt::Par gamma3ptPar;
+        gamma3ptPar.output = std::to_string(npt) + "pt/" + label;
+        gamma3ptPar.q1 = q1;
+        gamma3ptPar.q2 = q2;
+        gamma3ptPar.q3 = q3;
+        gamma3ptPar.tSnk = tSnk;
+        gamma3ptPar.gamma = gamma;
+        application.createModule<MContraction::Gamma3pt>(modName, gamma3ptPar);
+    }
+ }
+
+/*******************************************************************************
+ * Name: weakContraction[Eye,NonEye]
+ * Purpose: Create Weak Hamiltonian contraction module for Eye/NonEye topology
+ *          and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             npt         - specify n-point correlator (for labelling).
+ *             q1          - quark propagator 1.
+ *             q2          - quark propagator 2.
+ *             q3          - quark propagator 3.
+ *             q4          - quark propagator 4.
+ *             label       - unique label to construct module name.
+ *             tSnk        - time position of sink (for sink smearing).
+ * Returns: None.
+ ******************************************************************************/
+#define HW_CONTRACTION(top) \
+inline void weakContraction##top(Application &application, unsigned int npt,\
+                                 std::string &q1, std::string &q2, \
+                                 std::string &q3, std::string &q4, \
+                                 std::string &label, unsigned int tSnk = 0)\
+{\
+    std::string modName = std::to_string(npt) + "pt_" + label;\
+    if (!(VirtualMachine::getInstance().hasModule(modName)))\
+    {\
+        MContraction::WeakHamiltonian##top::Par weakPar;\
+        weakPar.output = std::to_string(npt) + "pt/" + label;\
+        weakPar.q1 = q1;\
+        weakPar.q2 = q2;\
+        weakPar.q3 = q3;\
+        weakPar.q4 = q4;\
+        weakPar.tSnk = tSnk;\
+        application.createModule<MContraction::WeakHamiltonian##top>(modName, weakPar);\
+    }\
+}
+HW_CONTRACTION(Eye)    // weakContractionEye
+HW_CONTRACTION(NonEye) // weakContractionNonEye
+
+/*******************************************************************************
+ * Name: disc0Contraction
+ * Purpose: Create contraction module for 4pt Weak Hamiltonian + current
+ *          disconnected topology for neutral mesons and add to application 
+ *          module.
+ * Parameters: application - main application that stores modules.
+ *             q1          - quark propagator 1.
+ *             q2          - quark propagator 2.
+ *             q3          - quark propagator 3.
+ *             q4          - quark propagator 4.
+ *             label       - unique label to construct module name.
+ * Returns: None.
+ ******************************************************************************/
+inline void disc0Contraction(Application &application, 
+                             std::string &q1, std::string &q2,
+                             std::string &q3, std::string &q4,
+                             std::string &label)
+{
+    std::string modName = "4pt_" + label;
+    if (!(VirtualMachine::getInstance().hasModule(modName)))
+    {
+        MContraction::WeakNeutral4ptDisc::Par disc0Par;
+        disc0Par.output = "4pt/" + label;
+        disc0Par.q1 = q1;
+        disc0Par.q2 = q2;
+        disc0Par.q3 = q3;
+        disc0Par.q4 = q4;
+        application.createModule<MContraction::WeakNeutral4ptDisc>(modName, disc0Par);
+    }
+ }
+
+/*******************************************************************************
+ * Name: discLoopContraction
+ * Purpose: Create contraction module for disconnected loop and add to
+ *          application module.
+ * Parameters: application - main application that stores modules.
+ *             q_loop      - loop quark propagator.
+ *             modName     - unique module name.
+ *             gamma       - gamma matrix to use in contraction.
+ * Returns: None.
+ ******************************************************************************/
+inline void discLoopContraction(Application &application,
+                                std::string &q_loop, std::string &modName,
+                                Gamma::Algebra gamma = Gamma::Algebra::Identity)
+{
+    if (!(VirtualMachine::getInstance().hasModule(modName)))
+    {
+        MContraction::DiscLoop::Par discPar;
+        discPar.output = "disc/" + modName;
+        discPar.q_loop = q_loop;
+        discPar.gamma  = gamma;
+        application.createModule<MContraction::DiscLoop>(modName, discPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeWITest
+ * Purpose: Create module to test Ward Identities for conserved current
+ *          contractions and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             modName     - name of module to create.
+ *             propName    - 4D quark propagator.
+ *             actionName  - action used to compute quark propagator.
+ *             mass        - mass of quark.
+ *             Ls          - length of 5th dimension (default = 1).
+ *             test_axial  - whether or not to check PCAC relation.
+ * Returns: None.
+ ******************************************************************************/
+inline void makeWITest(Application &application, std::string &modName,
+                       std::string &propName, std::string &actionName, 
+                       double mass, unsigned int Ls = 1, bool test_axial = false)
+{
+    if (!(VirtualMachine::getInstance().hasModule(modName)))
+    {
+        MContraction::WardIdentity::Par wiPar;
+        if (Ls > 1)
+        {
+            wiPar.q = LABEL_5D(propName);
+        }
+        else
+        {
+            wiPar.q = propName;
+        }
+        wiPar.action     = actionName;
+        wiPar.mass       = mass;
+        wiPar.test_axial = test_axial;
+        application.createModule<MContraction::WardIdentity>(modName, wiPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeSeqCurrComparison
+ * Purpose: Create module to compare sequential insertion of conserved current
+ *          against sink contraction and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             modName     - name of module to create.
+ *             propName    - quark propagator (point source), 5D if available.
+ *             seqName     - 4D quark propagator with sequential insertion of
+ *                           conserved current.
+ *             actionName  - action used to compute quark propagators.
+ *             origin      - origin of point source propagator.
+ *             t_J         - time at which sequential current is inserted.
+ *             mu          - Lorentz index of sequential current.
+ *             curr        - type of conserved current inserted.
+ * Returns: None.
+ ******************************************************************************/
+inline void makeSeqCurrComparison(Application &application, std::string &modName,
+                                 std::string &propName, std::string &seqName,
+                                 std::string &actionName, std::string &origin,
+                                 unsigned int t_J, unsigned int mu, Current curr)
+{
+    if (!(VirtualMachine::getInstance().hasModule(modName)))
+    {
+        MUtilities::TestSeqConserved::Par seqPar;
+        seqPar.q      = propName;
+        seqPar.qSeq   = seqName;
+        seqPar.action = actionName;
+        seqPar.origin = origin;
+        seqPar.t_J    = t_J;
+        seqPar.mu     = mu;
+        seqPar.curr   = curr;
+        application.createModule<MUtilities::TestSeqConserved>(modName, seqPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeSeqGamComparison
+ * Purpose: Create module to compare sequential insertion of gamma matrix
+ *          against sink contraction and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             modName     - name of module to create.
+ *             propName    - 4D quark propagator.
+ *             seqProp     - 4D quark propagator with sequential insertion of
+ *                           gamma matrix.
+ *             gamma       - Inserted gamma matrix.
+ *             t_g         - time at which gamma matrix is inserted 
+ *                           sequentially.
+ * Returns: None.
+ ******************************************************************************/
+inline void makeSeqGamComparison(Application &application, std::string &modName,
+                                 std::string &propName, std::string &seqProp,
+                                 std::string &origin, Gamma::Algebra gamma, 
+                                 unsigned int t_g)
+{
+    if (!(VirtualMachine::getInstance().hasModule(modName)))
+    {
+        MUtilities::TestSeqGamma::Par seqPar;
+        seqPar.q      = propName;
+        seqPar.qSeq   = seqProp;
+        seqPar.origin = origin;
+        seqPar.t_g    = t_g;
+        seqPar.gamma  = gamma;
+        application.createModule<MUtilities::TestSeqGamma>(modName, seqPar);
+    }
+}

tests/hadrons/Test_hadrons_3pt_contractions.cc

@@ -0,0 +1,122 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_3pt_contractions.cc
+ 
+ Copyright (C) 2017
+ 
+ Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+ 
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ 
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include "Test_hadrons.hpp"
+
+using namespace Grid;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    HADRONS_DEFAULT_INIT;
+
+    // run setup ///////////////////////////////////////////////////////////////
+    Application  application;
+    double       mass = 0.04;
+    double       M5   = 1.8;
+    unsigned int Ls   = 12;
+    unsigned int nt   = GridDefaultLatt()[Tp];
+    unsigned int t_i  = 0;
+    unsigned int t_f  = nt / 2;
+    std::string  mom  = "1. 0. 0. 0.";
+
+    // global parameters
+    HADRONS_DEFAULT_GLOBALS(application);
+
+    // gauge field
+    std::string gaugeField = "gauge";
+    application.createModule<MGauge::Unit>(gaugeField);
+
+    // Action & solver setup.
+    std::string action = "DWF";
+    std::string solver = "CG";
+    makeDWFAction(application, action, gaugeField, mass, M5, Ls);
+    makeRBPrecCGSolver(application, solver, action);
+
+    /***************************************************************************
+     * Weak Contraction test: Non-Eye class.
+     **************************************************************************/
+    // Make wall source propagators for each leg of 4-quark vertex.
+    std::string q_i_0 = "q_i_0";
+    std::string q_i_p = "q_i_p";
+    std::string q_f_0 = "q_f_0";
+    std::string q_f_p = "q_f_p";
+    MAKE_WALL_PROP(t_i, q_i_0, solver);
+    MAKE_WALL_PROP(t_f, q_f_0, solver);
+    MAKE_3MOM_WALL_PROP(t_i, mom, q_i_p, solver);
+    MAKE_3MOM_WALL_PROP(t_f, mom, q_f_p, solver);
+
+    // Perform contractions, zero and non-zero momentum.
+    std::string HW_CW_0 = LABEL_3PT("HW_CW_0", t_i, t_f);
+    std::string HW_CW_p = LABEL_3PT("HW_CW_p", t_i, t_f);
+    weakContractionNonEye(application, 3, q_i_0, q_i_0, q_f_0, q_f_0, HW_CW_0);
+    weakContractionNonEye(application, 3, q_i_0, q_i_p, q_f_p, q_f_0, HW_CW_p);
+
+    /***************************************************************************
+     * Weak Contraction test: Eye-class.
+     **************************************************************************/
+    // Create random propagator for loop.
+    std::string eta = "noise_source";
+    makeNoiseSource(application, eta, 0, nt - 1);
+    std::string loopProp = "loop";
+    std::string loopRes  = loopProp + "_res";
+    makePropagator(application, loopRes, eta, solver);
+    makeLoop(application, loopProp, eta, loopRes);
+
+    // Wall sink smear the propagator directly connecting the source & sink.
+    // (i.e. make point sink but smear before the contraction)
+    std::string wallSink = "wall_sink";
+    std::string qWall    = "q_wall";
+    makePointSink(application, wallSink);
+    sinkSmear(application, wallSink, q_i_0, qWall);
+
+    // Perform contractions, zero and non-zero momentum.
+    std::string HW_SE_0 = LABEL_3PT("HW_SE_0", t_i, t_f);
+    std::string HW_SE_p = LABEL_3PT("HW_SE_p", t_i, t_f);
+    weakContractionEye(application, 3, qWall, q_i_0, q_f_p, loopProp, HW_SE_0, t_f);
+    weakContractionEye(application, 3, qWall, q_i_p, q_f_p, loopProp, HW_SE_p, t_f);
+
+    /***************************************************************************
+     * Gamma insertion test.
+     **************************************************************************/
+    Gamma::Algebra gamma = Gamma::Algebra::GammaT;
+    std::string sd_0 = LABEL_3PT("sd_0", t_i, t_f);
+    std::string sd_p = LABEL_3PT("sd_p", t_i, t_f);
+    gamma3ptContraction(application, 3, qWall, q_i_0, q_f_0, sd_0, t_f, gamma);
+    gamma3ptContraction(application, 3, qWall, q_i_p, q_f_p, sd_p, t_f, gamma);
+
+    // execution
+    application.saveParameterFile("ContractionTest3pt.xml");
+    application.run();
+
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons_conserved_current.cc

@@ -0,0 +1,150 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_conserved_current.cc
+ 
+ Copyright (C) 2017
+ 
+ Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+ 
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+ 
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+ 
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ 
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include "Test_hadrons.hpp"
+
+using namespace Grid;
+using namespace Hadrons;
+
+inline void setupSeqCurrTests(Application &application, std::string modStem, 
+                              std::string &pointProp, std::string &seqStem,
+                              std::string &actionName, std::string &solverName,
+                              std::string &origin, Current curr, 
+                              unsigned int t_J, unsigned int mu,
+                              unsigned int Ls = 1)
+{
+    std::string modName = ADD_INDEX(modStem, mu);
+    std::string seqProp = ADD_INDEX(seqStem, mu);
+    std::string seqSrc  = seqProp + "_src";
+
+    // 5D actions require 5D propagator as input for conserved current
+    // insertions.
+    std::string propIn;
+    if (Ls > 1) 
+    {
+        propIn = LABEL_5D(pointProp); 
+    }
+    else
+    { 
+        propIn = pointProp; 
+    }
+
+    makeConservedSequentialSource(application, seqSrc, propIn,
+                                  actionName, t_J, curr, mu);
+    makePropagator(application, seqProp, seqSrc, solverName);
+    makeSeqCurrComparison(application, modName, propIn, seqProp, 
+                          actionName, origin, t_J, mu, curr);
+}
+
+inline void setupWardIdentityTests(Application &application,
+                                   std::string &actionName, 
+                                   double mass,
+                                   unsigned int Ls = 1, 
+                                   bool perform_axial_tests = false)
+{
+    // solver
+    std::string solverName = actionName + "_CG";
+    makeRBPrecCGSolver(application, solverName, actionName);
+
+    unsigned int nt = GridDefaultLatt()[Tp];
+    unsigned int t_J = nt/2;
+
+    /***************************************************************************
+     * Conserved current sink contractions: use a single point propagator for
+     * the Ward Identity test.
+     **************************************************************************/
+    std::string pointProp = actionName + "_q_0";
+    std::string origin    = "0 0 0 0";
+    std::string modName   = actionName + " Ward Identity Test";
+    MAKE_POINT_PROP(origin, pointProp, solverName);
+    makeWITest(application, modName, pointProp, actionName, mass, Ls,
+               perform_axial_tests);
+
+    /***************************************************************************
+     * Conserved current tests with sequential insertion of vector/axial 
+     * current. If above Ward Identity passes, sufficient to test sequential
+     * insertion of conserved current agrees with contracted version.
+     **************************************************************************/
+    // Compare sequential insertion to contraction. Should be enough to perform 
+    // for time and one space component.
+    std::string seqStem = ADD_INDEX(pointProp + "seq_V", t_J);
+    std::string modStem = actionName + " Vector Sequential Test mu";
+    setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName, 
+                      solverName, origin, Current::Vector, t_J, Tp, Ls);
+    setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName, 
+                      solverName, origin, Current::Vector, t_J, Xp, Ls);
+
+    // Perform axial tests only if partially-conserved axial current exists for
+    // the action.
+    if (perform_axial_tests)
+    {
+        seqStem = ADD_INDEX(pointProp + "seq_A", t_J);
+        modStem = actionName + " Axial Sequential Test mu";
+        setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName, 
+                          solverName, origin, Current::Axial, t_J, Tp, Ls);
+        setupSeqCurrTests(application, modStem, pointProp, seqStem, actionName, 
+                          solverName, origin, Current::Axial, t_J, Xp, Ls);
+    }
+}
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    HADRONS_DEFAULT_INIT;
+
+    // run setup ///////////////////////////////////////////////////////////////
+    Application  application;
+    double       mass = 0.04;
+    double       M5 = 1.8;
+    unsigned int Ls = 12;
+
+    // global parameters
+    HADRONS_DEFAULT_GLOBALS(application);
+
+    // gauge field
+    std::string gaugeField = "gauge";
+    application.createModule<MGauge::Unit>(gaugeField);
+
+    // Setup each action and the conserved current tests relevant to it.
+    std::string actionName = "DWF";
+    makeDWFAction(application, actionName, gaugeField, mass, M5, Ls);
+    setupWardIdentityTests(application, actionName, mass, Ls, true);
+
+    actionName = "Wilson";
+    makeWilsonAction(application, actionName, gaugeField, mass);
+    setupWardIdentityTests(application, actionName, mass);
+
+    // execution
+    application.saveParameterFile("ConservedCurrentTest.xml");
+    application.run();
+
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons_meson_conserved_3pt.cc

@@ -0,0 +1,115 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+
+ Source file: tests/hadrons/Test_hadrons_meson_conserved_3pt.cc
+
+ Copyright (C) 2017
+
+ Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include "Test_hadrons.hpp"
+
+using namespace Grid;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    HADRONS_DEFAULT_INIT;
+    
+    // run setup ///////////////////////////////////////////////////////////////
+    Application              application;
+    
+    // actions parameters
+    double        mass = 0.04;
+    unsigned int  Ls = 16;
+    double        M5 = 1.8;
+    
+    // kinematics
+    unsigned int  nt    = GridDefaultLatt()[Tp];
+    unsigned int  tSrc  = 0;
+    unsigned int  tJ    = nt / 4;
+    std::string   kmom  = "0. 0. 0. 0.";
+    std::string   pmom  = "1. 0. 0. 0.";
+
+    // Global parameters.
+    HADRONS_DEFAULT_GLOBALS(application);
+
+    // Unit gauge field.
+    std::string gaugeField = "Unit gauge";
+    application.createModule<MGauge::Unit>(gaugeField);
+
+    // DWF action
+    std::string actionName = "DWF";
+    makeDWFAction(application, actionName, gaugeField, mass, M5, Ls);
+
+    // Solver
+    std::string solver = "CG";
+    makeRBPrecCGSolver(application, solver, actionName);
+    
+    // main test body //////////////////////////////////////////////////////////
+    // Point sink modules.
+    std::string sink_0 = "sink_0";
+    std::string sink_p = "sink_p";
+    MSink::Point::Par sinkPar;
+    sinkPar.mom = kmom;
+    application.createModule<MSink::ScalarPoint>(sink_0, sinkPar);
+    sinkPar.mom = pmom;
+    application.createModule<MSink::ScalarPoint>(sink_p, sinkPar);
+
+    // 2pt pion contraction, zero momentum.
+    std::string q_0 = "Q_0";
+    MAKE_WALL_PROP(tSrc, q_0, solver);
+    std::string modName = INIT_INDEX("2pt_pion_WP", tSrc);
+    std::string output  = "2pt/pion_WP_0";
+    mesonContraction(application, modName, output, q_0, q_0, sink_0);
+
+    // 2pt pion contraction, with momentum p.
+    std::string q_p = "Q_p";
+    MAKE_3MOM_WALL_PROP(tSrc, pmom, q_p, solver);
+    modName = INIT_INDEX("2pt_pion_WP_p", tSrc);
+    output  = "2pt/pion_WP_p";
+    mesonContraction(application, modName, output, q_0, q_p, sink_p);
+
+    // 3pt pion(0) -> pion(p), with sequentially inserted vector current in 
+    // time direction.
+    std::string qSeq    = q_0 + INIT_INDEX("_seq_Vc3", tJ);
+    std::string q5d     = LABEL_5D(q_0); // Need 5D prop for DWF conserved current.
+    std::string srcName = qSeq + "_src";
+    modName = LABEL_3PT("3pt_pion_Vc3", tSrc, tJ);
+    output  = "3pt/pion_Vc3_p";
+    makeConservedSequentialSource(application, srcName, q5d, actionName,
+                                  tJ, Current::Vector, Tp, pmom);
+    makePropagator(application, qSeq, srcName, solver);
+    mesonContraction(application, modName, output, q_0, qSeq, sink_p);
+
+    std::string par_file_name = "conserved_3pt.xml";
+    application.saveParameterFile(par_file_name);
+    application.run();
+
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+
+    return EXIT_SUCCESS;
+}
+    
+    

tests/hadrons/Test_hadrons_quark.cc

@@ -0,0 +1,156 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+
+ Source file: tests/hadrons/Test_hadrons_quark.cc
+
+ Copyright (C) 2017
+
+ Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include "Test_hadrons.hpp"
+#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+/*******************************************************************************
+ * Unit test functions within Quark module.
+ ******************************************************************************/
+
+// Alternative 4D & 5D projections
+template<class vobj>
+inline void make_4D_with_gammas(Lattice<vobj> &in_5d, Lattice<vobj> &out_4d, int Ls)
+{
+    GridBase *_grid(out_4d._grid);
+    Lattice<vobj> tmp(_grid);
+    Gamma G5(Gamma::Algebra::Gamma5);
+
+    ExtractSlice(tmp, in_5d, 0, 0);
+    out_4d = 0.5 * (tmp - G5*tmp);
+    ExtractSlice(tmp, in_5d, Ls - 1, 0);
+    out_4d += 0.5 * (tmp + G5*tmp);
+}
+
+template<class vobj>
+inline void make_5D_with_gammas(Lattice<vobj> &in_4d, Lattice<vobj> &out_5d, int Ls)
+{
+    out_5d = zero;
+    Gamma G5(Gamma::Algebra::Gamma5);
+    GridBase *_grid(in_4d._grid);
+    Lattice<vobj> tmp(_grid);
+
+    tmp = 0.5 * (in_4d + G5*in_4d);
+    InsertSlice(tmp, out_5d, 0, 0);
+    tmp = 0.5 * (in_4d - G5*in_4d);
+    InsertSlice(tmp, out_5d, Ls - 1, 0);
+}
+
+int main(int argc, char **argv)
+{
+    /***************************************************************************
+     * Initialisation.
+     **************************************************************************/
+    Grid_init(&argc, &argv);
+
+    std::vector<int> latt_size   = GridDefaultLatt();
+    std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
+    std::vector<int> mpi_layout  = GridDefaultMpi();
+
+    const int Ls = 8;
+
+    GridCartesian   UGrid(latt_size,simd_layout,mpi_layout);
+    GridCartesian   *FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, &UGrid);
+    GridSerialRNG   sRNG;
+    GridParallelRNG pRNG(&UGrid);
+
+    std::vector<int> seeds4({1,2,3,4});
+    std::vector<int> seeds5({5,6,7,8});
+    GridParallelRNG  rng4(&UGrid);
+    GridParallelRNG  rng5(FGrid);
+    rng4.SeedFixedIntegers(seeds4);
+    rng5.SeedFixedIntegers(seeds5);
+
+    /***************************************************************************
+     * Build a 4D random source, and convert it to 5D.
+     **************************************************************************/
+    LatticeFermion test4(&UGrid);
+    LatticeFermion test5(FGrid);
+    LatticeFermion check5(FGrid);
+
+    gaussian(rng4, test4);
+    make_5D(test4, test5, Ls);
+    make_5D_with_gammas(test4, check5, Ls);
+    test5 -= check5;
+    std::cout << "4D -> 5D comparison, diff = " << Grid::sqrt(norm2(test5)) << std::endl;
+
+    /***************************************************************************
+     * Build a 5D random source, and project down to 4D.
+     **************************************************************************/
+    LatticeFermion check4(&UGrid);
+    gaussian(rng5, test5);
+    check5 = test5;
+
+    make_4D(test5, test4, Ls);
+    make_4D_with_gammas(check5, check4, Ls);
+    test4 -= check4;
+    std::cout << "5D -> 4D comparison, diff = " << Grid::sqrt(norm2(test4)) << std::endl;
+
+    /***************************************************************************
+     * Convert a propagator to a fermion & back.
+     **************************************************************************/
+    LatticeFermion ferm(&UGrid);
+    LatticePropagator prop(&UGrid), ref(&UGrid);
+    gaussian(rng4, prop);
+
+    // Define variables for sanity checking a single site.
+    typename SpinColourVector::scalar_object fermSite;
+    typename SpinColourMatrix::scalar_object propSite;
+    std::vector<int> site(Nd, 0);
+
+    for (int s = 0; s < Ns; ++s)
+    for (int c = 0; c < Nc; ++c)
+    {
+        ref = prop;
+        PropToFerm(ferm, prop, s, c);
+        FermToProp(prop, ferm, s, c);
+
+        std::cout << "Spin = " << s << ", Colour = " << c << std::endl;
+        ref -= prop;
+        std::cout << "Prop->Ferm->Prop test, diff = " << Grid::sqrt(norm2(ref)) << std::endl;
+
+        peekSite(fermSite, ferm, site);
+        peekSite(propSite, prop, site);
+        for (int s2 = 0; s2 < Ns; ++s2)
+        for (int c2 = 0; c2 < Nc; ++c2)
+        {
+            if (propSite()(s2, s)(c2, c) != fermSite()(s2)(c2))
+            {
+                std::cout << propSite()(s2, s)(c2, c) << " != "
+                          << fermSite()(s2)(c2) << " for spin = " << s2
+                          << ", col = " << c2 << std::endl;
+            }
+        }
+    }
+
+    Grid_finalize();
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons_seq_gamma.cc

@@ -0,0 +1,89 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+
+ Source file: tests/hadrons/Test_hadrons_seq_gamma.cc
+
+ Copyright (C) 2017
+
+ Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include "Test_hadrons.hpp"
+
+using namespace Grid;
+using namespace QCD;
+using namespace Hadrons;
+
+/*******************************************************************************
+ * Consistency test for sequential gamma insertion.
+ ******************************************************************************/
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    HADRONS_DEFAULT_INIT;
+
+    // run setup ///////////////////////////////////////////////////////////////
+    Application  application;
+    unsigned int nt   = GridDefaultLatt()[Tp];
+    unsigned int tS   = nt / 2;
+    unsigned int Ls   = 12;
+    double       mass = 0.04;
+    double       M5   = 1.8;
+
+    // global parameters
+    HADRONS_DEFAULT_GLOBALS(application);
+
+    // gauge field
+    std::string gaugeField = "gauge";
+    application.createModule<MGauge::Unit>(gaugeField);
+
+    // action
+    std::string actionName = "DWF";
+    makeDWFAction(application, actionName, gaugeField, mass, M5, Ls);
+
+    // solver
+    std::string solverName = "CG";
+    makeRBPrecCGSolver(application, solverName, actionName);
+
+    // test sequential propagator, with g5 insertion.
+    Gamma::Algebra g = Gamma::Algebra::Gamma5;
+    std::string pointProp = "q_0";
+    std::string point5d   = LABEL_5D(pointProp);
+    std::string origin    = "0 0 0 0";
+    MAKE_POINT_PROP(origin, pointProp, solverName);
+
+    std::string seqProp = ADD_INDEX(pointProp + "_seqg5", tS);
+    std::string seqSrc  = seqProp + "_src";
+    MAKE_SEQUENTIAL_PROP(tS, pointProp, ZERO_MOM, seqProp, solverName, g);
+
+    std::string modName = "Test g5 sequential insertion";
+    makeSeqGamComparison(application, modName, pointProp, seqProp, origin, g, tS);
+
+    // execution
+    application.saveParameterFile("SeqGamma5Test.xml");
+    application.run();
+
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+
+    return EXIT_SUCCESS;
+}