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

# Conflicts:
#	extras/Hadrons/modules.inc
#	lib/qcd/action/gauge/Photon.h

.travis.yml

@@ -19,6 +19,8 @@ before_install:
     - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install libmpc; fi
     
 install:
+    - export CWD=`pwd`
+    - echo $CWD
     - export CC=$CC$VERSION
     - export CXX=$CXX$VERSION
     - echo $PATH
@@ -36,11 +38,22 @@ script:
     - ./bootstrap.sh
     - mkdir build
     - cd build
-    - ../configure --enable-precision=single --enable-simd=SSE4 --enable-comms=none
+    - mkdir lime
+    - cd lime
+    - mkdir build
+    - cd build
+    - wget http://usqcd-software.github.io/downloads/c-lime/lime-1.3.2.tar.gz
+    - tar xf lime-1.3.2.tar.gz
+    - cd lime-1.3.2
+    - ./configure --prefix=$CWD/build/lime/install
+    - make -j4
+    - make install
+    - cd $CWD/build
+    - ../configure --enable-precision=single --enable-simd=SSE4 --enable-comms=none --with-lime=$CWD/build/lime/install
     - make -j4 
     - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals
     - echo make clean
-    - ../configure --enable-precision=double --enable-simd=SSE4 --enable-comms=none
+    - ../configure --enable-precision=double --enable-simd=SSE4 --enable-comms=none --with-lime=$CWD/build/lime/install
     - make -j4
     - ./benchmarks/Benchmark_dwf --threads 1 --debug-signals
     - make check

VERSION

@@ -1,4 +1,4 @@
-Version : 0.7.0
+Version : 0.8.0
 
 - Clang 3.5 and above, ICPC v16 and above, GCC 6.3 and above recommended
 - MPI and MPI3 comms optimisations for KNL and OPA finished

benchmarks/Benchmark_IO.cc

@@ -0,0 +1,108 @@
+#include <Grid/Grid.h>
+#ifdef HAVE_LIME
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+#define MSG cout << GridLogMessage
+#define SEP \
+"============================================================================="
+#ifndef BENCH_IO_LMAX
+#define BENCH_IO_LMAX 40
+#endif
+
+typedef function<void(const string, LatticeFermion &)> WriterFn;
+typedef function<void(LatticeFermion &, const string)> ReaderFn;
+
+string filestem(const int l)
+{
+  return "iobench_l" + to_string(l);
+}
+
+void limeWrite(const string filestem, LatticeFermion &vec)
+{
+  emptyUserRecord record;
+  ScidacWriter    binWriter(vec._grid->IsBoss());
+
+  binWriter.open(filestem + ".bin");
+  binWriter.writeScidacFieldRecord(vec, record);
+  binWriter.close();
+}
+
+void limeRead(LatticeFermion &vec, const string filestem)
+{
+  emptyUserRecord record;
+  ScidacReader    binReader;
+
+  binReader.open(filestem + ".bin");
+  binReader.readScidacFieldRecord(vec, record);
+  binReader.close();
+}
+
+void writeBenchmark(const int l, const WriterFn &write)
+{
+  auto                      mpi  = GridDefaultMpi();
+  auto                      simd = GridDefaultSimd(Nd, vComplex::Nsimd());
+  vector<int>               latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
+  unique_ptr<GridCartesian> gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
+  GridCartesian             *g = gPt.get();
+  GridParallelRNG           rng(g);
+  LatticeFermion            vec(g);
+  emptyUserRecord           record;
+  ScidacWriter              binWriter(g->IsBoss());
+
+  cout << "-- Local volume " << l << "^4" << endl;
+  random(rng, vec);
+  write(filestem(l), vec);
+}
+
+void readBenchmark(const int l, const ReaderFn &read)
+{
+  auto                      mpi  = GridDefaultMpi();
+  auto                      simd = GridDefaultSimd(Nd, vComplex::Nsimd());
+  vector<int>               latt = {l*mpi[0], l*mpi[1], l*mpi[2], l*mpi[3]};
+  unique_ptr<GridCartesian> gPt(SpaceTimeGrid::makeFourDimGrid(latt, simd, mpi));
+  GridCartesian             *g = gPt.get();
+  LatticeFermion            vec(g);
+  emptyUserRecord           record;
+  ScidacReader              binReader;
+
+  cout << "-- Local volume " << l << "^4" << endl;
+  read(vec, filestem(l));
+}
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  auto simd = GridDefaultSimd(Nd,vComplex::Nsimd());
+  auto mpi  = GridDefaultMpi();
+
+  int64_t threads = GridThread::GetThreads();
+  MSG << "Grid is setup to use " << threads << " threads" << endl;
+  MSG << SEP << endl;
+  MSG << "Benchmark Lime write" << endl;
+  MSG << SEP << endl;
+  for (int l = 4; l <= BENCH_IO_LMAX; l += 2)
+  {
+    writeBenchmark(l, limeWrite);
+  }
+
+  MSG << "Benchmark Lime read" << endl;
+  MSG << SEP << endl;
+  for (int l = 4; l <= BENCH_IO_LMAX; l += 2)
+  {
+    readBenchmark(l, limeRead);
+  }
+
+  Grid_finalize();
+
+  return EXIT_SUCCESS;
+}
+#else
+int main (int argc, char ** argv)
+{
+  return EXIT_SUCCESS;
+}
+#endif

benchmarks/Benchmark_ITT.cc

@@ -158,8 +158,10 @@ public:
 
 	  dbytes=0;
 	  ncomm=0;
-
-	  parallel_for(int dir=0;dir<8;dir++){
+#ifdef GRID_OMP
+#pragma omp parallel for num_threads(Grid::CartesianCommunicator::nCommThreads)
+#endif
+	  for(int dir=0;dir<8;dir++){
 
 	    double tbytes;
 	    int mu =dir % 4;
@@ -175,9 +177,14 @@ public:
 		int comm_proc = mpi_layout[mu]-1;
 		Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
 	      }
+#ifdef GRID_OMP
+	int tid = omp_get_thread_num(); 
+#else 
+        int tid = dir;
+#endif
 	      tbytes= Grid.StencilSendToRecvFrom((void *)&xbuf[dir][0], xmit_to_rank,
 						 (void *)&rbuf[dir][0], recv_from_rank,
-						 bytes,dir);
+						 bytes,tid);
 	  
 #ifdef GRID_OMP
 #pragma omp atomic

benchmarks/Benchmark_comms.cc

@@ -169,7 +169,11 @@ int main (int argc, char ** argv)
   for(int lat=4;lat<=maxlat;lat+=4){
     for(int Ls=8;Ls<=8;Ls*=2){
 
-      std::vector<int> latt_size  ({lat,lat,lat,lat});
+      std::vector<int> latt_size  ({lat*mpi_layout[0],
+                                    lat*mpi_layout[1],
+                                    lat*mpi_layout[2],
+                                    lat*mpi_layout[3]});
+
 
       GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
       RealD Nrank = Grid._Nprocessors;
@@ -446,7 +450,7 @@ int main (int argc, char ** argv)
   }    
 
 
-
+#ifdef GRID_OMP
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "= Benchmarking threaded STENCIL halo exchange in "<<nmu<<" dimensions"<<std::endl;
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
@@ -485,7 +489,8 @@ int main (int argc, char ** argv)
 	dbytes=0;
 	ncomm=0;
 
-	parallel_for(int dir=0;dir<8;dir++){
+#pragma omp parallel for num_threads(Grid::CartesianCommunicator::nCommThreads)
+	for(int dir=0;dir<8;dir++){
 
 	  double tbytes;
 	  int mu =dir % 4;
@@ -502,9 +507,9 @@ int main (int argc, char ** argv)
 	      int comm_proc = mpi_layout[mu]-1;
 	      Grid.ShiftedRanks(mu,comm_proc,xmit_to_rank,recv_from_rank);
 	    }
-
+            int tid = omp_get_thread_num();
 	    tbytes= Grid.StencilSendToRecvFrom((void *)&xbuf[dir][0], xmit_to_rank,
-					       (void *)&rbuf[dir][0], recv_from_rank, bytes,dir);
+					       (void *)&rbuf[dir][0], recv_from_rank, bytes,tid);
 
 #pragma omp atomic
 	    dbytes+=tbytes;
@@ -532,7 +537,7 @@ int main (int argc, char ** argv)
  
     }
   }    
-
+#endif
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "= All done; Bye Bye"<<std::endl;
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;

benchmarks/Benchmark_memory_bandwidth.cc

@@ -55,7 +55,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s"<<"\t\t"<<"Gflop/s"<<"\t\t seconds"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
-  uint64_t lmax=96;
+  uint64_t lmax=64;
 #define NLOOP (10*lmax*lmax*lmax*lmax/vol)
   for(int lat=8;lat<=lmax;lat+=8){
 

benchmarks/Benchmark_su3.cc

@@ -35,9 +35,11 @@ using namespace Grid::QCD;
 int main (int argc, char ** argv)
 {
   Grid_init(&argc,&argv);
-#define LMAX (64)
+#define LMAX (32)
+#define LMIN (16)
+#define LINC (4)
 
-  int64_t Nloop=20;
+  int64_t Nloop=2000;
 
   std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
   std::vector<int> mpi_layout  = GridDefaultMpi();
@@ -51,7 +53,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
-  for(int lat=2;lat<=LMAX;lat+=2){
+  for(int lat=LMIN;lat<=LMAX;lat+=LINC){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
       int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
@@ -83,7 +85,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
-  for(int lat=2;lat<=LMAX;lat+=2){
+  for(int lat=LMIN;lat<=LMAX;lat+=LINC){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
       int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
@@ -114,7 +116,7 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
-  for(int lat=2;lat<=LMAX;lat+=2){
+  for(int lat=LMIN;lat<=LMAX;lat+=LINC){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
       int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
@@ -145,7 +147,38 @@ int main (int argc, char ** argv)
   std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
   std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
-  for(int lat=2;lat<=LMAX;lat+=2){
+  for(int lat=LMIN;lat<=LMAX;lat+=LINC){
+    
+    std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
+    int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
+    
+    GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
+    GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+    
+    LatticeColourMatrix z(&Grid); random(pRNG,z);
+    LatticeColourMatrix x(&Grid); random(pRNG,x);
+    LatticeColourMatrix y(&Grid); random(pRNG,y);
+    
+    double start=usecond();
+    for(int64_t i=0;i<Nloop;i++){
+      mac(z,x,y);
+    }
+    double stop=usecond();
+    double time = (stop-start)/Nloop*1000.0;
+    
+    double bytes=3*vol*Nc*Nc*sizeof(Complex);
+    double flops=Nc*Nc*(6+8+8)*vol;
+    std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
+    
+  }
+
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "= Benchmarking SU3xSU3  CovShiftForward(z,x,y)"<<std::endl;
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
+  std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
+
+  for(int lat=LMIN;lat<=LMAX;lat+=LINC){
 
       std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
       int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
@@ -157,18 +190,64 @@ int main (int argc, char ** argv)
       LatticeColourMatrix x(&Grid); random(pRNG,x);
       LatticeColourMatrix y(&Grid); random(pRNG,y);
 
-      double start=usecond();
-      for(int64_t i=0;i<Nloop;i++){
-	mac(z,x,y);
+      for(int mu=0;mu<4;mu++){
+	      double start=usecond();
+	      for(int64_t i=0;i<Nloop;i++){
+	        z = PeriodicBC::CovShiftForward(x,mu,y);
+	    }
+	    double stop=usecond();
+	    double time = (stop-start)/Nloop*1000.0;
+	
+	
+	    double bytes=3*vol*Nc*Nc*sizeof(Complex);
+	    double flops=Nc*Nc*(6+8+8)*vol;
+	    std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
       }
-      double stop=usecond();
-      double time = (stop-start)/Nloop*1000.0;
-      
-      double bytes=3*vol*Nc*Nc*sizeof(Complex);
-      double flops=Nc*Nc*(8+8+8)*vol;
-      std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
+  }
+#if 1
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "= Benchmarking SU3xSU3  z= x * Cshift(y)"<<std::endl;
+  std::cout<<GridLogMessage << "===================================================================================================="<<std::endl;
+  std::cout<<GridLogMessage << "  L  "<<"\t\t"<<"bytes"<<"\t\t\t"<<"GB/s\t\t GFlop/s"<<std::endl;
+  std::cout<<GridLogMessage << "----------------------------------------------------------"<<std::endl;
 
+  for(int lat=LMIN;lat<=LMAX;lat+=LINC){
+      std::vector<int> latt_size  ({lat*mpi_layout[0],lat*mpi_layout[1],lat*mpi_layout[2],lat*mpi_layout[3]});
+      int64_t vol = latt_size[0]*latt_size[1]*latt_size[2]*latt_size[3];
+
+      GridCartesian     Grid(latt_size,simd_layout,mpi_layout);
+      GridParallelRNG          pRNG(&Grid);      pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
+
+      LatticeColourMatrix z(&Grid); random(pRNG,z);
+      LatticeColourMatrix x(&Grid); random(pRNG,x);
+      LatticeColourMatrix y(&Grid); random(pRNG,y);
+      LatticeColourMatrix tmp(&Grid);
+
+      for(int mu=0;mu<4;mu++){
+	double tshift=0;
+	double tmult =0;
+
+	double start=usecond();
+	for(int64_t i=0;i<Nloop;i++){
+	  tshift-=usecond();
+	  tmp = Cshift(y,mu,-1);
+	  tshift+=usecond();
+	  tmult-=usecond();
+	  z   = x*tmp;
+	  tmult+=usecond();
+	}
+	double stop=usecond();
+	double time = (stop-start)/Nloop;
+	tshift = tshift/Nloop;
+	tmult  = tmult /Nloop;
+	
+	double bytes=3*vol*Nc*Nc*sizeof(Complex);
+	double flops=Nc*Nc*(6+8+8)*vol;
+	std::cout<<GridLogMessage<<std::setprecision(3) << "total us "<<time<<" shift "<<tshift <<" mult "<<tmult<<std::endl;
+	time = time * 1000; // convert to NS for GB/s
+	std::cout<<GridLogMessage<<std::setprecision(3) << lat<<"\t\t"<<bytes<<"   \t\t"<<bytes/time<<"\t\t" << flops/time<<std::endl;
+      }
     }
-
+#endif
   Grid_finalize();
 }

configure.ac

@@ -340,7 +340,7 @@ case ${ac_PRECISION} in
 esac
 
 ######################  Shared memory allocation technique under MPI3
-AC_ARG_ENABLE([shm],[AC_HELP_STRING([--enable-shm=shmopen|hugetlbfs|shmnone],
+AC_ARG_ENABLE([shm],[AC_HELP_STRING([--enable-shm=shmopen|shmget|hugetlbfs|shmnone],
               [Select SHM allocation technique])],[ac_SHM=${enable_shm}],[ac_SHM=shmopen])
 
 case ${ac_SHM} in
@@ -349,6 +349,10 @@ case ${ac_SHM} in
      AC_DEFINE([GRID_MPI3_SHMOPEN],[1],[GRID_MPI3_SHMOPEN] )
      ;;
 
+     shmget)
+     AC_DEFINE([GRID_MPI3_SHMGET],[1],[GRID_MPI3_SHMGET] )
+     ;;
+
      shmnone)
      AC_DEFINE([GRID_MPI3_SHM_NONE],[1],[GRID_MPI3_SHM_NONE] )
      ;;
@@ -366,7 +370,7 @@ esac
 AC_ARG_ENABLE([shmpath],[AC_HELP_STRING([--enable-shmpath=path],
               [Select SHM mmap base path for hugetlbfs])],
 	      [ac_SHMPATH=${enable_shmpath}],
-	      [ac_SHMPATH=/var/lib/hugetlbfs/pagesize-2MB/])
+	      [ac_SHMPATH=/var/lib/hugetlbfs/global/pagesize-2MB/])
 AC_DEFINE_UNQUOTED([GRID_SHM_PATH],["$ac_SHMPATH"],[Path to a hugetlbfs filesystem for MMAPing])
 
 ############### communication type selection

extras/Hadrons/Application.cc

@@ -28,6 +28,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 
 #include <Grid/Hadrons/Application.hpp>
 #include <Grid/Hadrons/GeneticScheduler.hpp>
+#include <Grid/Hadrons/Modules.hpp>
 
 using namespace Grid;
 using namespace QCD;
@@ -113,12 +114,12 @@ void Application::parseParameterFile(const std::string parameterFileName)
     setPar(par);
     if (!push(reader, "modules"))
     {
-        HADRON_ERROR(Parsing, "Cannot open node 'modules' in parameter file '" 
+        HADRONS_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 '" 
+        HADRONS_ERROR(Parsing, "Cannot open node 'modules/module' in parameter file '" 
                               + parameterFileName + "'");
     }
     do
@@ -176,7 +177,7 @@ void Application::saveSchedule(const std::string filename)
         
         if (!scheduled_)
         {
-            HADRON_ERROR(Definition, "Computation not scheduled");
+            HADRONS_ERROR(Definition, "Computation not scheduled");
         }
 
         for (auto address: program_)
@@ -207,7 +208,7 @@ void Application::printSchedule(void)
 {
     if (!scheduled_)
     {
-        HADRON_ERROR(Definition, "Computation not scheduled");
+        HADRONS_ERROR(Definition, "Computation not scheduled");
     }
     auto peak = vm().memoryNeeded(program_);
     LOG(Message) << "Schedule (memory needed: " << sizeString(peak) << "):"

extras/Hadrons/Application.hpp

@@ -31,7 +31,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/VirtualMachine.hpp>
-#include <Grid/Hadrons/Modules.hpp>
+#include <Grid/Hadrons/Module.hpp>
 
 BEGIN_HADRONS_NAMESPACE
 

extras/Hadrons/EigenPack.hpp

@@ -44,9 +44,21 @@ class EigenPack
 {
 public:
     typedef F Field;
+    struct PackRecord
+    {
+        std::string operatorXml, solverXml;
+    };
+    struct VecRecord: Serializable
+    {
+        GRID_SERIALIZABLE_CLASS_MEMBERS(VecRecord,
+                                        unsigned int, index,
+                                        double,       eval);
+        VecRecord(void): index(0), eval(0.) {}
+    };
 public:
     std::vector<RealD> eval;
     std::vector<F>     evec;
+    PackRecord         record;
 public:
     EigenPack(void)          = default;
     virtual ~EigenPack(void) = default;
@@ -62,68 +74,139 @@ public:
         evec.resize(size, grid);
     }
 
-    virtual void read(const std::string fileStem, const int traj = -1)
+    virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        std::string evecFilename, evalFilename;
-
-        makeFilenames(evecFilename, evalFilename, fileStem, traj);
-        XmlReader xmlReader(evalFilename);
-        basicRead(evec, evecFilename, evec.size());
-        LOG(Message) << "Reading " << eval.size() << " eigenvalues from '" 
-                     << evalFilename << "'" << std::endl;
-        Grid::read(xmlReader, "evals", eval);
+        if (multiFile)
+        {
+            for(int k = 0; k < evec.size(); ++k)
+            {
+                basicReadSingle(evec[k], eval[k], evecFilename(fileStem, k, traj), k);
+            }
+        }
+        else
+        {
+            basicRead(evec, eval, evecFilename(fileStem, -1, traj), evec.size());
+        }
     }
 
-    virtual void write(const std::string fileStem, const int traj = -1)
+    virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        std::string evecFilename, evalFilename;
-
-        makeFilenames(evecFilename, evalFilename, fileStem, traj);
-        XmlWriter xmlWriter(evalFilename);
-        basicWrite(evecFilename, evec, evec.size());
-        LOG(Message) << "Writing " << eval.size() << " eigenvalues to '" 
-                     << evalFilename << "'" << std::endl;
-        Grid::write(xmlWriter, "evals", eval);
+        if (multiFile)
+        {
+            for(int k = 0; k < evec.size(); ++k)
+            {
+                basicWriteSingle(evecFilename(fileStem, k, traj), evec[k], eval[k], k);
+            }
+        }
+        else
+        {
+            basicWrite(evecFilename(fileStem, -1, traj), evec, eval, evec.size());
+        }
     }
 protected:
-    void makeFilenames(std::string &evecFilename, std::string &evalFilename,
-                       const std::string stem, const int traj = -1)
+    std::string evecFilename(const std::string stem, const int vec, const int traj)
     {
         std::string t = (traj < 0) ? "" : ("." + std::to_string(traj));
 
-        evecFilename = stem + "_evec" + t + ".bin";
-        evalFilename = stem + "_eval" + t + ".xml";
+        if (vec == -1)
+        {
+            return stem + t + ".bin";
+        }
+        else
+        {
+            return stem + t + "/v" + std::to_string(vec) + ".bin";
+        }
     }
 
     template <typename T>
-    static void basicRead(std::vector<T> &evec, const std::string filename, 
-                          const unsigned int size)
+    void basicRead(std::vector<T> &evec, std::vector<double> &eval,
+                   const std::string filename, const unsigned int size)
     {
-        emptyUserRecord record;
         ScidacReader    binReader;
 
         binReader.open(filename);
+        binReader.skipPastObjectRecord(SCIDAC_FILE_XML);
         for(int k = 0; k < size; ++k) 
         {
-            binReader.readScidacFieldRecord(evec[k], record);
+            VecRecord vecRecord;
+
+            LOG(Message) << "Reading eigenvector " << k << std::endl;
+            binReader.readScidacFieldRecord(evec[k], vecRecord);
+            if (vecRecord.index != k)
+            {
+                HADRONS_ERROR(Io, "Eigenvector " + std::to_string(k) + " has a"
+                              + " wrong index (expected " + std::to_string(vecRecord.index) 
+                              + ") in file '" + filename + "'");
+            }
+            eval[k] = vecRecord.eval;
         }
         binReader.close();
     }
 
     template <typename T>
-    static void basicWrite(const std::string filename, std::vector<T> &evec, 
-                           const unsigned int size)
+    void basicReadSingle(T &evec, double &eval, const std::string filename, 
+                         const unsigned int index)
     {
-        emptyUserRecord record;
-        ScidacWriter    binWriter;
+        ScidacReader binReader;
+        VecRecord    vecRecord;
 
+        binReader.open(filename);
+        binReader.skipPastObjectRecord(SCIDAC_FILE_XML);
+        LOG(Message) << "Reading eigenvector " << index << std::endl;
+        binReader.readScidacFieldRecord(evec, vecRecord);
+        if (vecRecord.index != index)
+        {
+            HADRONS_ERROR(Io, "Eigenvector " + std::to_string(index) + " has a"
+                          + " wrong index (expected " + std::to_string(vecRecord.index) 
+                          + ") in file '" + filename + "'");
+        }
+        eval = vecRecord.eval;
+        binReader.close();
+    }
+
+    template <typename T>
+    void basicWrite(const std::string filename, std::vector<T> &evec, 
+                    const std::vector<double> &eval, const unsigned int size)
+    {
+        ScidacWriter binWriter(evec[0]._grid->IsBoss());
+        XmlWriter    xmlWriter("", "eigenPackPar");
+
+        makeFileDir(filename, evec[0]._grid);
+        xmlWriter.pushXmlString(record.operatorXml);
+        xmlWriter.pushXmlString(record.solverXml);
         binWriter.open(filename);
+        binWriter.writeLimeObject(1, 1, xmlWriter, "parameters", SCIDAC_FILE_XML);
         for(int k = 0; k < size; ++k) 
         {
-            binWriter.writeScidacFieldRecord(evec[k], record);
+            VecRecord vecRecord;
+
+            vecRecord.index = k;
+            vecRecord.eval  = eval[k];
+            LOG(Message) << "Writing eigenvector " << k << std::endl;
+            binWriter.writeScidacFieldRecord(evec[k], vecRecord, DEFAULT_ASCII_PREC);
         }
         binWriter.close();
     }
+
+    template <typename T>
+    void basicWriteSingle(const std::string filename, T &evec, 
+                          const double eval, const unsigned int index)
+    {
+        ScidacWriter binWriter(evec._grid->IsBoss());
+        XmlWriter    xmlWriter("", "eigenPackPar");
+        VecRecord    vecRecord;
+
+        makeFileDir(filename, evec._grid);
+        xmlWriter.pushXmlString(record.operatorXml);
+        xmlWriter.pushXmlString(record.solverXml);
+        binWriter.open(filename);
+        binWriter.writeLimeObject(1, 1, xmlWriter, "parameters", SCIDAC_FILE_XML);
+        vecRecord.index = index;
+        vecRecord.eval  = eval;
+        LOG(Message) << "Writing eigenvector " << index << std::endl;
+        binWriter.writeScidacFieldRecord(evec, vecRecord, DEFAULT_ASCII_PREC);
+        binWriter.close();
+    }
 };
 
 template <typename FineF, typename CoarseF>
@@ -152,54 +235,76 @@ public:
         evecCoarse.resize(sizeCoarse, gridCoarse);
     }
 
-    virtual void read(const std::string fileStem, const int traj = -1)
+    void readFine(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        std::string evecFineFilename, evalFineFilename;
-        std::string evecCoarseFilename, evalCoarseFilename;
-
-        this->makeFilenames(evecFineFilename, evalFineFilename, 
-                            fileStem + "_fine", traj);
-        this->makeFilenames(evecCoarseFilename, evalCoarseFilename, 
-                            fileStem + "_coarse", traj);
-        XmlReader xmlFineReader(evalFineFilename);
-        XmlReader xmlCoarseReader(evalCoarseFilename);
-        LOG(Message) << "Reading " << this->evec.size() << " fine eigenvectors from '" 
-                     << evecFineFilename << "'" << std::endl;
-        this->basicRead(this->evec, evecFineFilename, this->evec.size());
-        LOG(Message) << "Reading " << evecCoarse.size() << " coarse eigenvectors from '" 
-                     << evecCoarseFilename << "'" << std::endl;
-        this->basicRead(evecCoarse, evecCoarseFilename, evecCoarse.size());
-        LOG(Message) << "Reading " << this->eval.size() << " fine eigenvalues from '" 
-                     << evalFineFilename << "'" << std::endl;
-        Grid::read(xmlFineReader, "evals", this->eval);
-        LOG(Message) << "Reading " << evalCoarse.size() << " coarse eigenvalues from '" 
-                     << evalCoarseFilename << "'" << std::endl;
-        Grid::read(xmlCoarseReader, "evals", evalCoarse);
+        if (multiFile)
+        {
+            for(int k = 0; k < this->evec.size(); ++k)
+            {
+                this->basicReadSingle(this->evec[k], this->eval[k], this->evecFilename(fileStem + "_fine", k, traj), k);
+            }
+        }
+        else
+        {
+            this->basicRead(this->evec, this->eval, this->evecFilename(fileStem + "_fine", -1, traj), this->evec.size());
+        }
     }
 
-    virtual void write(const std::string fileStem, const int traj = -1)
+    void readCoarse(const std::string fileStem, const bool multiFile, const int traj = -1)
     {
-        std::string evecFineFilename, evalFineFilename;
-        std::string evecCoarseFilename, evalCoarseFilename;
+        if (multiFile)
+        {
+            for(int k = 0; k < evecCoarse.size(); ++k)
+            {
+                this->basicReadSingle(evecCoarse[k], evalCoarse[k], this->evecFilename(fileStem + "_coarse", k, traj), k);
+            }
+        }
+        else
+        {
+            this->basicRead(evecCoarse, evalCoarse, this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse.size());
+        }
+    }
 
-        this->makeFilenames(evecFineFilename, evalFineFilename, 
-                            fileStem + "_fine", traj);
-        this->makeFilenames(evecCoarseFilename, evalCoarseFilename,
-                            fileStem + "_coarse", traj);
-        XmlWriter xmlFineWriter(evalFineFilename);
-        XmlWriter xmlCoarseWriter(evalCoarseFilename);
-        LOG(Message) << "Writing " << this->evec.size() << " fine eigenvectors to '" 
-                     << evecFineFilename << "'" << std::endl;
-        this->basicWrite(evecFineFilename, this->evec, this->evec.size());
-        LOG(Message) << "Writing " << evecCoarse.size() << " coarse eigenvectors to '" 
-                     << evecCoarseFilename << "'" << std::endl;
-        this->basicWrite(evecCoarseFilename, evecCoarse, evecCoarse.size());
-        LOG(Message) << "Writing " << this->eval.size() << " fine eigenvalues to '" 
-                     << evalFineFilename << "'" << std::endl;
-        Grid::write(xmlFineWriter, "evals", this->eval);
-        LOG(Message) << "Writing " << evalCoarse.size() << " coarse eigenvalues to '" 
-                     << evalCoarseFilename << "'" << std::endl;
-        Grid::write(xmlCoarseWriter, "evals", evalCoarse);
+    virtual void read(const std::string fileStem, const bool multiFile, const int traj = -1)
+    {
+        readFine(fileStem, multiFile, traj);
+        readCoarse(fileStem, multiFile, traj);
+    }
+
+    void writeFine(const std::string fileStem, const bool multiFile, const int traj = -1)
+    {
+        if (multiFile)
+        {
+            for(int k = 0; k < this->evec.size(); ++k)
+            {
+                this->basicWriteSingle(this->evecFilename(fileStem + "_fine", k, traj), this->evec[k], this->eval[k], k);
+            }
+        }
+        else
+        {
+            this->basicWrite(this->evecFilename(fileStem + "_fine", -1, traj), this->evec, this->eval, this->evec.size());
+        }
+    }
+
+    void writeCoarse(const std::string fileStem, const bool multiFile, const int traj = -1)
+    {
+        if (multiFile)
+        {
+            for(int k = 0; k < evecCoarse.size(); ++k)
+            {
+                this->basicWriteSingle(this->evecFilename(fileStem + "_coarse", k, traj), evecCoarse[k], evalCoarse[k], k);
+            }
+        }
+        else
+        {
+            this->basicWrite(this->evecFilename(fileStem + "_coarse", -1, traj), evecCoarse, evalCoarse, evecCoarse.size());
+        }
+    }
+    
+    virtual void write(const std::string fileStem, const bool multiFile, const int traj = -1)
+    {
+        writeFine(fileStem, multiFile, traj);
+        writeCoarse(fileStem, multiFile, traj);
     }
 };
 

extras/Hadrons/Environment.cc

@@ -35,7 +35,7 @@ using namespace QCD;
 using namespace Hadrons;
 
 #define ERROR_NO_ADDRESS(address)\
-HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
+HADRONS_ERROR(Definition, "no object with address " + std::to_string(address));
 
 /******************************************************************************
  *                       Environment implementation                           *
@@ -85,7 +85,7 @@ void Environment::createCoarseGrid(const std::vector<int> &blockSize,
         coarseDim[d] = fineDim[d]/blockSize[d];
         if (coarseDim[d]*blockSize[d] != fineDim[d])
         {
-            HADRON_ERROR(Size, "Fine dimension " + std::to_string(d) 
+            HADRONS_ERROR(Size, "Fine dimension " + std::to_string(d) 
                          + " (" + std::to_string(fineDim[d]) 
                          + ") not divisible by coarse dimension ("
                          + std::to_string(coarseDim[d]) + ")"); 
@@ -96,7 +96,7 @@ void Environment::createCoarseGrid(const std::vector<int> &blockSize,
         cLs = Ls/blockSize[nd];
         if (cLs*blockSize[nd] != Ls)
         {
-            HADRON_ERROR(Size, "Fine Ls (" + std::to_string(Ls) 
+            HADRONS_ERROR(Size, "Fine Ls (" + std::to_string(Ls) 
                          + ") not divisible by coarse Ls ("
                          + std::to_string(cLs) + ")");
         }
@@ -106,14 +106,10 @@ void Environment::createCoarseGrid(const std::vector<int> &blockSize,
     gridCoarse4d_[key4d].reset(
         SpaceTimeGrid::makeFourDimGrid(coarseDim, 
             GridDefaultSimd(nd, vComplex::Nsimd()), GridDefaultMpi()));
-    gridCoarseRb4d_[key4d].reset(
-        SpaceTimeGrid::makeFourDimRedBlackGrid(gridCoarse4d_[key4d].get()));
     if (Ls > 1)
     {
         gridCoarse5d_[key5d].reset(
             SpaceTimeGrid::makeFiveDimGrid(cLs, gridCoarse4d_[key4d].get()));
-        gridCoarseRb5d_[key5d].reset(
-            SpaceTimeGrid::makeFiveDimRedBlackGrid(cLs, gridCoarse4d_[key4d].get()));
     }
 }
 
@@ -132,7 +128,7 @@ GridCartesian * Environment::getGrid(const unsigned int Ls) const
     }
     catch(std::out_of_range &)
     {
-        HADRON_ERROR(Definition, "no grid with Ls= " + std::to_string(Ls));
+        HADRONS_ERROR(Definition, "no grid with Ls= " + std::to_string(Ls));
     }
 }
 
@@ -151,7 +147,7 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
     }
     catch(std::out_of_range &)
     {
-        HADRON_ERROR(Definition, "no red-black grid with Ls= " + std::to_string(Ls));
+        HADRONS_ERROR(Definition, "no red-black grid with Ls= " + std::to_string(Ls));
     }
 }
 
@@ -175,31 +171,7 @@ GridCartesian * Environment::getCoarseGrid(
     }
     catch(std::out_of_range &)
     {
-        HADRON_ERROR(Definition, "no coarse grid with Ls= " + std::to_string(Ls));
-    }
-}
-
-GridRedBlackCartesian * Environment::getRbCoarseGrid(
-    const std::vector<int> &blockSize, const unsigned int Ls) const
-{
-    auto key = blockSize;
-
-    try
-    {
-        if (Ls == 1)
-        {
-            key.resize(getNd());
-            return gridCoarseRb4d_.at(key).get();
-        }
-        else
-        {
-            key.push_back(Ls);
-            return gridCoarseRb5d_.at(key).get();
-        }
-    }
-    catch(std::out_of_range &)
-    {
-        HADRON_ERROR(Definition, "no coarse red-black grid with Ls= " + std::to_string(Ls));
+        HADRONS_ERROR(Definition, "no coarse grid with Ls= " + std::to_string(Ls));
     }
 }
 
@@ -249,7 +221,7 @@ void Environment::addObject(const std::string name, const int moduleAddress)
     }
     else
     {
-        HADRON_ERROR(Definition, "object '" + name + "' already exists");
+        HADRONS_ERROR(Definition, "object '" + name + "' already exists");
     }
 }
 
@@ -272,7 +244,7 @@ unsigned int Environment::getObjectAddress(const std::string name) const
     }
     else
     {
-        HADRON_ERROR(Definition, "no object with name '" + name + "'");
+        HADRONS_ERROR(Definition, "no object with name '" + name + "'");
     }
 }
 

extras/Hadrons/Environment.hpp

@@ -92,8 +92,6 @@ public:
     GridRedBlackCartesian * getRbGrid(const unsigned int Ls = 1) const;
     GridCartesian *         getCoarseGrid(const std::vector<int> &blockSize,
                                           const unsigned int Ls = 1) const;
-    GridRedBlackCartesian * getRbCoarseGrid(const std::vector<int> &blockSize,
-                                            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;
@@ -166,9 +164,7 @@ private:
     GridRbPt                               gridRb4d_;
     std::map<unsigned int, GridRbPt>       gridRb5d_;
     std::map<std::vector<int>, GridPt>     gridCoarse4d_;
-    std::map<std::vector<int>, GridRbPt>   gridCoarseRb4d_;
     std::map<std::vector<int>, GridPt>     gridCoarse5d_;
-    std::map<std::vector<int>, GridRbPt>   gridCoarseRb5d_;
     unsigned int                           nd_;
     // random number generator
     RngPt                                  rng4d_;
@@ -249,7 +245,7 @@ void Environment::createDerivedObject(const std::string name,
              (object_[address].type        != &typeid(B))     or
              (object_[address].derivedType != &typeid(T)))
     {
-        HADRON_ERROR(Definition, "object '" + name + "' already allocated");
+        HADRONS_ERROR(Definition, "object '" + name + "' already allocated");
     }
 }
 
@@ -283,7 +279,7 @@ T * Environment::getDerivedObject(const unsigned int address) const
                     }
                     else
                     {
-                        HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
+                        HADRONS_ERROR(Definition, "object with address " + std::to_string(address) +
                             " cannot be casted to '" + typeName(&typeid(T)) +
                             "' (has type '" + typeName(&typeid(h->get())) + "')");
                     }
@@ -291,20 +287,20 @@ T * Environment::getDerivedObject(const unsigned int address) const
             }
             else
             {
-                HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
+                HADRONS_ERROR(Definition, "object with address " + std::to_string(address) +
                             " does not have type '" + typeName(&typeid(B)) +
                             "' (has type '" + getObjectType(address) + "')");
             }
         }
         else
         {
-            HADRON_ERROR(Definition, "object with address " + std::to_string(address) +
+            HADRONS_ERROR(Definition, "object with address " + std::to_string(address) +
                          " is empty");
         }
     }
     else
     {
-        HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
+        HADRONS_ERROR(Definition, "no object with address " + std::to_string(address));
     }
 }
 
@@ -342,7 +338,7 @@ bool Environment::isObjectOfType(const unsigned int address) const
     }
     else
     {
-        HADRON_ERROR(Definition, "no object with address " + std::to_string(address));
+        HADRONS_ERROR(Definition, "no object with address " + std::to_string(address));
     }
 }
 

extras/Hadrons/Exceptions.cc

@@ -27,6 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 /*  END LEGAL */
 
 #include <Grid/Hadrons/Exceptions.hpp>
+#include <Grid/Hadrons/VirtualMachine.hpp>
+#include <Grid/Hadrons/Module.hpp>
 
 #ifndef ERR_SUFF
 #define ERR_SUFF " (" + loc + ")"
@@ -47,6 +49,7 @@ 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))
@@ -55,3 +58,24 @@ 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))
+
+// abort functions
+void Grid::Hadrons::Exceptions::abort(const std::exception& e)
+{
+    auto &vm = VirtualMachine::getInstance();
+    int  mod = vm.getCurrentModule();
+
+    LOG(Error) << "FATAL ERROR -- Exception " << typeName(&typeid(e)) 
+               << std::endl;
+    if (mod >= 0)
+    {
+        LOG(Error) << "During execution of module '"
+                    << vm.getModuleName(mod) << "' (address " << mod << ")"
+                    << std::endl;
+    }
+    LOG(Error) << e.what() << std::endl;
+    LOG(Error) << "Aborting program" << std::endl;
+    Grid_finalize();
+
+    exit(EXIT_FAILURE);
+}

extras/Hadrons/Exceptions.hpp

@@ -34,10 +34,10 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #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)\
-throw(Exceptions::exc(msg, SRC_LOC));
+#define HADRONS_SRC_LOC std::string(__FUNCTION__) + " at " \
+                        + std::string(__FILE__) + ":" + std::to_string(__LINE__)
+#define HADRONS_ERROR(exc, msg)\
+throw(Exceptions::exc(msg, HADRONS_SRC_LOC));
 
 #define DECL_EXC(name, base) \
 class name: public base\
@@ -56,6 +56,7 @@ namespace Exceptions
     DECL_EXC(Implementation, Logic);
     DECL_EXC(Range, Logic);
     DECL_EXC(Size, Logic);
+
     // runtime errors
     DECL_EXC(Runtime, std::runtime_error);
     DECL_EXC(Argument, Runtime);
@@ -64,6 +65,9 @@ namespace Exceptions
     DECL_EXC(Parsing, Runtime);
     DECL_EXC(Program, Runtime);
     DECL_EXC(System, Runtime);
+
+    // abort functions
+    void abort(const std::exception& e);
 }
 
 END_HADRONS_NAMESPACE

extras/Hadrons/Factory.hpp

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

extras/Hadrons/GeneticScheduler.hpp

@@ -57,7 +57,9 @@ public:
     virtual ~GeneticScheduler(void) = default;
     // access
     const Gene & getMinSchedule(void);
-    int          getMinValue(void);
+    V            getMinValue(void);
+    // reset population
+    void initPopulation(void);
     // breed a new generation
     void nextGeneration(void);
     // heuristic benchmarks
@@ -76,8 +78,6 @@ public:
         return out;
     }
 private:
-    // evolution steps
-    void initPopulation(void);
     void doCrossover(void);
     void doMutation(void);
     // genetic operators
@@ -116,7 +116,7 @@ GeneticScheduler<V, T>::getMinSchedule(void)
 }
 
 template <typename V, typename T>
-int GeneticScheduler<V, T>::getMinValue(void)
+V GeneticScheduler<V, T>::getMinValue(void)
 {
     return population_.begin()->first;
 }

extras/Hadrons/Global.cc

@@ -57,8 +57,8 @@ void Hadrons::initLogger(void)
     GridLogIterative.setChanWidth(cw);
     GridLogDebug.setChanWidth(cw);
     GridLogIRL.setChanWidth(cw);
-    HadronsLogError.Active(GridLogError.isActive());
-    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogError.Active(true);
+    HadronsLogWarning.Active(true);
     HadronsLogMessage.Active(GridLogMessage.isActive());
     HadronsLogIterative.Active(GridLogIterative.isActive());
     HadronsLogDebug.Active(GridLogDebug.isActive());
@@ -92,3 +92,84 @@ const std::string Hadrons::resultFileExt = "h5";
 #else
 const std::string Hadrons::resultFileExt = "xml";
 #endif
+
+// recursive mkdir /////////////////////////////////////////////////////////////
+int Hadrons::mkdir(const std::string dirName)
+{
+    if (!dirName.empty() and access(dirName.c_str(), R_OK|W_OK|X_OK))
+    {
+        mode_t mode755;
+        char   tmp[MAX_PATH_LENGTH];
+        char   *p = NULL;
+        size_t len;
+
+        mode755 = S_IRWXU|S_IRGRP|S_IXGRP|S_IROTH|S_IXOTH;
+
+        snprintf(tmp, sizeof(tmp), "%s", dirName.c_str());
+        len = strlen(tmp);
+        if(tmp[len - 1] == '/')
+        {
+            tmp[len - 1] = 0;
+        }
+        for(p = tmp + 1; *p; p++)
+        {
+            if(*p == '/')
+            {
+                *p = 0;
+                ::mkdir(tmp, mode755);
+                *p = '/';
+            }
+        }
+
+        return ::mkdir(tmp, mode755);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+std::string Hadrons::basename(const std::string &s)
+{
+    constexpr char sep = '/';
+    size_t         i   = s.rfind(sep, s.length());
+    
+    if (i != std::string::npos)
+    {
+        return s.substr(i+1, s.length() - i);
+    }
+    else
+    {
+        return s;
+    }
+}
+
+std::string Hadrons::dirname(const std::string &s)
+{
+    constexpr char sep = '/';
+    size_t         i   = s.rfind(sep, s.length());
+    
+    if (i != std::string::npos)
+    {
+        return s.substr(0, i);
+    }
+    else
+    {
+        return "";
+    }
+}
+
+void Hadrons::makeFileDir(const std::string filename, GridBase *g)
+{
+    if (g->IsBoss())
+    {
+        std::string dir    = dirname(filename);
+        int         status = mkdir(dir);
+
+        if (status)
+        {
+            HADRONS_ERROR(Io, "cannot create directory '" + dir
+                          + "' ( " + std::strerror(errno) + ")");
+        }
+    }
+}

extras/Hadrons/Global.hpp

@@ -39,6 +39,10 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define SITE_SIZE_TYPE size_t
 #endif
 
+#ifndef DEFAULT_ASCII_PREC
+#define DEFAULT_ASCII_PREC 16
+#endif
+
 #define BEGIN_HADRONS_NAMESPACE \
 namespace Grid {\
 using namespace QCD;\
@@ -110,7 +114,7 @@ public:
 };
 
 #define LOG(channel) std::cout << HadronsLog##channel
-#define DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl;
+#define HADRONS_DEBUG_VAR(var) LOG(Debug) << #var << "= " << (var) << std::endl;
 
 extern HadronsLogger HadronsLogError;
 extern HadronsLogger HadronsLogWarning;
@@ -187,8 +191,14 @@ typedef XmlWriter ResultWriter;
 #define RESULT_FILE_NAME(name) \
 name + "." + std::to_string(vm().getTrajectory()) + "." + resultFileExt
 
-// default Schur convention
+// recursive mkdir
+#define MAX_PATH_LENGTH 512u
+int         mkdir(const std::string dirName);
+std::string basename(const std::string &s);
+std::string dirname(const std::string &s);
+void        makeFileDir(const std::string filename, GridBase *g);
 
+// default Schur convention
 #ifndef HADRONS_DEFAULT_SCHUR 
 #define HADRONS_DEFAULT_SCHUR DiagTwo
 #endif

extras/Hadrons/Graph.hpp

@@ -184,7 +184,7 @@ void Graph<T>::removeVertex(const T &value)
     }
     else
     {
-        HADRON_ERROR(Range, "vertex does not exists");
+        HADRONS_ERROR(Range, "vertex does not exists");
     }
 
     // remove all edges containing the vertex
@@ -213,7 +213,7 @@ void Graph<T>::removeEdge(const Edge &e)
     }
     else
     {
-        HADRON_ERROR(Range, "edge does not exists");
+        HADRONS_ERROR(Range, "edge does not exists");
     }
 }
 
@@ -259,7 +259,7 @@ void Graph<T>::mark(const T &value, const bool doMark)
     }
     else
     {
-        HADRON_ERROR(Range, "vertex does not exists");
+        HADRONS_ERROR(Range, "vertex does not exists");
     }
 }
 
@@ -297,7 +297,7 @@ bool Graph<T>::isMarked(const T &value) const
     }
     else
     {
-        HADRON_ERROR(Range, "vertex does not exists");
+        HADRONS_ERROR(Range, "vertex does not exists");
         
         return false;
     }
@@ -543,7 +543,7 @@ std::vector<T> Graph<T>::topoSort(void)
     {
         if (tmpMarked.at(v))
         {
-            HADRON_ERROR(Range, "cannot topologically sort a cyclic graph");
+            HADRONS_ERROR(Range, "cannot topologically sort a cyclic graph");
         }
         if (!isMarked(v))
         {
@@ -602,7 +602,7 @@ std::vector<T> Graph<T>::topoSort(Gen &gen)
     {
         if (tmpMarked.at(v))
         {
-            HADRON_ERROR(Range, "cannot topologically sort a cyclic graph");
+            HADRONS_ERROR(Range, "cannot topologically sort a cyclic graph");
         }
         if (!isMarked(v))
         {

extras/Hadrons/HadronsXmlRun.cc

@@ -69,12 +69,7 @@ int main(int argc, char *argv[])
     }
     catch (const std::exception& e)
     {
-        LOG(Error) << "FATAL ERROR -- Exception " << typeName(&typeid(e)) << std::endl;
-        LOG(Error) << e.what() << std::endl;
-        LOG(Error) << "Aborting program" << std::endl;
-        Grid_finalize();
-
-        return EXIT_FAILURE;
+        Exceptions::abort(e);
     }
     
     // epilogue

extras/Hadrons/Module.cc

@@ -49,7 +49,7 @@ std::string ModuleBase::getName(void) const
 // get factory registration name if available
 std::string ModuleBase::getRegisteredName(void)
 {
-    HADRON_ERROR(Definition, "module '" + getName() + "' has no registered type"
+    HADRONS_ERROR(Definition, "module '" + getName() + "' has no registered type"
                  + " in the factory");
 }
 

extras/Hadrons/Module.hpp

@@ -35,32 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 BEGIN_HADRONS_NAMESPACE
 
 // module registration macros
-#define MODULE_REGISTER(mod, base)\
-class mod: public base\
-{\
-public:\
-    typedef base Base;\
-    using Base::Base;\
-    virtual std::string getRegisteredName(void)\
-    {\
-        return std::string(#mod);\
-    }\
-};\
-class mod##ModuleRegistrar\
-{\
-public:\
-    mod##ModuleRegistrar(void)\
-    {\
-        ModuleFactory &modFac = ModuleFactory::getInstance();\
-        modFac.registerBuilder(#mod, [&](const std::string name)\
-                              {\
-                                  return std::unique_ptr<mod>(new mod(name));\
-                              });\
-    }\
-};\
-static mod##ModuleRegistrar mod##ModuleRegistrarInstance;
-
-#define MODULE_REGISTER_NS(mod, base, ns)\
+#define MODULE_REGISTER(mod, base, ns)\
 class mod: public base\
 {\
 public:\
@@ -85,6 +60,10 @@ public:\
 };\
 static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance;
 
+#define MODULE_REGISTER_TMP(mod, base, ns)\
+extern template class base;\
+MODULE_REGISTER(mod, ARG(base), ns);
+
 #define ARG(...) __VA_ARGS__
 #define MACRO_REDIRECT(arg1, arg2, arg3, macro, ...) macro
 
@@ -141,10 +120,13 @@ envTmp(type, name, Ls, env().getGrid(Ls))
 MACRO_REDIRECT(__VA_ARGS__, envTmpLat5, envTmpLat4)(__VA_ARGS__)
 
 #define saveResult(ioStem, name, result)\
-if (env().getGrid()->IsBoss())\
+if (env().getGrid()->IsBoss() and !ioStem.empty())\
 {\
-    ResultWriter _writer(RESULT_FILE_NAME(ioStem));\
-    write(_writer, name, result);\
+    makeFileDir(ioStem, env().getGrid());\
+    {\
+        ResultWriter _writer(RESULT_FILE_NAME(ioStem));\
+        write(_writer, name, result);\
+    }\
 }
 
 /******************************************************************************
@@ -172,6 +154,8 @@ public:
     // parse parameters
     virtual void parseParameters(XmlReader &reader, const std::string name) = 0;
     virtual void saveParameters(XmlWriter &writer, const std::string name) = 0;
+    // parameter string
+    virtual std::string parString(void) const = 0;
     // setup
     virtual void setup(void) {};
     virtual void execute(void) = 0;
@@ -200,9 +184,11 @@ public:
     // parse parameters
     virtual void parseParameters(XmlReader &reader, const std::string name);
     virtual void saveParameters(XmlWriter &writer, const std::string name);
+    // parameter string
+    virtual std::string parString(void) const;
     // parameter access
-    const P & par(void) const;
-    void      setPar(const P &par);
+    const P &   par(void) const;
+    void        setPar(const P &par);
 private:
     P par_;
 };
@@ -225,6 +211,8 @@ public:
         push(writer, "options");
         pop(writer);
     };
+    // parameter string (empty)
+    virtual std::string parString(void) const {return "";};
 };
 
 /******************************************************************************
@@ -247,6 +235,16 @@ void Module<P>::saveParameters(XmlWriter &writer, const std::string name)
     write(writer, name, par_);
 }
 
+template <typename P>
+std::string Module<P>::parString(void) const
+{
+    XmlWriter writer("", "");
+
+    write(writer, par_.SerialisableClassName(), par_);
+
+    return writer.string();
+}
+
 template <typename P>
 const P & Module<P>::par(void) const
 {

extras/Hadrons/Modules.hpp

@@ -1,33 +1,3 @@
-/*************************************************************************************
-
-Grid physics library, www.github.com/paboyle/Grid 
-
-Source file: extras/Hadrons/Modules.hpp
-
-Copyright (C) 2015-2018
-
-Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Guido Cossu <guido.cossu@ed.ac.uk>
-Author: Lanny91 <andrew.lawson@gmail.com>
-Author: pretidav <david.preti@csic.es>
-
-This program is free software; you can redistribute it and/or modify
-it under the terms of the GNU General Public License as published by
-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/Meson.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
@@ -64,14 +34,18 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
 #include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
 #include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
 #include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
 #include <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
 #include <Grid/Hadrons/Modules/MIO/LoadNersc.hpp>

extras/Hadrons/Modules/MAction/DWF.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/DWF.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MAction/DWF.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MAction;
+
+template class Grid::Hadrons::MAction::TDWF<FIMPL>;
+

extras/Hadrons/Modules/MAction/DWF.hpp

@@ -61,7 +61,7 @@ public:
     // constructor
     TDWF(const std::string name);
     // destructor
-    virtual ~TDWF(void) = default;
+    virtual ~TDWF(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -72,7 +72,8 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(DWF, TDWF<FIMPL>, MAction);
+extern template class TDWF<FIMPL>;
+MODULE_REGISTER_TMP(DWF, TDWF<FIMPL>, MAction);
 
 /******************************************************************************
  *                        DWF template implementation                         *

extras/Hadrons/Modules/MAction/Wilson.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/Wilson.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MAction;
+
+template class Grid::Hadrons::MAction::TWilson<FIMPL>;
+

extras/Hadrons/Modules/MAction/Wilson.hpp

@@ -59,7 +59,7 @@ public:
     // constructor
     TWilson(const std::string name);
     // destructor
-    virtual ~TWilson(void) = default;
+    virtual ~TWilson(void) {};
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -70,7 +70,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(Wilson, TWilson<FIMPL>, MAction);
+MODULE_REGISTER_TMP(Wilson, TWilson<FIMPL>, MAction);
 
 /******************************************************************************
  *                     TWilson template implementation                        *

extras/Hadrons/Modules/MAction/WilsonClover.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/WilsonClover.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MAction/WilsonClover.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MAction;
+
+template class Grid::Hadrons::MAction::TWilsonClover<FIMPL>;
+

extras/Hadrons/Modules/MAction/WilsonClover.hpp

@@ -64,7 +64,7 @@ public:
     // constructor
     TWilsonClover(const std::string name);
     // destructor
-    virtual ~TWilsonClover(void) = default;
+    virtual ~TWilsonClover(void) {};
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -74,7 +74,7 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(WilsonClover, TWilsonClover<FIMPL>, MAction);
+MODULE_REGISTER_TMP(WilsonClover, TWilsonClover<FIMPL>, MAction);
 
 /******************************************************************************
  *                     TWilsonClover template implementation                        *

extras/Hadrons/Modules/MAction/ZMobiusDWF.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MAction/ZMobiusDWF.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MAction/ZMobiusDWF.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MAction;
+
+template class Grid::Hadrons::MAction::TZMobiusDWF<ZFIMPL>;
+

extras/Hadrons/Modules/MAction/ZMobiusDWF.hpp

@@ -62,7 +62,7 @@ public:
     // constructor
     TZMobiusDWF(const std::string name);
     // destructor
-    virtual ~TZMobiusDWF(void) = default;
+    virtual ~TZMobiusDWF(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -72,7 +72,7 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(ZMobiusDWF, TZMobiusDWF<ZFIMPL>, MAction);
+MODULE_REGISTER_TMP(ZMobiusDWF, TZMobiusDWF<ZFIMPL>, MAction);
 
 /******************************************************************************
  *                 TZMobiusDWF implementation                             *

extras/Hadrons/Modules/MContraction/Baryon.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/Baryon.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TBaryon<FIMPL,FIMPL,FIMPL>;
+

extras/Hadrons/Modules/MContraction/Baryon.hpp

@@ -68,7 +68,7 @@ public:
     // constructor
     TBaryon(const std::string name);
     // destructor
-    virtual ~TBaryon(void) = default;
+    virtual ~TBaryon(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -79,7 +79,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
+MODULE_REGISTER_TMP(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
 
 /******************************************************************************
  *                         TBaryon implementation                             *

extras/Hadrons/Modules/MContraction/DiscLoop.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/DiscLoop.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TDiscLoop<FIMPL>;
+

extras/Hadrons/Modules/MContraction/DiscLoop.hpp

@@ -65,7 +65,7 @@ public:
     // constructor
     TDiscLoop(const std::string name);
     // destructor
-    virtual ~TDiscLoop(void) = default;
+    virtual ~TDiscLoop(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -76,7 +76,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(DiscLoop, TDiscLoop<FIMPL>, MContraction);
+MODULE_REGISTER_TMP(DiscLoop, TDiscLoop<FIMPL>, MContraction);
 
 /******************************************************************************
  *                       TDiscLoop implementation                             *

extras/Hadrons/Modules/MContraction/Gamma3pt.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/Gamma3pt.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TGamma3pt<FIMPL,FIMPL,FIMPL>;
+

extras/Hadrons/Modules/MContraction/Gamma3pt.hpp

@@ -96,7 +96,7 @@ public:
     // constructor
     TGamma3pt(const std::string name);
     // destructor
-    virtual ~TGamma3pt(void) = default;
+    virtual ~TGamma3pt(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -107,7 +107,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(Gamma3pt, ARG(TGamma3pt<FIMPL, FIMPL, FIMPL>), MContraction);
+MODULE_REGISTER_TMP(Gamma3pt, ARG(TGamma3pt<FIMPL, FIMPL, FIMPL>), MContraction);
 
 /******************************************************************************
  *                       TGamma3pt implementation                             *

extras/Hadrons/Modules/MContraction/Meson.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/Meson.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TMeson<FIMPL,FIMPL>;
+

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -90,7 +90,7 @@ public:
     // constructor
     TMeson(const std::string name);
     // destructor
-    virtual ~TMeson(void) = default;
+    virtual ~TMeson(void) {};
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -102,7 +102,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(Meson, ARG(TMeson<FIMPL, FIMPL>), MContraction);
+MODULE_REGISTER_TMP(Meson, ARG(TMeson<FIMPL, FIMPL>), MContraction);
 
 /******************************************************************************
  *                           TMeson implementation                            *

extras/Hadrons/Modules/MContraction/WardIdentity.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WardIdentity.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TWardIdentity<FIMPL>;
+

extras/Hadrons/Modules/MContraction/WardIdentity.hpp

@@ -71,7 +71,7 @@ public:
     // constructor
     TWardIdentity(const std::string name);
     // destructor
-    virtual ~TWardIdentity(void) = default;
+    virtual ~TWardIdentity(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -84,7 +84,7 @@ private:
     unsigned int Ls_;
 };
 
-MODULE_REGISTER_NS(WardIdentity, TWardIdentity<FIMPL>, MContraction);
+MODULE_REGISTER_TMP(WardIdentity, TWardIdentity<FIMPL>, MContraction);
 
 /******************************************************************************
  *                     TWardIdentity implementation                           *
@@ -119,7 +119,7 @@ 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");
+        HADRONS_ERROR(Size, "Ls mismatch between quark action and propagator");
     }
     envTmpLat(PropagatorField, "tmp");
     envTmpLat(PropagatorField, "vector_WI");

extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp

@@ -97,7 +97,7 @@ public:\
     /* constructor */ \
     T##modname(const std::string name);\
     /* destructor */ \
-    virtual ~T##modname(void) = default;\
+    virtual ~T##modname(void) {};\
     /* dependency relation */ \
     virtual std::vector<std::string> getInput(void);\
     virtual std::vector<std::string> getOutput(void);\
@@ -109,7 +109,7 @@ protected:\
     /* execution */ \
     virtual void execute(void);\
 };\
-MODULE_REGISTER_NS(modname, T##modname, MContraction);
+MODULE_REGISTER(modname, T##modname, MContraction);
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MFermion/GaugeProp.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MFermion/GaugeProp.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MFermion;
+
+template class Grid::Hadrons::MFermion::TGaugeProp<FIMPL>;
+

extras/Hadrons/Modules/MFermion/GaugeProp.hpp

@@ -81,7 +81,7 @@ public:
     // constructor
     TGaugeProp(const std::string name);
     // destructor
-    virtual ~TGaugeProp(void) = default;
+    virtual ~TGaugeProp(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -95,7 +95,7 @@ private:
     SolverFn     *solver_{nullptr};
 };
 
-MODULE_REGISTER_NS(GaugeProp, TGaugeProp<FIMPL>, MFermion);
+MODULE_REGISTER_TMP(GaugeProp, TGaugeProp<FIMPL>, MFermion);
 /******************************************************************************
  *                      TGaugeProp implementation                             *
  ******************************************************************************/
@@ -177,7 +177,7 @@ void TGaugeProp<FImpl>::execute(void)
         {
             if (Ls_ != env().getObjectLs(par().source))
             {
-                HADRON_ERROR(Size, "Ls mismatch between quark action and source");
+                HADRONS_ERROR(Size, "Ls mismatch between quark action and source");
             }
             else
             {

extras/Hadrons/Modules/MGauge/FundtoHirep.cc

@@ -4,9 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid
 
 Source file: extras/Hadrons/Modules/MGauge/FundtoHirep.cc
 
-Copyright (C) 2015
-Copyright (C) 2016
+Copyright (C) 2015-2018
 
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
+Author: pretidav <david.preti@csic.es>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by

extras/Hadrons/Modules/MGauge/FundtoHirep.hpp

@@ -55,7 +55,7 @@ public:
     // constructor
     TFundtoHirep(const std::string name);
     // destructor
-    virtual ~TFundtoHirep(void) = default;
+    virtual ~TFundtoHirep(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -65,9 +65,9 @@ public:
     void execute(void);
 };
 
-//MODULE_REGISTER_NS(FundtoAdjoint,   TFundtoHirep<AdjointRepresentation>, MGauge);
-//MODULE_REGISTER_NS(FundtoTwoIndexSym, TFundtoHirep<TwoIndexSymmetricRepresentation>, MGauge);
-//MODULE_REGISTER_NS(FundtoTwoIndexAsym, TFundtoHirep<TwoIndexAntiSymmetricRepresentation>, MGauge);
+//MODULE_REGISTER_TMP(FundtoAdjoint,   TFundtoHirep<AdjointRepresentation>, MGauge);
+//MODULE_REGISTER_TMP(FundtoTwoIndexSym, TFundtoHirep<TwoIndexSymmetricRepresentation>, MGauge);
+//MODULE_REGISTER_TMP(FundtoTwoIndexAsym, TFundtoHirep<TwoIndexAntiSymmetricRepresentation>, MGauge);
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MGauge/Random.hpp

@@ -46,7 +46,7 @@ public:
     // constructor
     TRandom(const std::string name);
     // destructor
-    virtual ~TRandom(void) = default;
+    virtual ~TRandom(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -57,7 +57,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(Random, TRandom, MGauge);
+MODULE_REGISTER(Random, TRandom, MGauge);
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MGauge/StochEm.cc

@@ -7,6 +7,8 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.cc
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: James Harrison <j.harrison@soton.ac.uk>
+Author: Vera Guelpers <vmg1n14@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by

extras/Hadrons/Modules/MGauge/StochEm.hpp

@@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: James Harrison <j.harrison@soton.ac.uk>
 Author: Vera Guelpers <vmg1n14@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
@@ -59,7 +60,7 @@ public:
     // constructor
     TStochEm(const std::string name);
     // destructor
-    virtual ~TStochEm(void) = default;
+    virtual ~TStochEm(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -72,7 +73,7 @@ private:
     bool    weightDone_;
 };
 
-MODULE_REGISTER_NS(StochEm, TStochEm, MGauge);
+MODULE_REGISTER(StochEm, TStochEm, MGauge);
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MGauge/Unit.hpp

@@ -46,7 +46,7 @@ public:
     // constructor
     TUnit(const std::string name);
     // destructor
-    virtual ~TUnit(void) = default;
+    virtual ~TUnit(void) {};
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -57,7 +57,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(Unit, TUnit, MGauge);
+MODULE_REGISTER(Unit, TUnit, MGauge);
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MGauge/UnitEm.cc

@@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.cc
 Copyright (C) 2015
 Copyright (C) 2016
 
+Author: James Harrison <j.harrison@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by

extras/Hadrons/Modules/MGauge/UnitEm.hpp

@@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MGauge/StochEm.hpp
 Copyright (C) 2015
 Copyright (C) 2016
 
+Author: James Harrison <j.harrison@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by

extras/Hadrons/Modules/MIO/LoadBinary.cc

@@ -0,0 +1,40 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MIO/LoadBinary.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MIO/LoadBinary.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MIO;
+
+template class Grid::Hadrons::MIO::TLoadBinary<GIMPL>;
+template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MIO::TLoadBinary<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MIO/LoadBinary.hpp

@@ -61,7 +61,7 @@ public:
     // constructor
     TLoadBinary(const std::string name);
     // destructor
-    virtual ~TLoadBinary(void) = default;
+    virtual ~TLoadBinary(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -71,12 +71,12 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(LoadBinary, TLoadBinary<GIMPL>, MIO);
-MODULE_REGISTER_NS(LoadBinaryScalarSU2, TLoadBinary<ScalarNxNAdjImplR<2>>, MIO);
-MODULE_REGISTER_NS(LoadBinaryScalarSU3, TLoadBinary<ScalarNxNAdjImplR<3>>, MIO);
-MODULE_REGISTER_NS(LoadBinaryScalarSU4, TLoadBinary<ScalarNxNAdjImplR<4>>, MIO);
-MODULE_REGISTER_NS(LoadBinaryScalarSU5, TLoadBinary<ScalarNxNAdjImplR<5>>, MIO);
-MODULE_REGISTER_NS(LoadBinaryScalarSU6, TLoadBinary<ScalarNxNAdjImplR<6>>, MIO);
+MODULE_REGISTER_TMP(LoadBinary, TLoadBinary<GIMPL>, MIO);
+MODULE_REGISTER_TMP(LoadBinaryScalarSU2, TLoadBinary<ScalarNxNAdjImplR<2>>, MIO);
+MODULE_REGISTER_TMP(LoadBinaryScalarSU3, TLoadBinary<ScalarNxNAdjImplR<3>>, MIO);
+MODULE_REGISTER_TMP(LoadBinaryScalarSU4, TLoadBinary<ScalarNxNAdjImplR<4>>, MIO);
+MODULE_REGISTER_TMP(LoadBinaryScalarSU5, TLoadBinary<ScalarNxNAdjImplR<5>>, MIO);
+MODULE_REGISTER_TMP(LoadBinaryScalarSU6, TLoadBinary<ScalarNxNAdjImplR<6>>, MIO);
 
 /******************************************************************************
  *                         TLoadBinary implementation                         *

extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MIO;
+
+template class Grid::Hadrons::MIO::TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>>;
+

extras/Hadrons/Modules/MIO/LoadCoarseEigenPack.hpp

@@ -45,6 +45,7 @@ class LoadCoarseEigenPackPar: Serializable
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(LoadCoarseEigenPackPar,
                                     std::string, filestem,
+                                    bool,         multiFile,
                                     unsigned int, sizeFine,
                                     unsigned int, sizeCoarse,
                                     unsigned int, Ls,
@@ -56,11 +57,14 @@ class TLoadCoarseEigenPack: public Module<LoadCoarseEigenPackPar>
 {
 public:
     typedef CoarseEigenPack<typename Pack::Field, typename Pack::CoarseField> BasePack;
+    template <typename vtype> 
+    using iImplScalar = iScalar<iScalar<iScalar<vtype>>>;
+    typedef iImplScalar<typename Pack::Field::vector_type> SiteComplex;
 public:
     // constructor
     TLoadCoarseEigenPack(const std::string name);
     // destructor
-    virtual ~TLoadCoarseEigenPack(void) = default;
+    virtual ~TLoadCoarseEigenPack(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -70,8 +74,7 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(LoadCoarseFermionEigenPack, 
-    ARG(TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>>), MIO);
+MODULE_REGISTER_TMP(LoadCoarseFermionEigenPack, ARG(TLoadCoarseEigenPack<CoarseFermionEigenPack<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>>), MIO);
 
 /******************************************************************************
  *                 TLoadCoarseEigenPack implementation                             *
@@ -114,9 +117,15 @@ void TLoadCoarseEigenPack<Pack>::setup(void)
 template <typename Pack>
 void TLoadCoarseEigenPack<Pack>::execute(void)
 {
-    auto &epack = envGetDerived(BasePack, Pack, getName());
+    auto                 cg     = env().getCoarseGrid(par().blockSize, par().Ls);
+    auto                 &epack = envGetDerived(BasePack, Pack, getName());
+    Lattice<SiteComplex> dummy(cg);
 
-    epack.read(par().filestem, vm().getTrajectory());
+    epack.read(par().filestem, par().multiFile, vm().getTrajectory());
+    LOG(Message) << "Block Gramm-Schmidt pass 1"<< std::endl;
+    blockOrthogonalise(dummy, epack.evec);
+    LOG(Message) << "Block Gramm-Schmidt pass 2"<< std::endl;
+    blockOrthogonalise(dummy, epack.evec);
 }
 
 END_MODULE_NAMESPACE

extras/Hadrons/Modules/MIO/LoadEigenPack.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MIO/LoadEigenPack.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MIO/LoadEigenPack.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MIO;
+
+template class Grid::Hadrons::MIO::TLoadEigenPack<FermionEigenPack<FIMPL>>;
+

extras/Hadrons/Modules/MIO/LoadEigenPack.hpp

@@ -45,6 +45,7 @@ class LoadEigenPackPar: Serializable
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(LoadEigenPackPar,
                                     std::string, filestem,
+                                    bool, multiFile,
                                     unsigned int, size,
                                     unsigned int, Ls);
 };
@@ -58,7 +59,7 @@ public:
     // constructor
     TLoadEigenPack(const std::string name);
     // destructor
-    virtual ~TLoadEigenPack(void) = default;
+    virtual ~TLoadEigenPack(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -68,7 +69,7 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(LoadFermionEigenPack, TLoadEigenPack<FermionEigenPack<FIMPL>>, MIO);
+MODULE_REGISTER_TMP(LoadFermionEigenPack, TLoadEigenPack<FermionEigenPack<FIMPL>>, MIO);
 
 /******************************************************************************
  *                    TLoadEigenPack implementation                           *
@@ -111,7 +112,8 @@ void TLoadEigenPack<Pack>::execute(void)
 {
     auto &epack = envGetDerived(BasePack, Pack, getName());
 
-    epack.read(par().filestem, vm().getTrajectory());
+    epack.read(par().filestem, par().multiFile, vm().getTrajectory());
+    epack.eval.resize(par().size);
 }
 
 END_MODULE_NAMESPACE

extras/Hadrons/Modules/MIO/LoadNersc.hpp

@@ -52,7 +52,7 @@ public:
     // constructor
     TLoadNersc(const std::string name);
     // destructor
-    virtual ~TLoadNersc(void) = default;
+    virtual ~TLoadNersc(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -62,7 +62,7 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(LoadNersc, TLoadNersc, MIO);
+MODULE_REGISTER(LoadNersc, TLoadNersc, MIO);
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MLoop/NoiseLoop.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MLoop/NoiseLoop.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MLoop;
+
+template class Grid::Hadrons::MLoop::TNoiseLoop<FIMPL>;
+

extras/Hadrons/Modules/MLoop/NoiseLoop.hpp

@@ -71,7 +71,7 @@ public:
     // constructor
     TNoiseLoop(const std::string name);
     // destructor
-    virtual ~TNoiseLoop(void) = default;
+    virtual ~TNoiseLoop(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -82,7 +82,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(NoiseLoop, TNoiseLoop<FIMPL>, MLoop);
+MODULE_REGISTER_TMP(NoiseLoop, TNoiseLoop<FIMPL>, MLoop);
 
 /******************************************************************************
  *                 TNoiseLoop implementation                                  *

extras/Hadrons/Modules/MScalar/ChargedProp.hpp

@@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MScalar/ChargedProp.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: James Harrison <j.harrison@soton.ac.uk>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -81,7 +82,7 @@ public:
     // constructor
     TChargedProp(const std::string name);
     // destructor
-    virtual ~TChargedProp(void) = default;
+    virtual ~TChargedProp(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -103,7 +104,7 @@ private:
     std::vector<ScalarField *> phase_;
 };
 
-MODULE_REGISTER_NS(ChargedProp, TChargedProp, MScalar);
+MODULE_REGISTER(ChargedProp, TChargedProp, MScalar);
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MScalar/FreeProp.hpp

@@ -56,7 +56,7 @@ public:
     // constructor
     TFreeProp(const std::string name);
     // destructor
-    virtual ~TFreeProp(void) = default;
+    virtual ~TFreeProp(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -70,7 +70,7 @@ private:
     bool        freePropDone_;
 };
 
-MODULE_REGISTER_NS(FreeProp, TFreeProp, MScalar);
+MODULE_REGISTER(FreeProp, TFreeProp, MScalar);
 
 END_MODULE_NAMESPACE
 

extras/Hadrons/Modules/MScalarSUN/Div.cc

@@ -0,0 +1,39 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/Div.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/Div.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TDiv<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MScalarSUN/Div.hpp

@@ -49,24 +49,25 @@ public:
                                     std::string,              output);
 };
 
+class DivResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(DivResult,
+                                    DiffType, type,
+                                    Complex,  value);
+};
+
 template <typename SImpl>
 class TDiv: public Module<DivPar>
 {
 public:
     typedef typename SImpl::Field        Field;
     typedef typename SImpl::ComplexField ComplexField;
-    class Result: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        DiffType, type,
-                                        Complex,  value);
-    };
 public:
     // constructor
     TDiv(const std::string name);
     // destructor
-    virtual ~TDiv(void) = default;
+    virtual ~TDiv(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -76,11 +77,11 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(DivSU2, TDiv<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_NS(DivSU3, TDiv<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_NS(DivSU4, TDiv<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_NS(DivSU5, TDiv<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_NS(DivSU6, TDiv<ScalarNxNAdjImplR<6>>, MScalarSUN);
+MODULE_REGISTER_TMP(DivSU2, TDiv<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(DivSU3, TDiv<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(DivSU4, TDiv<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(DivSU5, TDiv<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(DivSU6, TDiv<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
  *                           TDiv implementation                              *
@@ -112,7 +113,7 @@ void TDiv<SImpl>::setup(void)
 {
     if (par().op.size() != env().getNd())
     {
-        HADRON_ERROR(Size, "the number of components differs from number of dimensions");
+        HADRONS_ERROR(Size, "the number of components differs from number of dimensions");
     }
     envCreateLat(ComplexField, getName());
 }
@@ -126,7 +127,7 @@ void TDiv<SImpl>::execute(void)
     LOG(Message) << "Computing the " << par().type << " divergence of [";
     for (unsigned int mu = 0; mu < nd; ++mu)
     {
-        std::cout << par().op[mu] << ((mu == nd - 1) ? "]" : ", ");
+        std::cout << "'" << par().op[mu] << ((mu == nd - 1) ? "']" : "', ");
     }
     std::cout << std::endl;
 
@@ -139,7 +140,7 @@ void TDiv<SImpl>::execute(void)
     }
     if (!par().output.empty())
     {
-        Result       r;
+        DivResult r;
 
         r.type  = par().type;
         r.value = TensorRemove(sum(div));

extras/Hadrons/Modules/MScalarSUN/EMT.cc

@@ -0,0 +1,39 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/EMT.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/EMT.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TEMT<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MScalarSUN/EMT.hpp

@@ -54,6 +54,17 @@ public:
                                     std::string, output);
 };
 
+class EMTResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(EMTResult,
+                                    std::vector<std::vector<Complex>>, value,
+                                    double,                            m2,
+                                    double,                            lambda,
+                                    double,                            g,
+                                    double,                            xi);
+};
+
 template <typename SImpl>
 class TEMT: public Module<EMTPar>
 {
@@ -64,7 +75,7 @@ public:
     // constructor
     TEMT(const std::string name);
     // destructor
-    virtual ~TEMT(void) = default;
+    virtual ~TEMT(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -74,11 +85,11 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(EMTSU2, TEMT<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_NS(EMTSU3, TEMT<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_NS(EMTSU4, TEMT<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_NS(EMTSU5, TEMT<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_NS(EMTSU6, TEMT<ScalarNxNAdjImplR<6>>, MScalarSUN);
+MODULE_REGISTER_TMP(EMTSU2, TEMT<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(EMTSU3, TEMT<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(EMTSU4, TEMT<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(EMTSU5, TEMT<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(EMTSU6, TEMT<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
  *                           TEMT implementation                              *
@@ -99,8 +110,12 @@ std::vector<std::string> TEMT<SImpl>::getInput(void)
     for (unsigned int nu = mu; nu < env().getNd(); ++nu)
     {
         in.push_back(varName(par().kinetic, mu, nu));
-        in.push_back(varName(par().improvement, mu, nu));
+        if (!par().improvement.empty())
+        {
+            in.push_back(varName(par().improvement, mu, nu));
+        }
     }
+    in.push_back(varName(par().kinetic, "sum"));
     in.push_back(varName(par().phiPow, 2));
     in.push_back(varName(par().phiPow, 4));
 
@@ -130,7 +145,6 @@ void TEMT<SImpl>::setup(void)
     {
         envCreateLat(ComplexField, varName(getName(), mu, nu));
     }
-    envTmpLat(ComplexField, "sumkin");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
@@ -140,37 +154,59 @@ void TEMT<SImpl>::execute(void)
     LOG(Message) << "Computing energy-momentum tensor" << std::endl;
     LOG(Message) << "  kinetic terms: '" << par().kinetic << "'" << std::endl;
     LOG(Message) << "      tr(phi^n): '" << par().phiPow << "'" << std::endl;
-    LOG(Message) << "    improvement: '" << par().improvement << "'" << std::endl;
+    if (!par().improvement.empty())
+    {
+        LOG(Message) << "    improvement: '" << par().improvement << "'" << std::endl;
+    }
     LOG(Message) << "            m^2= " << par().m2 << std::endl;
     LOG(Message) << "         lambda= " << par().lambda << std::endl;
     LOG(Message) << "              g= " << par().g << std::endl;
-    LOG(Message) << "             xi= " << par().xi << std::endl;
+    if (!par().improvement.empty())
+    {
+        LOG(Message) << "             xi= " << par().xi << std::endl;
+    }
 
-    const unsigned int N       = SImpl::Group::Dimension;
+    const unsigned int N = SImpl::Group::Dimension, nd = env().getNd();
     auto               &trphi2 = envGet(ComplexField, varName(par().phiPow, 2));
     auto               &trphi4 = envGet(ComplexField, varName(par().phiPow, 4));
+    auto               &sumkin = envGet(ComplexField, varName(par().kinetic, "sum"));
+    EMTResult          result;
 
-    envGetTmp(ComplexField, sumkin);
-    sumkin = zero;
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    if (!par().output.empty())
     {
-        auto &trkin = envGet(ComplexField, varName(par().kinetic, mu, mu));
-
-        sumkin += trkin;
+        result.m2     = par().m2;
+        result.g      = par().g;
+        result.lambda = par().lambda;
+        result.xi     = par().xi;
+        result.value.resize(nd, std::vector<Complex>(nd));
     }
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    for (unsigned int mu = 0; mu < nd; ++mu)
+    for (unsigned int nu = mu; nu < nd; ++nu)
     {
         auto &out   = envGet(ComplexField, varName(getName(), mu, nu));
         auto &trkin = envGet(ComplexField, varName(par().kinetic, mu, nu));
-        auto &imp   = envGet(ComplexField, varName(par().improvement, mu, nu));
+        
+        out = 2.*trkin;
+        if (!par().improvement.empty())
+        {
+            auto &imp = envGet(ComplexField, varName(par().improvement, mu, nu));
 
-        out = 2.*trkin + par().xi*imp;
+            out += par().xi*imp;
+        }
         if (mu == nu)
         {
             out -= sumkin + par().m2*trphi2 + par().lambda*trphi4;
         }
         out *= N/par().g;
+        if (!par().output.empty())
+        {
+            result.value[mu][nu] = TensorRemove(sum(out));
+            result.value[mu][nu] = result.value[nu][mu];
+        }
+    }
+    if (!par().output.empty())
+    {
+        saveResult(par().output, "emt", result);
     }
 }
 

extras/Hadrons/Modules/MScalarSUN/Grad.cc

@@ -0,0 +1,38 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/Grad.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/Grad.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TGrad<ScalarNxNAdjImplR<6>>;

extras/Hadrons/Modules/MScalarSUN/Grad.hpp

@@ -0,0 +1,166 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/Grad.hpp
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef Hadrons_MScalarSUN_Grad_hpp_
+#define Hadrons_MScalarSUN_Grad_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                       Gradient of a complex field                          *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+class GradPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(GradPar,
+                                    std::string, op,
+                                    DiffType,    type,
+                                    std::string, output);
+};
+
+class GradResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(GradResult,
+                                    DiffType,              type,
+                                    std::vector<Complex>,  value);
+};
+
+template <typename SImpl>
+class TGrad: public Module<GradPar>
+{
+public:
+    typedef typename SImpl::Field        Field;
+    typedef typename SImpl::ComplexField ComplexField;
+public:
+    // constructor
+    TGrad(const std::string name);
+    // destructor
+    virtual ~TGrad(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(GradSU2, TGrad<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(GradSU3, TGrad<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(GradSU4, TGrad<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(GradSU5, TGrad<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(GradSU6, TGrad<ScalarNxNAdjImplR<6>>, MScalarSUN);
+
+/******************************************************************************
+ *                         TGrad implementation                               *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename SImpl>
+TGrad<SImpl>::TGrad(const std::string name)
+: Module<GradPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename SImpl>
+std::vector<std::string> TGrad<SImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().op};
+    
+    return in;
+}
+
+template <typename SImpl>
+std::vector<std::string> TGrad<SImpl>::getOutput(void)
+{
+    std::vector<std::string> out;
+    const auto               nd = env().getNd();
+
+    for (unsigned int mu = 0; mu < nd; ++mu)
+    {
+        out.push_back(varName(getName(), mu));
+    }
+
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TGrad<SImpl>::setup(void)
+{
+    const auto nd = env().getNd();
+
+    for (unsigned int mu = 0; mu < nd; ++mu)
+    {
+        envCreateLat(ComplexField, varName(getName(), mu));
+    }
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TGrad<SImpl>::execute(void)
+{
+    LOG(Message) << "Computing the " << par().type << " gradient of '"
+                 << par().op << "'" << std::endl;
+
+    const unsigned int nd = env().getNd();
+    GradResult         result;
+    auto               &op = envGet(ComplexField, par().op);
+
+    if (!par().output.empty())
+    {
+        result.type = par().type;
+        result.value.resize(nd);
+    }
+    for (unsigned int mu = 0; mu < nd; ++mu)
+    {
+        auto &der = envGet(ComplexField, varName(getName(), mu));
+
+        dmu(der, op, mu, par().type);
+        if (!par().output.empty())
+        {
+            result.value[mu] = TensorRemove(sum(der));
+        }
+    }
+    if (!par().output.empty())
+    {
+        saveResult(par().output, "grad", result);
+    }
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_Grad_hpp_

extras/Hadrons/Modules/MScalarSUN/ShiftProbe.cc

@@ -0,0 +1,39 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/ShiftProbe.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TShiftProbe<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MScalarSUN/ShiftProbe.hpp

@@ -40,7 +40,7 @@ BEGIN_HADRONS_NAMESPACE
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MScalarSUN)
 
-typedef std::pair<unsigned int, unsigned int> ShiftPair;
+typedef std::pair<int, int> ShiftPair;
 
 class ShiftProbePar: Serializable
 {
@@ -51,25 +51,25 @@ public:
                                     std::string, output);
 };
 
+class ShiftProbeResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ShiftProbeResult,
+                                    std::string, shifts,
+                                    Complex,     value);
+};
+
 template <typename SImpl>
 class TShiftProbe: public Module<ShiftProbePar>
 {
 public:
-    
     typedef typename SImpl::Field                          Field;
     typedef typename SImpl::ComplexField                   ComplexField;
-    class Result: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        std::string, op,
-                                        Complex    , value);
-    };
 public:
     // constructor
     TShiftProbe(const std::string name);
     // destructor
-    virtual ~TShiftProbe(void) = default;
+    virtual ~TShiftProbe(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -79,11 +79,11 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(ShiftProbeSU2, TShiftProbe<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_NS(ShiftProbeSU3, TShiftProbe<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_NS(ShiftProbeSU4, TShiftProbe<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_NS(ShiftProbeSU5, TShiftProbe<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_NS(ShiftProbeSU6, TShiftProbe<ScalarNxNAdjImplR<6>>, MScalarSUN);
+MODULE_REGISTER_TMP(ShiftProbeSU2, TShiftProbe<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(ShiftProbeSU3, TShiftProbe<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(ShiftProbeSU4, TShiftProbe<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(ShiftProbeSU5, TShiftProbe<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(ShiftProbeSU6, TShiftProbe<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
  *                        TShiftProbe implementation                          *
@@ -127,27 +127,27 @@ void TShiftProbe<SImpl>::execute(void)
                  << std::endl;
 
     std::vector<ShiftPair> shift;
-    unsigned int           sign;
+    double                 sign;
     auto                   &phi   = envGet(Field, par().field);
     auto                   &probe = envGet(ComplexField, getName());
 
     shift = strToVec<ShiftPair>(par().shifts);
     if (shift.size() % 2 != 0)
     {
-        HADRON_ERROR(Size, "the number of shifts is odd");
+        HADRONS_ERROR(Size, "the number of shifts is odd");
     }
-    sign = (shift.size() % 4 == 0) ? 1 : -1;
+    sign = (shift.size() % 4 == 0) ? 1. : -1.;
     for (auto &s: shift)
     {
         if (s.first >= env().getNd())
         {
-            HADRON_ERROR(Size, "dimension to large for shift <" 
+            HADRONS_ERROR(Size, "dimension to large for shift <" 
                                + std::to_string(s.first) + " " 
                                + std::to_string(s.second) + ">" );
         }
     }
     envGetTmp(Field, acc);
-    acc   = 1.;
+    acc = 1.;
     for (unsigned int i = 0; i < shift.size(); ++i)
     {
         if (shift[i].second == 0)
@@ -160,6 +160,14 @@ void TShiftProbe<SImpl>::execute(void)
         }
     }
     probe = sign*trace(acc);
+    if (!par().output.empty())
+    {
+        ShiftProbeResult r;
+
+        r.shifts = par().shifts;
+        r.value  = TensorRemove(sum(probe));
+        saveResult(par().output, "probe", r);
+    }
 }
 
 END_MODULE_NAMESPACE

extras/Hadrons/Modules/MScalarSUN/StochFreeField.cc

@@ -0,0 +1,11 @@
+#include <Grid/Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TStochFreeField<ScalarNxNAdjImplR<6>>;

extras/Hadrons/Modules/MScalarSUN/StochFreeField.hpp

@@ -0,0 +1,151 @@
+#ifndef Hadrons_MScalarSUN_StochFreeField_hpp_
+#define Hadrons_MScalarSUN_StochFreeField_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                      stochastic free SU(N) scalar field                    *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+class StochFreeFieldPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(StochFreeFieldPar,
+                                    double, m2,
+                                    double, g);
+};
+
+template <typename SImpl>
+class TStochFreeField: public Module<StochFreeFieldPar>
+{
+public:
+    typedef typename SImpl::Field                    Field;
+    typedef typename SImpl::ComplexField             ComplexField;
+    typedef typename SImpl::Group                    Group;
+    typedef typename SImpl::SiteField::scalar_object Site;
+public:
+    // constructor
+    TStochFreeField(const std::string name);
+    // destructor
+    virtual ~TStochFreeField(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+private:
+    bool create_weight;
+};
+
+MODULE_REGISTER_TMP(StochFreeFieldSU2, TStochFreeField<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(StochFreeFieldSU3, TStochFreeField<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(StochFreeFieldSU4, TStochFreeField<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(StochFreeFieldSU5, TStochFreeField<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(StochFreeFieldSU6, TStochFreeField<ScalarNxNAdjImplR<6>>, MScalarSUN);
+
+/******************************************************************************
+ *                 TStochFreeField implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename SImpl>
+TStochFreeField<SImpl>::TStochFreeField(const std::string name)
+: Module<StochFreeFieldPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename SImpl>
+std::vector<std::string> TStochFreeField<SImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename SImpl>
+std::vector<std::string> TStochFreeField<SImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TStochFreeField<SImpl>::setup(void)
+{
+    create_weight = false; 
+    if (!env().hasCreatedObject("_" + getName() + "_weight"))
+    {
+        envCacheLat(ComplexField, "_" + getName() + "_weight");
+        create_weight = true;
+    }
+    envTmpLat(Field, "phift");
+    envTmpLat(ComplexField, "ca");
+    envCreateLat(Field, getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TStochFreeField<SImpl>::execute(void)
+{
+    LOG(Message) << "Generating stochastic scalar field" << std::endl;
+    
+    const   unsigned int N    = Group::Dimension;
+    const   unsigned int Nadj = Group::AdjointDimension;
+    auto    &phi              = envGet(Field, getName());
+    auto    &w                = envGet(ComplexField, "_" + getName() + "_weight");
+    auto    &rng              = *env().get4dRng();
+    double  trphi2;
+    FFT     fft(env().getGrid());
+    Integer vol;
+
+    vol = 1;
+    for(int d = 0; d < env().getNd(); d++)
+    {
+        vol = vol*env().getDim(d);
+    }
+    if (create_weight)
+    {
+        LOG(Message) << "Caching momentum-space scalar action" << std::endl;
+        
+        SImpl::MomentumSpacePropagator(w, sqrt(par().m2));
+        w *= par().g/N;
+        w  = sqrt(vol)*sqrt(w);
+    }
+    LOG(Message) << "Generating random momentum-space field" << std::endl;
+    envGetTmp(Field, phift);
+    envGetTmp(ComplexField, ca);
+    phift = zero;
+    for (int a = 0; a < Nadj; ++a) 
+    {
+        Site ta;
+
+        gaussian(rng, ca);
+        Group::generator(a, ta);
+        phift += ca*ta;
+    }
+    phift *= w;
+    LOG(Message) << "Field Fourier transform" << std::endl;
+    fft.FFT_all_dim(phi, phift, FFT::backward);
+    phi = 0.5*(phi - adj(phi));
+    trphi2 = -TensorRemove(sum(trace(phi*phi))).real()/vol;
+    LOG(Message) << "tr(phi^2)= " << trphi2 << std::endl;
+
+    // ComplexField phi2(env().getGrid());
+
+    // phi2=trace(phi*phi);
+    // std::cout << phi2 << std::endl;
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_StochFreeField_hpp_

extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.cc

@@ -0,0 +1,11 @@
+#include <Grid/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TTimeMomProbe<ScalarNxNAdjImplR<6>>;

extras/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp

@@ -0,0 +1,241 @@
+#ifndef Hadrons_MScalarSUN_TimeMomProbe_hpp_
+#define Hadrons_MScalarSUN_TimeMomProbe_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *          n-point functions O(t,p)*tr(phi(t_1,p_1)*...*phi(t_n,p_n))        *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+class TimeMomProbePar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbePar,
+                                    std::string,              field,
+                                    std::vector<std::string>, op,
+                                    std::vector<std::vector<std::string>>, timeMom,
+                                    std::string,              output);
+};
+
+class TimeMomProbeResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TimeMomProbeResult,
+                                    std::string,                   op,
+                                    std::vector<std::vector<int>>, timeMom,
+                                    std::vector<Complex>,          data);
+};
+
+template <typename SImpl>
+class TTimeMomProbe: public Module<TimeMomProbePar>
+{
+public:
+    typedef typename SImpl::Field                    Field;
+    typedef typename SImpl::SiteField::scalar_object Site;
+    typedef typename SImpl::ComplexField             ComplexField;
+    typedef          std::vector<Complex>            SlicedOp;
+public:
+    // constructor
+    TTimeMomProbe(const std::string name);
+    // destructor
+    virtual ~TTimeMomProbe(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+private:
+    void vectorModulo(std::vector<int> &v);
+};
+
+MODULE_REGISTER_TMP(TimeMomProbeSU2, TTimeMomProbe<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(TimeMomProbeSU3, TTimeMomProbe<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(TimeMomProbeSU4, TTimeMomProbe<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(TimeMomProbeSU5, TTimeMomProbe<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(TimeMomProbeSU6, TTimeMomProbe<ScalarNxNAdjImplR<6>>, MScalarSUN);
+
+/******************************************************************************
+ *                        TTimeMomProbe implementation                        *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename SImpl>
+TTimeMomProbe<SImpl>::TTimeMomProbe(const std::string name)
+: Module<TimeMomProbePar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename SImpl>
+std::vector<std::string> TTimeMomProbe<SImpl>::getInput(void)
+{
+    std::vector<std::string> in = par().op;
+    
+    in.push_back(par().field);
+
+    return in;
+}
+
+template <typename SImpl>
+std::vector<std::string> TTimeMomProbe<SImpl>::getOutput(void)
+{
+    std::vector<std::string> out;
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TTimeMomProbe<SImpl>::setup(void)
+{
+    envTmpLat(ComplexField, "ftBuf");
+    envTmpLat(Field, "ftMatBuf");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+// NB: time is direction 0
+template <typename SImpl>
+void TTimeMomProbe<SImpl>::vectorModulo(std::vector<int> &v)
+{
+    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    {
+        auto d = env().getDim(mu);
+        v[mu] = ((v[mu] % d) + d) % d;
+    }
+}
+
+template <typename SImpl>
+void TTimeMomProbe<SImpl>::execute(void)
+{
+    const unsigned int                           nd = env().getNd();
+    const unsigned int                           nt = env().getDim(0);
+    double                                       partVol = 1.;
+    std::set<std::vector<int>>                   timeMomSet;
+    std::vector<std::vector<std::vector<int>>>   timeMom;
+    std::vector<std::vector<int>>                transferMom;
+    FFT                                          fft(env().getGrid());
+    std::vector<int>                             dMask(nd, 1);
+    std::vector<TimeMomProbeResult>              result;
+    std::map<std::string, std::vector<SlicedOp>> slicedOp;
+    std::vector<SlicedOp>                        slicedProbe;
+    auto                                         &phi = envGet(Field, par().field);
+
+    envGetTmp(ComplexField, ftBuf);
+    envGetTmp(Field, ftMatBuf);
+    dMask[0] = 0;
+    for (unsigned int mu = 1; mu < nd; ++mu)
+    {
+        partVol *= env().getDim(mu);
+    }
+    timeMom.resize(par().timeMom.size());
+    for (unsigned int p = 0; p < timeMom.size(); ++p)
+    {
+        for (auto &tms: par().timeMom[p])
+        {
+            std::vector<int> tm = strToVec<int>(tms);
+            
+            timeMom[p].push_back(tm);
+            timeMomSet.insert(tm);
+        }
+        transferMom.push_back(std::vector<int>(nd - 1, 0));
+        for (auto &tm: timeMom[p])
+        {
+            for (unsigned int j = 1; j < nd; ++j)
+            {
+                transferMom[p][j - 1] -= tm[j];
+            }
+        }
+        LOG(Message) << "Probe " << p << " (" << timeMom[p].size() << " points) : " << std::endl;
+        LOG(Message) << "  phi(t_i, p_i) for (t_i, p_i) in " << timeMom[p] << std::endl;
+        LOG(Message) << "  operator with momentum " << transferMom[p] << std::endl;
+    }
+    LOG(Message) << "FFT: field '" << par().field << "'" << std::endl;
+    fft.FFT_dim_mask(ftMatBuf, phi, dMask, FFT::forward);
+    slicedProbe.resize(timeMom.size());
+    for (unsigned int p = 0; p < timeMom.size(); ++p)
+    {
+        std::vector<int> qt;
+
+        LOG(Message) << "Making probe " << p << std::endl;
+        slicedProbe[p].resize(nt);
+        for (unsigned int t = 0; t < nt; ++t)
+        {
+            Site acc;
+            
+            for (unsigned int i = 0; i < timeMom[p].size(); ++i)
+            {
+                Site buf;
+
+                qt     = timeMom[p][i];
+                qt[0] += t;
+                vectorModulo(qt);
+                peekSite(buf, ftMatBuf, qt);
+                if (i == 0)
+                {
+                    acc = buf;
+                }
+                else
+                {
+                    acc *= buf;
+                }
+            }
+            slicedProbe[p][t] = TensorRemove(trace(acc));
+        }
+        //std::cout << slicedProbe[p]<< std::endl;
+    }
+    for (auto &o: par().op)
+    {
+        auto &op = envGet(ComplexField, o);
+
+        slicedOp[o].resize(transferMom.size());
+        LOG(Message) << "FFT: operator '" << o << "'" << std::endl;
+        fft.FFT_dim_mask(ftBuf, op, dMask, FFT::forward);
+        //std::cout << ftBuf << std::endl;
+        for (unsigned int p = 0; p < transferMom.size(); ++p)
+        {
+            std::vector<int> qt(nd, 0);
+
+            for (unsigned int j = 1; j < nd; ++j)
+            {
+                qt[j] = transferMom[p][j - 1];
+            }
+            slicedOp[o][p].resize(nt);
+            for (unsigned int t = 0; t < nt; ++t)
+            {
+                TComplex buf;
+
+                qt[0] = t;
+                vectorModulo(qt);
+                peekSite(buf, ftBuf, qt);
+                slicedOp[o][p][t] = TensorRemove(buf);
+            }
+            //std::cout << ftBuf << std::endl;
+            //std::cout << slicedOp[o][p] << std::endl;
+        }
+    }
+    LOG(Message) << "Making correlators" << std::endl;
+    for (auto &o: par().op)
+    for (unsigned int p = 0; p < timeMom.size(); ++p)
+    {
+        TimeMomProbeResult r;
+
+        LOG(Message) << "  <" << o << " probe_" << p << ">" << std::endl;
+        r.op      = o;
+        r.timeMom = timeMom[p];
+        r.data    = makeTwoPoint(slicedOp[o][p], slicedProbe[p], 1./partVol);
+        result.push_back(r);
+    }
+    saveResult(par().output, "timemomprobe", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_TimeMomProbe_hpp_

extras/Hadrons/Modules/MScalarSUN/TrKinetic.cc

@@ -0,0 +1,39 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/TrKinetic.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TTrKinetic<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TTrKinetic<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TTrKinetic<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TTrKinetic<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TTrKinetic<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MScalarSUN/TrKinetic.hpp

@@ -49,24 +49,25 @@ public:
                                     std::string,  output);
 };
 
+class TrKineticResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TrKineticResult,
+                                    std::vector<std::vector<Complex>>, value,
+                                    DiffType,                          type);
+};
+
 template <typename SImpl>
 class TTrKinetic: public Module<TrKineticPar>
 {
 public:
     typedef typename SImpl::Field        Field;
     typedef typename SImpl::ComplexField ComplexField;
-    class Result: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        std::string, op,
-                                        Complex    , value);
-    };
 public:
     // constructor
     TTrKinetic(const std::string name);
     // destructor
-    virtual ~TTrKinetic(void) = default;
+    virtual ~TTrKinetic(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -76,11 +77,11 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(TrKineticSU2, TTrKinetic<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_NS(TrKineticSU3, TTrKinetic<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_NS(TrKineticSU4, TTrKinetic<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_NS(TrKineticSU5, TTrKinetic<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_NS(TrKineticSU6, TTrKinetic<ScalarNxNAdjImplR<6>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrKineticSU2, TTrKinetic<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrKineticSU3, TTrKinetic<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrKineticSU4, TTrKinetic<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrKineticSU5, TTrKinetic<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrKineticSU6, TTrKinetic<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
  *                      TTrKinetic implementation                             *
@@ -110,6 +111,7 @@ std::vector<std::string> TTrKinetic<SImpl>::getOutput(void)
     {
         out.push_back(varName(getName(), mu, nu));
     }
+    out.push_back(varName(getName(), "sum"));
     
     return out;
 }
@@ -123,6 +125,7 @@ void TTrKinetic<SImpl>::setup(void)
     {
         envCreateLat(ComplexField, varName(getName(), mu, nu));
     }
+    envCreateLat(ComplexField, varName(getName(), "sum"));
     envTmp(std::vector<Field>, "der", 1, env().getNd(), env().getGrid());
 }
 
@@ -133,34 +136,39 @@ void TTrKinetic<SImpl>::execute(void)
     LOG(Message) << "Computing tr(d_mu phi*d_nu phi) using " << par().type
                  << " derivative" << std::endl; 
 
-    std::vector<Result> result;
-    auto                &phi = envGet(Field, par().field);
+    const unsigned int nd = env().getNd();
+    TrKineticResult    result;
+    auto               &phi    = envGet(Field, par().field);
+    auto               &sumkin = envGet(ComplexField, varName(getName(), "sum"));
 
     envGetTmp(std::vector<Field>, der);
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    sumkin = zero;
+    if (!par().output.empty())
+    {
+        result.type = par().type;
+        result.value.resize(nd, std::vector<Complex>(nd));
+    }
+    for (unsigned int mu = 0; mu < nd; ++mu)
     {
         dmu(der[mu], phi, mu, par().type);
     }
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    for (unsigned int mu = 0; mu < nd; ++mu)
+    for (unsigned int nu = mu; nu < nd; ++nu)
     {
         auto &out = envGet(ComplexField, varName(getName(), mu, nu));
 
         out = -trace(der[mu]*der[nu]);
+        if (mu == nu)
+        {
+            sumkin += out;
+        }
         if (!par().output.empty())
         {
-            Result r;
-
-            r.op    = "tr(d_" + std::to_string(mu) + "phi*d_" 
-                      + std::to_string(nu) + "phi)";
-            r.value = TensorRemove(sum(out));
-            result.push_back(r);
+            result.value[mu][nu] = TensorRemove(sum(out));
+            result.value[mu][nu] = result.value[nu][mu];
         }
     }
-    if (result.size() > 0)
-    {
-        saveResult(par().output, "trkinetic", result);
-    }
+    saveResult(par().output, "trkinetic", result);
 }
 
 END_MODULE_NAMESPACE

extras/Hadrons/Modules/MScalarSUN/TrMag.cc

@@ -0,0 +1,39 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/TrMag.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/TrMag.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TTrMag<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TTrMag<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TTrMag<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TTrMag<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TTrMag<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MScalarSUN/TrMag.hpp

@@ -49,24 +49,25 @@ public:
                                     std::string,  output);
 };
 
+class TrMagResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TrMagResult,
+                                    std::string, op,
+                                    Real,        value);
+};
+
 template <typename SImpl>
 class TTrMag: public Module<TrMagPar>
 {
 public:
     typedef typename SImpl::Field        Field;
     typedef typename SImpl::ComplexField ComplexField;
-    class Result: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        std::string, op,
-                                        Real,        value);
-    };
 public:
     // constructor
     TTrMag(const std::string name);
     // destructor
-    virtual ~TTrMag(void) = default;
+    virtual ~TTrMag(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -76,11 +77,11 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(TrMagSU2, TTrMag<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_NS(TrMagSU3, TTrMag<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_NS(TrMagSU4, TTrMag<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_NS(TrMagSU5, TTrMag<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_NS(TrMagSU6, TTrMag<ScalarNxNAdjImplR<6>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrMagSU2, TTrMag<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrMagSU3, TTrMag<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrMagSU4, TTrMag<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrMagSU5, TTrMag<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrMagSU6, TTrMag<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
  *                         TTrMag implementation                              *
@@ -120,8 +121,8 @@ void TTrMag<SImpl>::execute(void)
     LOG(Message) << "Computing tr(mag^n) for n even up to " << par().maxPow
                  << std::endl;
 
-    std::vector<Result> result;
-    auto                &phi = envGet(Field, par().field);
+    std::vector<TrMagResult> result;
+    auto                     &phi = envGet(Field, par().field);
 
     auto m2 = sum(phi), mn = m2;
 
@@ -129,7 +130,7 @@ void TTrMag<SImpl>::execute(void)
     mn = 1.;
     for (unsigned int n = 2; n <= par().maxPow; n += 2)
     {
-        Result r;
+        TrMagResult r;
 
         mn = mn*m2;
         r.op    = "tr(mag^" + std::to_string(n) + ")";

extras/Hadrons/Modules/MScalarSUN/TrPhi.cc

@@ -0,0 +1,39 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/TrPhi.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/TrPhi.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TTrPhi<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TTrPhi<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TTrPhi<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TTrPhi<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TTrPhi<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MScalarSUN/TrPhi.hpp

@@ -49,24 +49,26 @@ public:
                                     std::string,  output);
 };
 
+class TrPhiResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TrPhiResult,
+                                    std::string, op,
+                                    Real,        value);
+};
+
 template <typename SImpl>
 class TTrPhi: public Module<TrPhiPar>
 {
 public:
     typedef typename SImpl::Field        Field;
     typedef typename SImpl::ComplexField ComplexField;
-    class Result: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        std::string, op,
-                                        Real,        value);
-    };
+
 public:
     // constructor
     TTrPhi(const std::string name);
     // destructor
-    virtual ~TTrPhi(void) = default;
+    virtual ~TTrPhi(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -76,11 +78,11 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(TrPhiSU2, TTrPhi<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_NS(TrPhiSU3, TTrPhi<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_NS(TrPhiSU4, TTrPhi<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_NS(TrPhiSU5, TTrPhi<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_NS(TrPhiSU6, TTrPhi<ScalarNxNAdjImplR<6>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrPhiSU2, TTrPhi<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrPhiSU3, TTrPhi<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrPhiSU4, TTrPhi<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrPhiSU5, TTrPhi<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(TrPhiSU6, TTrPhi<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
  *                          TTrPhi implementation                             *
@@ -119,7 +121,7 @@ void TTrPhi<SImpl>::setup(void)
 {
     if (par().maxPow < 2)
     {
-        HADRON_ERROR(Size, "'maxPow' should be at least equal to 2");
+        HADRONS_ERROR(Size, "'maxPow' should be at least equal to 2");
     }
     envTmpLat(Field, "phi2");
     envTmpLat(Field, "buf");
@@ -136,8 +138,8 @@ void TTrPhi<SImpl>::execute(void)
     LOG(Message) << "Computing tr(phi^n) for n even up to " << par().maxPow
                  << std::endl; 
 
-    std::vector<Result> result;
-    auto                &phi = envGet(Field, par().field);
+    std::vector<TrPhiResult> result;
+    auto                     &phi = envGet(Field, par().field);
 
     envGetTmp(Field, phi2);
     envGetTmp(Field, buf);
@@ -151,7 +153,7 @@ void TTrPhi<SImpl>::execute(void)
         phin = trace(buf);
         if (!par().output.empty())
         {
-            Result r;
+            TrPhiResult r;
 
             r.op    = "tr(phi^" + std::to_string(n) + ")";
             r.value = TensorRemove(sum(phin)).real();

extras/Hadrons/Modules/MScalarSUN/TransProj.cc

@@ -0,0 +1,39 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/TransProj.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/TransProj.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TTransProj<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TTransProj<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TTransProj<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TTransProj<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TTransProj<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MScalarSUN/TransProj.hpp

@@ -49,24 +49,25 @@ public:
                                     std::string,  output);
 };
 
+class TransProjResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TransProjResult,
+                                    std::vector<std::vector<Complex>>, value,
+                                    DiffType,                          type);
+};
+
 template <typename SImpl>
 class TTransProj: public Module<TransProjPar>
 {
 public:
     typedef typename SImpl::Field        Field;
     typedef typename SImpl::ComplexField ComplexField;
-    class Result: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        std::string, op,
-                                        Complex    , value);
-    };
 public:
     // constructor
     TTransProj(const std::string name);
     // destructor
-    virtual ~TTransProj(void) = default;
+    virtual ~TTransProj(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -77,11 +78,11 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(TransProjSU2, TTransProj<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_NS(TransProjSU3, TTransProj<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_NS(TransProjSU4, TTransProj<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_NS(TransProjSU5, TTransProj<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_NS(TransProjSU6, TTransProj<ScalarNxNAdjImplR<6>>, MScalarSUN);
+MODULE_REGISTER_TMP(TransProjSU2, TTransProj<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(TransProjSU3, TTransProj<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(TransProjSU4, TTransProj<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(TransProjSU5, TTransProj<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(TransProjSU6, TTransProj<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
  *                        TTransProj implementation                           *
@@ -137,21 +138,27 @@ void TTransProj<SImpl>::execute(void)
                  << par().type << " derivatives and op= '" << par().op 
                  << "'" << std::endl; 
 
-    std::vector<Result> result;
-    auto                &op = envGet(ComplexField, par().op);
+    const unsigned int nd = env().getNd();
+    TransProjResult    result;
+    auto               &op = envGet(ComplexField, par().op);
 
     envGetTmp(ComplexField, buf1);
     envGetTmp(ComplexField, buf2);
     envGetTmp(ComplexField, lap);
     lap = zero;
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
+    if (!par().output.empty())
+    {
+        result.type = par().type;
+        result.value.resize(nd, std::vector<Complex>(nd));
+    }
+    for (unsigned int mu = 0; mu < nd; ++mu)
     {
         dmu(buf1, op, mu, par().type);
         dmu(buf2, buf1, mu, par().type);
         lap += buf2;
     }
-    for (unsigned int mu = 0; mu < env().getNd(); ++mu)
-    for (unsigned int nu = mu; nu < env().getNd(); ++nu)
+    for (unsigned int mu = 0; mu < nd; ++mu)
+    for (unsigned int nu = mu; nu < nd; ++nu)
     {
         auto &out = envGet(ComplexField, varName(getName(), mu, nu));
         dmu(buf1, op, mu, par().type);
@@ -163,16 +170,11 @@ void TTransProj<SImpl>::execute(void)
         }
         if (!par().output.empty())
         {
-            Result r;
-
-            r.op    = "(delta_" + std::to_string(mu) + "," + std::to_string(nu)
-                      + " d^2 - d_" + std::to_string(mu) + "*d_" 
-                      + std::to_string(nu) + ")*op";
-            r.value = TensorRemove(sum(out));
-            result.push_back(r);
+            result.value[mu][nu] = TensorRemove(sum(out));
+            result.value[mu][nu] = result.value[nu][mu];
         }
     }
-    if (result.size() > 0)
+    if (!par().output.empty())
     {
         saveResult(par().output, "transproj", result);
     }

extras/Hadrons/Modules/MScalarSUN/TwoPoint.cc

@@ -0,0 +1,39 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MScalarSUN/TwoPoint.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TTwoPoint<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TTwoPoint<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TTwoPoint<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TTwoPoint<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TTwoPoint<ScalarNxNAdjImplR<6>>;
+

extras/Hadrons/Modules/MScalarSUN/TwoPoint.hpp

@@ -43,30 +43,35 @@ BEGIN_MODULE_NAMESPACE(MScalarSUN)
 class TwoPointPar: Serializable
 {
 public:
+    typedef std::pair<std::string, std::string> OpPair;
     GRID_SERIALIZABLE_CLASS_MEMBERS(TwoPointPar,
-                                    std::vector<std::string>, op,
+                                    std::vector<OpPair>,      op,
+                                    std::vector<std::string>, mom,
                                     std::string,              output);
 };
 
+class TwoPointResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TwoPointResult,
+                                    std::string, sink,
+                                    std::string, source,
+                                    std::vector<int>, mom,
+                                    std::vector<Complex>, data);
+};
+
 template <typename SImpl>
 class TTwoPoint: public Module<TwoPointPar>
 {
 public:
-    typedef typename SImpl::Field        Field;
-    typedef typename SImpl::ComplexField ComplexField;
-    class Result: Serializable
-    {
-    public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        std::string, sink,
-                                        std::string, source,
-                                        std::vector<Complex>, data);
-    };
+    typedef typename SImpl::Field         Field;
+    typedef typename SImpl::ComplexField  ComplexField;
+    typedef          std::vector<Complex> SlicedOp;
 public:
     // constructor
     TTwoPoint(const std::string name);
     // destructor
-    virtual ~TTwoPoint(void) = default;
+    virtual ~TTwoPoint(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -75,17 +80,14 @@ public:
     // execution
     virtual void execute(void);
 private:
-    // make 2-pt function
-    template <class SinkSite, class SourceSite>
-    std::vector<Complex> makeTwoPoint(const std::vector<SinkSite>   &sink,
-                                      const std::vector<SourceSite> &source);
+    std::vector<std::vector<int>> mom_;
 };
 
-MODULE_REGISTER_NS(TwoPointSU2, TTwoPoint<ScalarNxNAdjImplR<2>>, MScalarSUN);
-MODULE_REGISTER_NS(TwoPointSU3, TTwoPoint<ScalarNxNAdjImplR<3>>, MScalarSUN);
-MODULE_REGISTER_NS(TwoPointSU4, TTwoPoint<ScalarNxNAdjImplR<4>>, MScalarSUN);
-MODULE_REGISTER_NS(TwoPointSU5, TTwoPoint<ScalarNxNAdjImplR<5>>, MScalarSUN);
-MODULE_REGISTER_NS(TwoPointSU6, TTwoPoint<ScalarNxNAdjImplR<6>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointSU2, TTwoPoint<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointSU3, TTwoPoint<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointSU4, TTwoPoint<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointSU5, TTwoPoint<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointSU6, TTwoPoint<ScalarNxNAdjImplR<6>>, MScalarSUN);
 
 /******************************************************************************
  *                       TTwoPoint implementation                             *
@@ -100,7 +102,20 @@ TTwoPoint<SImpl>::TTwoPoint(const std::string name)
 template <typename SImpl>
 std::vector<std::string> TTwoPoint<SImpl>::getInput(void)
 {   
-    return par().op;
+    std::vector<std::string> in;
+    std::set<std::string>    ops;
+
+    for (auto &p: par().op)
+    {
+        ops.insert(p.first);
+        ops.insert(p.second);
+    }
+    for (auto &o: ops)
+    {
+        in.push_back(o);
+    }
+
+    return in;
 }
 
 template <typename SImpl>
@@ -115,68 +130,91 @@ std::vector<std::string> TTwoPoint<SImpl>::getOutput(void)
 template <typename SImpl>
 void TTwoPoint<SImpl>::setup(void)
 {
-    const unsigned int nt = env().getDim().back();
-    envTmp(std::vector<std::vector<TComplex>>, "slicedOp", 1, par().op.size(), 
-           std::vector<TComplex>(nt));
+    const unsigned int nd = env().getDim().size();
+
+    mom_.resize(par().mom.size());
+    for (unsigned int i = 0; i < mom_.size(); ++i)
+    {
+        mom_[i] = strToVec<int>(par().mom[i]);
+        if (mom_[i].size() != nd - 1)
+        {
+            HADRONS_ERROR(Size, "momentum number of components different from " 
+                               + std::to_string(nd-1));
+        }
+    }
+    envTmpLat(ComplexField, "ftBuf");
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename SImpl>
 void TTwoPoint<SImpl>::execute(void)
 {
-    LOG(Message) << "Computing 2-point functions for operators:" << std::endl;
-    for (auto &o: par().op)
+    LOG(Message) << "Computing 2-point functions" << std::endl;
+    for (auto &p: par().op)
     {
-        LOG(Message) << "  '" << o << "'" << std::endl;
+        LOG(Message) << "  <" << p.first << " " << p.second << ">" << std::endl;
     }
 
-    const unsigned int  nd = env().getDim().size();
-    std::vector<Result> result;
-    
-    envGetTmp(std::vector<std::vector<TComplex>>, slicedOp);
-    for (unsigned int i = 0; i < par().op.size(); ++i)
-    {
-        auto &op = envGet(ComplexField, par().op[i]);
+    const unsigned int                           nd      = env().getNd();
+    const unsigned int                           nt      = env().getDim().back();
+    const unsigned int                           nop     = par().op.size();
+    const unsigned int                           nmom    = mom_.size();
+    double                                       partVol = 1.;
+    std::vector<int>                             dMask(nd, 1);
+    std::set<std::string>                        ops;
+    std::vector<TwoPointResult>                  result;
+    std::map<std::string, std::vector<SlicedOp>> slicedOp;
+    FFT                                          fft(env().getGrid());
+    TComplex                                     buf;
 
-        sliceSum(op, slicedOp[i], nd - 1);
+    envGetTmp(ComplexField, ftBuf);
+    dMask[nd - 1] = 0;
+    for (unsigned int mu = 0; mu < nd - 1; ++mu)
+    {
+        partVol *= env().getDim()[mu];
     }
-    for (unsigned int i = 0; i < par().op.size(); ++i)
-    for (unsigned int j = 0; j < par().op.size(); ++j)
+    for (auto &p: par().op)
     {
-        Result r;
+        ops.insert(p.first);
+        ops.insert(p.second);
+    }
+    for (auto &o: ops)
+    {
+        auto &op = envGet(ComplexField, o);
 
-        r.sink   = par().op[i];
-        r.source = par().op[j];
-        r.data   = makeTwoPoint(slicedOp[i], slicedOp[j]);
+        slicedOp[o].resize(nmom);
+        LOG(Message) << "Operator '" << o << "' FFT" << std::endl;
+        fft.FFT_dim_mask(ftBuf, op, dMask, FFT::forward);
+        for (unsigned int m = 0; m < nmom; ++m)
+        {
+            auto qt = mom_[m];
+
+            qt.resize(nd);
+            slicedOp[o][m].resize(nt);
+            for (unsigned int t = 0; t < nt; ++t)
+            {
+                qt[nd - 1] = t;
+                peekSite(buf, ftBuf, qt);
+                slicedOp[o][m][t] = TensorRemove(buf);
+            }
+        }
+    }
+    LOG(Message) << "Making contractions" << std::endl;
+    for (unsigned int m = 0; m < nmom; ++m)
+    for (auto &p: par().op)
+    {
+        TwoPointResult r;
+
+        r.sink   = p.first;
+        r.source = p.second;
+        r.mom    = mom_[m];
+        r.data   = makeTwoPoint(slicedOp[p.first][m], slicedOp[p.second][m], 
+                                1./partVol);
         result.push_back(r);
     }
     saveResult(par().output, "twopt", result);
 }
 
-// make 2-pt function //////////////////////////////////////////////////////////
-template <class SImpl>
-template <class SinkSite, class SourceSite>
-std::vector<Complex> TTwoPoint<SImpl>::makeTwoPoint(
-                                  const std::vector<SinkSite>   &sink,
-                                  const std::vector<SourceSite> &source)
-{
-    assert(sink.size() == source.size());
-    
-    unsigned int         nt = sink.size();
-    std::vector<Complex> res(nt, 0.);
-    
-    for (unsigned int dt = 0; dt < nt; ++dt)
-    {
-        for (unsigned int t  = 0; t < nt; ++t)
-        {
-            res[dt] += TensorRemove(trace(sink[(t+dt)%nt]*source[t]));
-        }
-        res[dt] *= 1./static_cast<double>(nt);
-    }
-    
-    return res;
-}
-
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE

extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.cc

@@ -0,0 +1,11 @@
+#include <Grid/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MScalarSUN;
+
+template class Grid::Hadrons::MScalarSUN::TTwoPointNPR<ScalarNxNAdjImplR<2>>;
+template class Grid::Hadrons::MScalarSUN::TTwoPointNPR<ScalarNxNAdjImplR<3>>;
+template class Grid::Hadrons::MScalarSUN::TTwoPointNPR<ScalarNxNAdjImplR<4>>;
+template class Grid::Hadrons::MScalarSUN::TTwoPointNPR<ScalarNxNAdjImplR<5>>;
+template class Grid::Hadrons::MScalarSUN::TTwoPointNPR<ScalarNxNAdjImplR<6>>;

extras/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp

@@ -0,0 +1,182 @@
+#ifndef Hadrons_MScalarSUN_TwoPointNPR_hpp_
+#define Hadrons_MScalarSUN_TwoPointNPR_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+#include <Grid/Hadrons/Modules/MScalarSUN/Utils.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         TwoPointNPR                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MScalarSUN)
+
+class TwoPointNPRPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TwoPointNPRPar,
+                                    std::vector<std::string>, op,
+                                    std::string,              field,
+                                    std::string,              output);
+};
+
+class TwoPointNPRResult: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(TwoPointNPRResult,
+                                    std::string,          op,
+                                    std::vector<Complex>, data);
+};
+
+template <typename SImpl>
+class TTwoPointNPR: public Module<TwoPointNPRPar>
+{
+public:
+    typedef typename SImpl::Field                    Field;
+    typedef typename SImpl::SiteField::scalar_object Site;
+    typedef typename SImpl::ComplexField             ComplexField;
+public:
+    // constructor
+    TTwoPointNPR(const std::string name);
+    // destructor
+    virtual ~TTwoPointNPR(void) {};
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_TMP(TwoPointNPRSU2, TTwoPointNPR<ScalarNxNAdjImplR<2>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointNPRSU3, TTwoPointNPR<ScalarNxNAdjImplR<3>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointNPRSU4, TTwoPointNPR<ScalarNxNAdjImplR<4>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointNPRSU5, TTwoPointNPR<ScalarNxNAdjImplR<5>>, MScalarSUN);
+MODULE_REGISTER_TMP(TwoPointNPRSU6, TTwoPointNPR<ScalarNxNAdjImplR<6>>, MScalarSUN);
+
+/******************************************************************************
+ *                 TTwoPointNPR implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename SImpl>
+TTwoPointNPR<SImpl>::TTwoPointNPR(const std::string name)
+: Module<TwoPointNPRPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename SImpl>
+std::vector<std::string> TTwoPointNPR<SImpl>::getInput(void)
+{
+    std::vector<std::string> in = par().op;
+
+    in.push_back(par().field);
+
+    return in;
+}
+
+template <typename SImpl>
+std::vector<std::string> TTwoPointNPR<SImpl>::getOutput(void)
+{
+    std::vector<std::string> out;
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TTwoPointNPR<SImpl>::setup(void)
+{
+    const unsigned int nl = env().getDim(0);
+
+    for (unsigned int mu = 1; mu < env().getNd(); ++mu)
+    {
+        if (nl != env().getDim(mu))
+        {
+            HADRONS_ERROR(Size, "non-cubic grid");
+        }
+    }
+    envTmpLat(ComplexField, "ftBuf");
+    envTmpLat(Field, "ftMatBuf");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename SImpl>
+void TTwoPointNPR<SImpl>::execute(void)
+{
+    const unsigned int             nd = env().getNd();
+    const unsigned int             nl = env().getDim(0);
+    FFT                            fft(env().getGrid());
+    std::vector<TwoPointNPRResult> result;
+    TwoPointNPRResult              twoPtp1, twoPtp2;
+    auto                           &phi    = envGet(Field, par().field);
+    bool                           doTwoPt = true;
+
+    envGetTmp(ComplexField, ftBuf);
+    envGetTmp(Field, ftMatBuf);
+    LOG(Message) << "FFT: field '" << par().field << "'" << std::endl;
+    fft.FFT_all_dim(ftMatBuf, phi, FFT::forward);
+    for (auto &opName: par().op)
+    {
+        auto              &op = envGet(ComplexField, opName);
+        std::vector<int>  p1, p2, p;
+        Site              phip1, phip2;
+        TComplex          opp;
+        TwoPointNPRResult r;
+
+        LOG(Message) << "FFT: operator '" << opName << "'" << std::endl;
+        fft.FFT_all_dim(ftBuf, op, FFT::forward);
+        LOG(Message) << "Generating vertex function" << std::endl;
+        r.op = opName;
+        r.data.resize(nl);
+        if (doTwoPt)
+        {
+            twoPtp1.op = "phi_prop_p1";
+            twoPtp1.data.resize(nl);
+            twoPtp2.op = "phi_prop_p2";
+            twoPtp2.data.resize(nl);
+        }
+        for (unsigned int n = 0; n < nl; ++n)
+        {
+            p1.assign(nd, 0);
+            p2.assign(nd, 0);
+            p.assign(nd, 0);
+            // non-exceptional RI/SMOM kinematic
+            // p1 = mu*(1,1,0): in mom
+            // p2 = mu*(0,1,1): out mom
+            // p  = p1 - p2 = mu*(1,0,-1)
+            // mu = 2*n*pi/L
+            p1[0] = n;
+            p1[1] = n;
+            p2[1] = n;
+            p2[2] = n;
+            p[0]  = n;
+            p[2]  = (nl - n) % nl;
+            peekSite(phip1, ftMatBuf, p1);
+            peekSite(phip2, ftMatBuf, p2);
+            peekSite(opp, ftBuf, p);
+            if (doTwoPt)
+            {
+                twoPtp1.data[n] = TensorRemove(trace(phip1*adj(phip1)));
+                twoPtp2.data[n] = TensorRemove(trace(phip2*adj(phip2)));
+            }
+            r.data[n]  = TensorRemove(trace(phip2*adj(phip1))*opp);
+        }
+        if (doTwoPt)
+        {
+            result.push_back(twoPtp1);
+            result.push_back(twoPtp2);
+        }
+        result.push_back(r);
+        doTwoPt = false;
+    }
+    saveResult(par().output, "twoptnpr", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MScalarSUN_TwoPointNPR_hpp_

extras/Hadrons/Modules/MScalarSUN/Utils.hpp

@@ -44,7 +44,7 @@ inline void dmu(Field &out, const Field &in, const unsigned int mu, const DiffTy
 
     if (mu >= env.getNd())
     {
-        HADRON_ERROR(Range, "Derivative direction out of range");
+        HADRONS_ERROR(Range, "Derivative direction out of range");
     }
     switch(type)
     {
@@ -58,7 +58,7 @@ inline void dmu(Field &out, const Field &in, const unsigned int mu, const DiffTy
             out = 0.5*(Cshift(in, mu, 1) - Cshift(in, mu, -1));
             break;
         default:
-            HADRON_ERROR(Argument, "Derivative type invalid");
+            HADRONS_ERROR(Argument, "Derivative type invalid");
             break;
     }
 }
@@ -70,7 +70,7 @@ inline void dmuAcc(Field &out, const Field &in, const unsigned int mu, const Dif
 
     if (mu >= env.getNd())
     {
-        HADRON_ERROR(Range, "Derivative direction out of range");
+        HADRONS_ERROR(Range, "Derivative direction out of range");
     }
     switch(type)
     {
@@ -84,20 +84,47 @@ inline void dmuAcc(Field &out, const Field &in, const unsigned int mu, const Dif
             out += 0.5*(Cshift(in, mu, 1) - Cshift(in, mu, -1));
             break;
         default:
-            HADRON_ERROR(Argument, "Derivative type invalid");
+            HADRONS_ERROR(Argument, "Derivative type invalid");
             break;
     }
 }
 
+template <class SinkSite, class SourceSite>
+std::vector<Complex> makeTwoPoint(const std::vector<SinkSite>   &sink,
+                                  const std::vector<SourceSite> &source,
+                                  const double factor = 1.)
+{
+    assert(sink.size() == source.size());
+    
+    unsigned int         nt = sink.size();
+    std::vector<Complex> res(nt, 0.);
+    
+    for (unsigned int dt = 0; dt < nt; ++dt)
+    {
+        for (unsigned int t  = 0; t < nt; ++t)
+        {
+            res[dt] += trace(sink[(t+dt)%nt]*source[t]);
+        }
+        res[dt] *= factor/static_cast<double>(nt);
+    }
+    
+    return res;
+}
+
+inline std::string varName(const std::string name, const std::string suf)
+{
+    return name + "_" + suf;
+}
+
 inline std::string varName(const std::string name, const unsigned int mu)
 {
-    return name + "_" + std::to_string(mu);
+    return varName(name, std::to_string(mu));
 }
 
 inline std::string varName(const std::string name, const unsigned int mu, 
                            const unsigned int nu)
 {
-    return name + "_" + std::to_string(mu) + "_" + std::to_string(nu);
+    return varName(name, std::to_string(mu) + "_" + std::to_string(nu));
 }
 
 END_MODULE_NAMESPACE

extras/Hadrons/Modules/MSink/Point.cc

@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSink/Point.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MSink/Point.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSink;
+
+template class Grid::Hadrons::MSink::TPoint<FIMPL>;
+template class Grid::Hadrons::MSink::TPoint<ScalarImplCR>;
+

extras/Hadrons/Modules/MSink/Point.hpp

@@ -7,6 +7,7 @@ Source file: extras/Hadrons/Modules/MSink/Point.hpp
 Copyright (C) 2015-2018
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Guido Cossu <guido.cossu@ed.ac.uk>
 Author: Lanny91 <andrew.lawson@gmail.com>
 
 This program is free software; you can redistribute it and/or modify
@@ -58,7 +59,7 @@ public:
     // constructor
     TPoint(const std::string name);
     // destructor
-    virtual ~TPoint(void) = default;
+    virtual ~TPoint(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -72,8 +73,8 @@ private:
     std::string momphName_;
 };
 
-MODULE_REGISTER_NS(Point,       TPoint<FIMPL>,        MSink);
-MODULE_REGISTER_NS(ScalarPoint, TPoint<ScalarImplCR>, MSink);
+MODULE_REGISTER_TMP(Point,       TPoint<FIMPL>,        MSink);
+MODULE_REGISTER_TMP(ScalarPoint, TPoint<ScalarImplCR>, MSink);
 
 /******************************************************************************
  *                          TPoint implementation                             *

extras/Hadrons/Modules/MSink/Smear.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSink/Smear.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MSink/Smear.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSink;
+
+template class Grid::Hadrons::MSink::TSmear<FIMPL>;
+

extras/Hadrons/Modules/MSink/Smear.hpp

@@ -59,7 +59,7 @@ public:
     // constructor
     TSmear(const std::string name);
     // destructor
-    virtual ~TSmear(void) = default;
+    virtual ~TSmear(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -70,7 +70,7 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(Smear, TSmear<FIMPL>, MSink);
+MODULE_REGISTER_TMP(Smear, TSmear<FIMPL>, MSink);
 
 /******************************************************************************
  *                          TSmear implementation                             *

extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc

@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSolver;
+
+template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos<FIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>;
+template class Grid::Hadrons::MSolver::TLocalCoherenceLanczos<ZFIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>;
+

extras/Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp

@@ -51,7 +51,8 @@ public:
                                     ChebyParams,   smoother,
                                     RealD,         coarseRelaxTol,
                                     std::string,   blockSize,
-                                    std::string,   output);
+                                    std::string,   output,
+                                    bool,          multiFile);
 };
 
 template <typename FImpl, int nBasis>
@@ -69,7 +70,7 @@ public:
     // constructor
     TLocalCoherenceLanczos(const std::string name);
     // destructor
-    virtual ~TLocalCoherenceLanczos(void) = default;
+    virtual ~TLocalCoherenceLanczos(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -79,12 +80,8 @@ public:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(LocalCoherenceLanczos, 
-    ARG(TLocalCoherenceLanczos<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), 
-    MSolver);
-MODULE_REGISTER_NS(ZLocalCoherenceLanczos, 
-    ARG(TLocalCoherenceLanczos<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), 
-    MSolver);
+MODULE_REGISTER_TMP(LocalCoherenceLanczos, ARG(TLocalCoherenceLanczos<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
+MODULE_REGISTER_TMP(ZLocalCoherenceLanczos, ARG(TLocalCoherenceLanczos<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
 
 /******************************************************************************
  *                 TLocalCoherenceLanczos implementation                      *
@@ -125,19 +122,18 @@ void TLocalCoherenceLanczos<FImpl, nBasis>::setup(void)
 
     env().createCoarseGrid(blockSize, Ls);
 
-    auto cg   = env().getCoarseGrid(blockSize, Ls);
-    auto cgrb = env().getRbCoarseGrid(blockSize, Ls);
-    int  cNm  = (par().doCoarse) ? par().coarseParams.Nm : 0;
+    auto cg  = env().getCoarseGrid(blockSize, Ls);
+    int  cNm = (par().doCoarse) ? par().coarseParams.Nm : 0;
 
     LOG(Message) << "Coarse grid: " << cg->GlobalDimensions() << std::endl;
     envCreateDerived(BasePack, CoarsePack, getName(), Ls,
-                     par().fineParams.Nm, cNm, env().getRbGrid(Ls), cgrb);
+                     par().fineParams.Nm, cNm, env().getRbGrid(Ls), cg);
 
     auto &epack = envGetDerived(BasePack, CoarsePack, getName());
 
     envTmp(SchurFMat, "mat", Ls, envGet(FMat, par().action));
     envGetTmp(SchurFMat, mat);
-    envTmp(LCL, "solver", Ls, env().getRbGrid(Ls), cgrb, mat, 
+    envTmp(LCL, "solver", Ls, env().getRbGrid(Ls), cg, mat, 
            Odd, epack.evec, epack.evecCoarse, epack.eval, epack.evalCoarse);
 }
 
@@ -149,6 +145,8 @@ void TLocalCoherenceLanczos<FImpl, nBasis>::execute(void)
     auto &coarsePar = par().coarseParams;
     auto &epack     = envGetDerived(BasePack, CoarsePack, getName());
 
+    epack.record.operatorXml = vm().getModule(par().action)->parString();
+    epack.record.solverXml   = parString();
     envGetTmp(LCL, solver);
     LOG(Message) << "Performing fine grid IRL -- Nstop= " 
                  << finePar.Nstop << ", Nk= " << finePar.Nk << ", Nm= " 
@@ -157,6 +155,10 @@ void TLocalCoherenceLanczos<FImpl, nBasis>::execute(void)
                     finePar.resid,finePar.MaxIt, finePar.betastp, 
                     finePar.MinRes);
     solver.testFine(finePar.resid*100.0);
+    if (!par().output.empty())
+    {
+        epack.writeFine(par().output, par().multiFile, vm().getTrajectory());
+    }
     if (par().doCoarse)
     {
         LOG(Message) << "Orthogonalising" << std::endl;
@@ -170,10 +172,10 @@ void TLocalCoherenceLanczos<FImpl, nBasis>::execute(void)
                           coarsePar.MinRes);
         solver.testCoarse(coarsePar.resid*100.0, par().smoother, 
                         par().coarseRelaxTol);
-    }
-    if (!par().output.empty())
-    {
-        epack.write(par().output, vm().getTrajectory());
+        if (!par().output.empty())
+        {
+            epack.writeCoarse(par().output, par().multiFile, vm().getTrajectory());
+        }
     }
 }
 

extras/Hadrons/Modules/MSolver/RBPrecCG.cc

@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSolver/RBPrecCG.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSolver;
+
+template class Grid::Hadrons::MSolver::TRBPrecCG<FIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>;
+template class Grid::Hadrons::MSolver::TRBPrecCG<ZFIMPL,HADRONS_DEFAULT_LANCZOS_NBASIS>;
+

extras/Hadrons/Modules/MSolver/RBPrecCG.hpp

@@ -67,7 +67,7 @@ public:
     // constructor
     TRBPrecCG(const std::string name);
     // destructor
-    virtual ~TRBPrecCG(void) = default;
+    virtual ~TRBPrecCG(void) {};
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getReference(void);
@@ -79,10 +79,8 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(RBPrecCG,  
-    ARG(TRBPrecCG<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
-MODULE_REGISTER_NS(ZRBPrecCG, 
-    ARG(TRBPrecCG<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
+MODULE_REGISTER_TMP(RBPrecCG, ARG(TRBPrecCG<FIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
+MODULE_REGISTER_TMP(ZRBPrecCG, ARG(TRBPrecCG<ZFIMPL, HADRONS_DEFAULT_LANCZOS_NBASIS>), MSolver);
 
 /******************************************************************************
  *                      TRBPrecCG template implementation                     *
@@ -107,6 +105,11 @@ std::vector<std::string> TRBPrecCG<FImpl, nBasis>::getReference(void)
 {
     std::vector<std::string> ref = {par().action};
     
+    if (!par().eigenPack.empty())
+    {
+        ref.push_back(par().eigenPack);
+    }
+
     return ref;
 }
 
@@ -124,7 +127,7 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
 {
     if (par().maxIteration == 0)
     {
-        HADRON_ERROR(Argument, "zero maximum iteration");
+        HADRONS_ERROR(Argument, "zero maximum iteration");
     }
 
     LOG(Message) << "setting up Schur red-black preconditioned CG for"
@@ -147,6 +150,9 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
         {
             auto &epack = envGetDerived(EPack, CoarseEPack, par().eigenPack);
             
+            LOG(Message) << "using low-mode deflation with coarse eigenpack '"
+                         << par().eigenPack << "' (" 
+                         << epack.evecCoarse.size() << " modes)" << std::endl;
             guesser.reset(new CoarseGuesser(epack.evec, epack.evecCoarse,
                                             epack.evalCoarse));
         }
@@ -154,6 +160,9 @@ void TRBPrecCG<FImpl, nBasis>::setup(void)
         {
             auto &epack = envGet(EPack, par().eigenPack);
 
+            LOG(Message) << "using low-mode deflation with eigenpack '"
+                         << par().eigenPack << "' (" 
+                         << epack.evec.size() << " modes)" << std::endl;
             guesser.reset(new FineGuesser(epack.evec, epack.eval));
         }
     }

extras/Hadrons/Modules/MSource/Point.cc

@@ -0,0 +1,36 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/Point.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MSource/Point.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSource;
+
+template class Grid::Hadrons::MSource::TPoint<FIMPL>;
+template class Grid::Hadrons::MSource::TPoint<ScalarImplCR>;
+

extras/Hadrons/Modules/MSource/Point.hpp

@@ -68,7 +68,7 @@ public:
     // constructor
     TPoint(const std::string name);
     // destructor
-    virtual ~TPoint(void) = default;
+    virtual ~TPoint(void) {};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
@@ -79,8 +79,8 @@ protected:
     virtual void execute(void);
 };
 
-MODULE_REGISTER_NS(Point,       TPoint<FIMPL>,        MSource);
-MODULE_REGISTER_NS(ScalarPoint, TPoint<ScalarImplCR>, MSource);
+MODULE_REGISTER_TMP(Point,       TPoint<FIMPL>,        MSource);
+MODULE_REGISTER_TMP(ScalarPoint, TPoint<ScalarImplCR>, MSource);
 
 /******************************************************************************
  *                       TPoint template implementation                       *

extras/Hadrons/Modules/MSource/SeqConserved.cc

@@ -0,0 +1,35 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/SeqConserved.cc
+
+Copyright (C) 2015-2018
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MSource;
+
+template class Grid::Hadrons::MSource::TSeqConserved<FIMPL>;
+