Merge remote-tracking branch 'upstream/develop' into feature/kl2QED

Conflicts: Hadrons/Modules.hpp Hadrons/modules.inc
2025-11-04 05:54:32 +00:00 · 2019-03-06 11:28:59 +00:00
parent 447b772136 acd25d0d01
commit 0b426bf9f6
207 changed files with 1714 additions and 2172 deletions
--- a/Grid/Grid.h
+++ b/Grid/Grid.h
@@ -42,6 +42,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #include <Grid/GridQCDcore.h>
 #include <Grid/qcd/action/Action.h>
 #include <Grid/qcd/utils/GaugeFix.h>
+#include <Grid/qcd/utils/CovariantSmearing.h>
 #include <Grid/qcd/smearing/Smearing.h>
 #include <Grid/parallelIO/MetaData.h>
 #include <Grid/qcd/hmc/HMC_aggregate.h>
--- a/Grid/communicator/Communicator_mpi3.cc
+++ b/Grid/communicator/Communicator_mpi3.cc
@@ -53,10 +53,12 @@ void CartesianCommunicator::Init(int *argc, char ***argv)
  // Never clean up as done once.
  MPI_Comm_dup (MPI_COMM_WORLD,&communicator_world);

+  Grid_quiesce_nodes();
  GlobalSharedMemory::Init(communicator_world);
  GlobalSharedMemory::SharedMemoryAllocate(
 		   GlobalSharedMemory::MAX_MPI_SHM_BYTES,
 		   GlobalSharedMemory::Hugepages);
+  Grid_unquiesce_nodes();
 }

 ///////////////////////////////////////////////////////////////////////////
--- a/Grid/parallelIO/BinaryIO.h
+++ b/Grid/parallelIO/BinaryIO.h
@@ -210,10 +210,10 @@ PARALLEL_CRITICAL
  static inline void le32toh_v(void *file_object,uint64_t bytes)
  {
    uint32_t *fp = (uint32_t *)file_object;
-    uint32_t f;

    uint64_t count = bytes/sizeof(uint32_t);
    parallel_for(uint64_t i=0;i<count;i++){  
+      uint32_t f;
      f = fp[i];
      // got network order and the network to host
      f = ((f&0xFF)<<24) | ((f&0xFF00)<<8) | ((f&0xFF0000)>>8) | ((f&0xFF000000UL)>>24) ; 
@@ -235,10 +235,9 @@ PARALLEL_CRITICAL
  static inline void le64toh_v(void *file_object,uint64_t bytes)
  {
    uint64_t *fp = (uint64_t *)file_object;
-    uint64_t f,g;
-    
    uint64_t count = bytes/sizeof(uint64_t);
    parallel_for(uint64_t i=0;i<count;i++){  
+      uint64_t f,g;
      f = fp[i];
      // got network order and the network to host
      g = ((f&0xFF)<<24) | ((f&0xFF00)<<8) | ((f&0xFF0000)>>8) | ((f&0xFF000000UL)>>24) ; 
@@ -349,7 +348,8 @@ PARALLEL_CRITICAL
    int ieee32    = (format == std::string("IEEE32"));
    int ieee64big = (format == std::string("IEEE64BIG"));
    int ieee64    = (format == std::string("IEEE64"));
-
+    assert(ieee64||ieee32|ieee64big||ieee32big);
+    assert((ieee64+ieee32+ieee64big+ieee32big)==1);
    //////////////////////////////////////////////////////////////////////////////
    // Do the I/O
    //////////////////////////////////////////////////////////////////////////////
--- a/Grid/parallelIO/IldgIO.h
+++ b/Grid/parallelIO/IldgIO.h
@@ -46,6 +46,12 @@ extern "C" {
 namespace Grid {
 namespace QCD {

+#define GRID_FIELD_NORM "FieldNormMetaData"
+#define GRID_FIELD_NORM_CALC(FieldNormMetaData_, n2ck) \
+0.5*fabs(FieldNormMetaData_.norm2 - n2ck)/(FieldNormMetaData_.norm2 + n2ck)
+#define GRID_FIELD_NORM_CHECK(FieldNormMetaData_, n2ck) \
+assert(GRID_FIELD_NORM_CALC(FieldNormMetaData_, n2ck) < 1.0e-5);
+
  /////////////////////////////////
  // Encode word types as strings
  /////////////////////////////////
@@ -205,6 +211,7 @@ class GridLimeReader : public BinaryIO {
  {
    typedef typename vobj::scalar_object sobj;
    scidacChecksum scidacChecksum_;
+    FieldNormMetaData  FieldNormMetaData_;
    uint32_t nersc_csum,scidac_csuma,scidac_csumb;

    std::string format = getFormatString<vobj>();
@@ -238,16 +245,47 @@ class GridLimeReader : public BinaryIO {
 	/////////////////////////////////////////////
 	// Insist checksum is next record
 	/////////////////////////////////////////////
-	readLimeObject(scidacChecksum_,std::string("scidacChecksum"),std::string(SCIDAC_CHECKSUM));
-
+	readScidacChecksum(scidacChecksum_,FieldNormMetaData_);
 	/////////////////////////////////////////////
 	// Verify checksums
 	/////////////////////////////////////////////
+	if(FieldNormMetaData_.norm2 != 0.0){ 
+	  RealD n2ck = norm2(field);
+	  std::cout << GridLogMessage << "Field norm: metadata= " << FieldNormMetaData_.norm2 
+              << " / field= " << n2ck << " / rdiff= " << GRID_FIELD_NORM_CALC(FieldNormMetaData_,n2ck) << std::endl;
+	  GRID_FIELD_NORM_CHECK(FieldNormMetaData_,n2ck);
+	}
 	assert(scidacChecksumVerify(scidacChecksum_,scidac_csuma,scidac_csumb)==1);
+
+	// find out if next field is a GridFieldNorm
 	return;
      }
    }
  }
+  void readScidacChecksum(scidacChecksum     &scidacChecksum_,
+			  FieldNormMetaData  &FieldNormMetaData_)
+  {
+    FieldNormMetaData_.norm2 =0.0;
+    std::string scidac_str(SCIDAC_CHECKSUM);
+    std::string field_norm_str(GRID_FIELD_NORM);
+    while ( limeReaderNextRecord(LimeR) == LIME_SUCCESS ) { 
+      uint64_t nbytes = limeReaderBytes(LimeR);//size of this record (configuration)
+      std::vector<char> xmlc(nbytes+1,'\0');
+      limeReaderReadData((void *)&xmlc[0], &nbytes, LimeR);    
+      std::string xmlstring = std::string(&xmlc[0]);
+      XmlReader RD(xmlstring, true, "");
+      if ( !strncmp(limeReaderType(LimeR), field_norm_str.c_str(),strlen(field_norm_str.c_str()) )  ) {
+	//	std::cout << "FieldNormMetaData "<<xmlstring<<std::endl;
+	read(RD,field_norm_str,FieldNormMetaData_);
+      }
+      if ( !strncmp(limeReaderType(LimeR), scidac_str.c_str(),strlen(scidac_str.c_str()) )  ) {
+	//	std::cout << SCIDAC_CHECKSUM << " " <<xmlstring<<std::endl;
+	read(RD,std::string("scidacChecksum"),scidacChecksum_);
+	return;
+      }      
+    }
+    assert(0);
+  }
  ////////////////////////////////////////////
  // Read a generic serialisable object
  ////////////////////////////////////////////
@@ -390,6 +428,8 @@ class GridLimeWriter : public BinaryIO
    GridBase *grid = field._grid;
    assert(boss_node == field._grid->IsBoss() );

+    FieldNormMetaData FNMD; FNMD.norm2 = norm2(field);
+
    ////////////////////////////////////////////
    // Create record header
    ////////////////////////////////////////////
@@ -448,6 +488,7 @@ class GridLimeWriter : public BinaryIO
    checksum.suma= streama.str();
    checksum.sumb= streamb.str();
    if ( boss_node ) { 
+      writeLimeObject(0,0,FNMD,std::string(GRID_FIELD_NORM),std::string(GRID_FIELD_NORM));
      writeLimeObject(0,1,checksum,std::string("scidacChecksum"),std::string(SCIDAC_CHECKSUM));
    }
  }
@@ -625,6 +666,12 @@ class IldgWriter : public ScidacWriter {
    assert(header.nd==4);
    assert(header.nd==header.dimension.size());

+    //////////////////////////////////////////////////////////////////////////////
+    // Field norm tests
+    //////////////////////////////////////////////////////////////////////////////
+    FieldNormMetaData FieldNormMetaData_;
+    FieldNormMetaData_.norm2 = norm2(Umu);
+
    //////////////////////////////////////////////////////////////////////////////
    // Fill the USQCD info field
    //////////////////////////////////////////////////////////////////////////////
@@ -633,11 +680,12 @@ class IldgWriter : public ScidacWriter {
    info.plaq   = header.plaquette;
    info.linktr = header.link_trace;

-    std::cout << GridLogMessage << " Writing config; IldgIO "<<std::endl;
+    //    std::cout << GridLogMessage << " Writing config; IldgIO n2 "<< FieldNormMetaData_.norm2<<std::endl;
    //////////////////////////////////////////////
    // Fill the Lime file record by record
    //////////////////////////////////////////////
    writeLimeObject(1,0,header ,std::string("FieldMetaData"),std::string(GRID_FORMAT)); // Open message 
+    writeLimeObject(0,0,FieldNormMetaData_,FieldNormMetaData_.SerialisableClassName(),std::string(GRID_FIELD_NORM));
    writeLimeObject(0,0,_scidacFile,_scidacFile.SerialisableClassName(),std::string(SCIDAC_PRIVATE_FILE_XML));
    writeLimeObject(0,1,info,info.SerialisableClassName(),std::string(SCIDAC_FILE_XML));
    writeLimeObject(1,0,_scidacRecord,_scidacRecord.SerialisableClassName(),std::string(SCIDAC_PRIVATE_RECORD_XML));
@@ -680,6 +728,7 @@ class IldgReader : public GridLimeReader {
    std::string    ildgLFN_       ;
    scidacChecksum scidacChecksum_; 
    usqcdInfo      usqcdInfo_     ;
+    FieldNormMetaData FieldNormMetaData_;

    // track what we read from file
    int found_ildgFormat    =0;
@@ -688,7 +737,7 @@ class IldgReader : public GridLimeReader {
    int found_usqcdInfo     =0;
    int found_ildgBinary =0;
    int found_FieldMetaData =0;
-
+    int found_FieldNormMetaData =0;
    uint32_t nersc_csum;
    uint32_t scidac_csuma;
    uint32_t scidac_csumb;
@@ -775,11 +824,17 @@ class IldgReader : public GridLimeReader {
 	  found_scidacChecksum = 1;
 	}

+	if ( !strncmp(limeReaderType(LimeR), GRID_FIELD_NORM,strlen(GRID_FIELD_NORM)) ) { 
+	  XmlReader RD(xmlstring, true, "");
+	  read(RD,GRID_FIELD_NORM,FieldNormMetaData_);
+	  found_FieldNormMetaData = 1;
+	}
+
      } else {  
 	/////////////////////////////////
 	// Binary data
 	/////////////////////////////////
-	std::cout << GridLogMessage << "ILDG Binary record found : "  ILDG_BINARY_DATA << std::endl;
+	//	std::cout << GridLogMessage << "ILDG Binary record found : "  ILDG_BINARY_DATA << std::endl;
 	uint64_t offset= ftello(File);
 	if ( format == std::string("IEEE64BIG") ) {
 	  GaugeSimpleMunger<dobj, sobj> munge;
@@ -846,6 +901,13 @@ class IldgReader : public GridLimeReader {
    ////////////////////////////////////////////////////////////
    // Really really want to mandate a scidac checksum
    ////////////////////////////////////////////////////////////
+    if ( found_FieldNormMetaData ) { 
+      RealD nn = norm2(Umu);
+      GRID_FIELD_NORM_CHECK(FieldNormMetaData_,nn);
+      std::cout << GridLogMessage<<"FieldNormMetaData matches " << std::endl;
+    }  else { 
+      std::cout << GridLogWarning<<"FieldNormMetaData not found. " << std::endl;
+    }
    if ( found_scidacChecksum ) {
      FieldMetaData_.scidac_checksuma = stoull(scidacChecksum_.suma,0,16);
      FieldMetaData_.scidac_checksumb = stoull(scidacChecksum_.sumb,0,16);
--- a/Grid/parallelIO/MetaData.h
+++ b/Grid/parallelIO/MetaData.h
@@ -56,6 +56,10 @@ namespace Grid {
  ////////////////////////////////////////////////////////////////////////////////
  // header specification/interpretation
  ////////////////////////////////////////////////////////////////////////////////
+    class FieldNormMetaData : Serializable {
+    public:
+      GRID_SERIALIZABLE_CLASS_MEMBERS(FieldNormMetaData, double, norm2);
+    };
    class FieldMetaData : Serializable {
    public:

--- a/Grid/qcd/action/fermion/FermionOperator.h
+++ b/Grid/qcd/action/fermion/FermionOperator.h
@@ -136,6 +136,7 @@ namespace Grid {
      //////////////////////////////////////////////////////////////////////
      // Conserved currents, either contract at sink or insert sequentially.
      //////////////////////////////////////////////////////////////////////
+
      virtual void ContractConservedCurrent(PropagatorField &q_in_1,
                                            PropagatorField &q_in_2,
                                            PropagatorField &q_out,
@@ -148,6 +149,12 @@ namespace Grid {
                                       unsigned int tmin, 
                                       unsigned int tmax,
                                       ComplexField &lattice_cmplx)=0;
+
+      // Only reimplemented in Wilson5D 
+      // Default to just a zero correlation function
+      virtual void ContractJ5q(FermionField &q_in   ,ComplexField &J5q) { J5q=zero; };
+      virtual void ContractJ5q(PropagatorField &q_in,ComplexField &J5q) { J5q=zero; };
+
      ///////////////////////////////////////////////
      // Physical field import/export
      ///////////////////////////////////////////////
--- a/Grid/qcd/action/fermion/WilsonFermion5D.cc
+++ b/Grid/qcd/action/fermion/WilsonFermion5D.cc
@@ -939,6 +939,75 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const Fe
    merge(qSiteRev, qSiteVec); \
 }

+//          psi = chiralProjectPlus(Result_s[Ls/2-1]);
+//          psi+= chiralProjectMinus(Result_s[Ls/2]);
+//         PJ5q+=localInnerProduct(psi,psi);
+
+template<class vobj> 
+Lattice<vobj> spProj5p(const Lattice<vobj> & in)
+{
+  GridBase *grid=in._grid;
+  Gamma G5(Gamma::Algebra::Gamma5);
+  Lattice<vobj> ret(grid);
+  parallel_for(int ss=0;ss<grid->oSites();ss++){
+    ret._odata[ss] = in._odata[ss] + G5*in._odata[ss];
+  }
+  return ret;
+}
+template<class vobj> 
+Lattice<vobj> spProj5m(const Lattice<vobj> & in)
+{
+  Gamma G5(Gamma::Algebra::Gamma5);
+  GridBase *grid=in._grid;
+  Lattice<vobj> ret(grid);
+  parallel_for(int ss=0;ss<grid->oSites();ss++){
+    ret._odata[ss] = in._odata[ss] - G5*in._odata[ss];
+  }
+  return ret;
+}
+
+template <class Impl>
+void WilsonFermion5D<Impl>::ContractJ5q(FermionField &q_in,ComplexField &J5q)
+{
+  conformable(GaugeGrid(), J5q._grid);
+  conformable(q_in._grid, FermionGrid());
+
+  // 4d field
+  int Ls = this->Ls;
+  FermionField psi(GaugeGrid());
+  FermionField p_plus (GaugeGrid());
+  FermionField p_minus(GaugeGrid());
+  FermionField p(GaugeGrid());
+
+  ExtractSlice(p_plus , q_in, Ls/2   , 0);
+  ExtractSlice(p_minus, q_in, Ls/2-1 , 0);
+  p_plus = spProj5p(p_plus );
+  p_minus= spProj5m(p_minus);
+  p=p_plus+p_minus;
+  J5q = localInnerProduct(p,p);
+}
+
+template <class Impl>
+void WilsonFermion5D<Impl>::ContractJ5q(PropagatorField &q_in,ComplexField &J5q)
+{
+  conformable(GaugeGrid(), J5q._grid);
+  conformable(q_in._grid, FermionGrid());
+
+  // 4d field
+  int Ls = this->Ls;
+  PropagatorField psi(GaugeGrid());
+  PropagatorField p_plus (GaugeGrid());
+  PropagatorField p_minus(GaugeGrid());
+  PropagatorField p(GaugeGrid());
+
+  ExtractSlice(p_plus , q_in, Ls/2   , 0);
+  ExtractSlice(p_minus, q_in, Ls/2-1 , 0);
+  p_plus = spProj5p(p_plus );
+  p_minus= spProj5m(p_minus);
+  p=p_plus+p_minus;
+  J5q = localInnerProduct(p,p);
+}
+
 template <class Impl>
 void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
                                                     PropagatorField &q_in_2,
@@ -949,6 +1018,7 @@ void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
    conformable(q_in_1._grid, FermionGrid());
    conformable(q_in_1._grid, q_in_2._grid);
    conformable(_FourDimGrid, q_out._grid);
+
    PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid());
    unsigned int LLs = q_in_1._grid->_rdimensions[0];
    q_out = zero;
@@ -995,7 +1065,6 @@ void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
 }


-
 template <class Impl>
 void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in, 
                                                PropagatorField &q_out,
--- a/Grid/qcd/action/fermion/WilsonFermion5D.h
+++ b/Grid/qcd/action/fermion/WilsonFermion5D.h
@@ -230,6 +230,10 @@ namespace QCD {
                             unsigned int tmin, 
                             unsigned int tmax,
 			     ComplexField &lattice_cmplx);
+
+    void ContractJ5q(PropagatorField &q_in,ComplexField &J5q);
+    void ContractJ5q(FermionField &q_in,ComplexField &J5q);
+
  };

 }}
--- a/Grid/qcd/spin/Gamma.cc
+++ b/Grid/qcd/spin/Gamma.cc
@@ -11,6 +11,24 @@ const std::array<const Gamma, 4> Gamma::gmu = {{
  Gamma(Gamma::Algebra::GammaZ),
  Gamma(Gamma::Algebra::GammaT)}};

+const std::array<const Gamma, 16> Gamma::gall = {{
+  Gamma(Gamma::Algebra::Identity),
+  Gamma(Gamma::Algebra::Gamma5),
+  Gamma(Gamma::Algebra::GammaX),
+  Gamma(Gamma::Algebra::GammaY),
+  Gamma(Gamma::Algebra::GammaZ),
+  Gamma(Gamma::Algebra::GammaT),
+  Gamma(Gamma::Algebra::GammaXGamma5),
+  Gamma(Gamma::Algebra::GammaYGamma5),
+  Gamma(Gamma::Algebra::GammaZGamma5),
+  Gamma(Gamma::Algebra::GammaTGamma5),
+  Gamma(Gamma::Algebra::SigmaXT),      
+  Gamma(Gamma::Algebra::SigmaXY),      
+  Gamma(Gamma::Algebra::SigmaXZ),      
+  Gamma(Gamma::Algebra::SigmaYT),
+  Gamma(Gamma::Algebra::SigmaYZ),
+  Gamma(Gamma::Algebra::SigmaZT)}};
+
 const std::array<const char *, Gamma::nGamma> Gamma::name = {{
  "-Gamma5      ",
  "Gamma5       ",
--- a/Grid/qcd/spin/Gamma.h
+++ b/Grid/qcd/spin/Gamma.h
@@ -48,6 +48,7 @@ class Gamma {
    static const std::array<std::array<Algebra, nGamma>, nGamma> mul;
    static const std::array<Algebra, nGamma>                     adj;
    static const std::array<const Gamma, 4>                      gmu;
+    static const std::array<const Gamma, 16>                     gall;
    Algebra                                                      g;
  public:
    Gamma(Algebra initg): g(initg) {}  
--- a/Grid/qcd/spin/gamma-gen/gamma-gen.nb
+++ b/Grid/qcd/spin/gamma-gen/gamma-gen.nb
@@ -10,10 +10,10 @@
 NotebookFileLineBreakTest
 NotebookFileLineBreakTest
 NotebookDataPosition[       158,          7]
-NotebookDataLength[     75090,       1956]
-NotebookOptionsPosition[     69536,       1867]
-NotebookOutlinePosition[     69898,       1883]
-CellTagsIndexPosition[     69855,       1880]
+NotebookDataLength[     67118,       1714]
+NotebookOptionsPosition[     63485,       1652]
+NotebookOutlinePosition[     63842,       1668]
+CellTagsIndexPosition[     63799,       1665]
 WindowFrame->Normal*)

 (* Beginning of Notebook Content *)
@@ -76,234 +76,6 @@ Cell[BoxData["\<\"/Users/antonin/Development/Grid/lib/qcd/spin/gamma-gen\"\>"]\

 Cell[CellGroupData[{

-Cell[BoxData[
- RowBox[{"FactorInteger", "[", "3152", "]"}]], "Input",
- CellChangeTimes->{{3.7432347536316767`*^9, 3.7432347764739027`*^9}, {
-  3.743234833567358*^9, 
-  3.743234862146022*^9}},ExpressionUUID->"d1a0fd03-85e1-43af-ba80-\
-3ca4235675d8"],
-
-Cell[BoxData[
- RowBox[{"{", 
-  RowBox[{
-   RowBox[{"{", 
-    RowBox[{"2", ",", "4"}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{"197", ",", "1"}], "}"}]}], "}"}]], "Output",
- CellChangeTimes->{{3.743234836792224*^9, 
-  3.743234862493619*^9}},ExpressionUUID->"16d3f953-4b24-4ed2-ae62-\
-306dcab66ca7"]
-}, Open  ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"sol", "=", 
-  RowBox[{"Solve", "[", 
-   RowBox[{
-    RowBox[{
-     RowBox[{
-      SuperscriptBox["x", "2"], "+", 
-      SuperscriptBox["y", "2"], "+", 
-      SuperscriptBox["z", "2"]}], "\[Equal]", "2"}], ",", 
-    RowBox[{"{", 
-     RowBox[{"x", ",", "y", ",", "z"}], "}"}], ",", "Integers"}], 
-   "]"}]}]], "Input",
- CellChangeTimes->{{3.743235304127721*^9, 
-  3.7432353087929983`*^9}},ExpressionUUID->"f0fa2a5c-3d81-4d75-a447-\
-50c7ca3459ff"],
-
-Cell[BoxData[
- RowBox[{"{", 
-  RowBox[{
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", 
-      RowBox[{"-", "1"}]}], ",", 
-     RowBox[{"y", "\[Rule]", 
-      RowBox[{"-", "1"}]}], ",", 
-     RowBox[{"z", "\[Rule]", "0"}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", 
-      RowBox[{"-", "1"}]}], ",", 
-     RowBox[{"y", "\[Rule]", "0"}], ",", 
-     RowBox[{"z", "\[Rule]", 
-      RowBox[{"-", "1"}]}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", 
-      RowBox[{"-", "1"}]}], ",", 
-     RowBox[{"y", "\[Rule]", "0"}], ",", 
-     RowBox[{"z", "\[Rule]", "1"}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", 
-      RowBox[{"-", "1"}]}], ",", 
-     RowBox[{"y", "\[Rule]", "1"}], ",", 
-     RowBox[{"z", "\[Rule]", "0"}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", "0"}], ",", 
-     RowBox[{"y", "\[Rule]", 
-      RowBox[{"-", "1"}]}], ",", 
-     RowBox[{"z", "\[Rule]", 
-      RowBox[{"-", "1"}]}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", "0"}], ",", 
-     RowBox[{"y", "\[Rule]", 
-      RowBox[{"-", "1"}]}], ",", 
-     RowBox[{"z", "\[Rule]", "1"}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", "0"}], ",", 
-     RowBox[{"y", "\[Rule]", "1"}], ",", 
-     RowBox[{"z", "\[Rule]", 
-      RowBox[{"-", "1"}]}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", "0"}], ",", 
-     RowBox[{"y", "\[Rule]", "1"}], ",", 
-     RowBox[{"z", "\[Rule]", "1"}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", "1"}], ",", 
-     RowBox[{"y", "\[Rule]", 
-      RowBox[{"-", "1"}]}], ",", 
-     RowBox[{"z", "\[Rule]", "0"}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", "1"}], ",", 
-     RowBox[{"y", "\[Rule]", "0"}], ",", 
-     RowBox[{"z", "\[Rule]", 
-      RowBox[{"-", "1"}]}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", "1"}], ",", 
-     RowBox[{"y", "\[Rule]", "0"}], ",", 
-     RowBox[{"z", "\[Rule]", "1"}]}], "}"}], ",", 
-   RowBox[{"{", 
-    RowBox[{
-     RowBox[{"x", "\[Rule]", "1"}], ",", 
-     RowBox[{"y", "\[Rule]", "1"}], ",", 
-     RowBox[{"z", "\[Rule]", "0"}]}], "}"}]}], "}"}]], "Output",
- CellChangeTimes->{{3.743235305220907*^9, 
-  3.743235309139554*^9}},ExpressionUUID->"d9825c95-24bb-442a-8734-\
-4c0f47e99dfc"]
-}, Open  ]],
-
-Cell[BoxData[
- RowBox[{
-  RowBox[{"xmlElem", "[", "x_", "]"}], ":=", 
-  RowBox[{"Print", "[", 
-   RowBox[{"\"\<<elem>\>\"", "<>", 
-    RowBox[{"ToString", "[", 
-     RowBox[{"x", "[", 
-      RowBox[{"[", "1", "]"}], "]"}], "]"}], "<>", "\"\< \>\"", "<>", 
-    RowBox[{"ToString", "[", 
-     RowBox[{"x", "[", 
-      RowBox[{"[", "2", "]"}], "]"}], "]"}], "<>", "\"\< \>\"", "<>", 
-    RowBox[{"ToString", "[", 
-     RowBox[{"x", "[", 
-      RowBox[{"[", "3", "]"}], "]"}], "]"}], "<>", "\"\<</elem>\>\""}], 
-   "]"}]}]], "Input",
- CellChangeTimes->{{3.74323534002862*^9, 3.743235351000985*^9}, {
-  3.743235403233039*^9, 3.743235413488028*^9}, {3.743235473169856*^9, 
-  3.7432354747126904`*^9}},ExpressionUUID->"aea76313-c89e-45e8-b429-\
-3f454091666d"],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{
-  RowBox[{
-   RowBox[{"xmlElem", "[", 
-    RowBox[{
-     RowBox[{"{", 
-      RowBox[{"x", ",", "y", ",", "z"}], "}"}], "/.", "#"}], "]"}], "&"}], "/@",
-   "sol"}]], "Input",
- CellChangeTimes->{{3.743235415820318*^9, 
-  3.743235467025091*^9}},ExpressionUUID->"07da3998-8eab-40ba-8c0b-\
-ac6b130cb4fb"],
-
-Cell[CellGroupData[{
-
-Cell[BoxData["\<\"<elem>-1 -1 0</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476581676*^9},ExpressionUUID->"c577ba06-b67a-405a-9ff5-\
-2bf7dc898d03"],
-
-Cell[BoxData["\<\"<elem>-1 0 -1</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476588011*^9},ExpressionUUID->"d041aa36-0cea-457c-9d4b-\
-1fe9be66e2ab"],
-
-Cell[BoxData["\<\"<elem>-1 0 1</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476596887*^9},ExpressionUUID->"bf141b55-86b2-4430-a994-\
-5c03d5a19441"],
-
-Cell[BoxData["\<\"<elem>-1 1 0</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476605785*^9},ExpressionUUID->"4968a660-4ecf-4b66-9071-\
-8bd798c18d21"],
-
-Cell[BoxData["\<\"<elem>0 -1 -1</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476613523*^9},ExpressionUUID->"4e22d943-2680-416b-a1d7-\
-a16ca20b781f"],
-
-Cell[BoxData["\<\"<elem>0 -1 1</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.7432354766218576`*^9},ExpressionUUID->"6dd38385-08b3-4dd9-932f-\
-98a00c6db1b2"],
-
-Cell[BoxData["\<\"<elem>0 1 -1</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476629427*^9},ExpressionUUID->"ef3baad3-91d1-4735-9a22-\
-53495a624c15"],
-
-Cell[BoxData["\<\"<elem>0 1 1</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476638257*^9},ExpressionUUID->"413fbb68-5017-4272-a62a-\
-fa234e6daaea"],
-
-Cell[BoxData["\<\"<elem>1 -1 0</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476646203*^9},ExpressionUUID->"3a832a60-ae00-414b-a9ac-\
-f5e86e67e917"],
-
-Cell[BoxData["\<\"<elem>1 0 -1</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476653907*^9},ExpressionUUID->"bfc79ef6-f6c7-4f1e-88e8-\
-005ac314be9c"],
-
-Cell[BoxData["\<\"<elem>1 0 1</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.743235476662575*^9},ExpressionUUID->"0f892891-f885-489c-9925-\
-ddef4d698410"],
-
-Cell[BoxData["\<\"<elem>1 1 0</elem>\"\>"], "Print",
- CellChangeTimes->{
-  3.7432354766702337`*^9},ExpressionUUID->"2906f190-e673-4f33-9c34-\
-e8e56efe7a27"]
-}, Open  ]],
-
-Cell[BoxData[
- RowBox[{"{", 
-  RowBox[{
-  "Null", ",", "Null", ",", "Null", ",", "Null", ",", "Null", ",", "Null", 
-   ",", "Null", ",", "Null", ",", "Null", ",", "Null", ",", "Null", ",", 
-   "Null"}], "}"}]], "Output",
- CellChangeTimes->{
-  3.7432354246225967`*^9, {3.7432354674878073`*^9, 
-   3.743235476678007*^9}},ExpressionUUID->"500ca3c1-88d8-46e5-a1a1-\
-86a7878e5638"]
-}, Open  ]],
-
-Cell[CellGroupData[{
-
 Cell["Clifford algebra generation", "Section",
 CellChangeTimes->{{3.6942089434583883`*^9, 
  3.694208978559093*^9}},ExpressionUUID->"a5b064b3-3011-4922-8559-\
@@ -1048,8 +820,9 @@ generated by the Mathematica notebook gamma-gen/gamma-gen.nb\n\n#include \
       "\"\<    static const std::array<const char *, nGamma>                \
 name;\n    static const std::array<std::array<Algebra, nGamma>, nGamma> mul;\n\
    static const std::array<Algebra, nGamma>                     adj;\n    \
-static const std::array<const Gamma, 4>                      gmu;\n    \
-Algebra                                                      g;\n  public:\n  \
+static const std::array<const Gamma, 4>                      gmu;\n    static \
+const std::array<const Gamma, 16>                     gall;\n    Algebra      \
+                                                g;\n  public:\n    \
 Gamma(Algebra initg): g(initg) {}  \n};\n\n\>\""}]}], ";", 
     "\[IndentingNewLine]", 
     RowBox[{"out", " ", "=", 
@@ -1076,7 +849,8 @@ Algebra                                                      g;\n  public:\n  \
   3.694963343265525*^9}, {3.694964367519239*^9, 3.69496439461199*^9}, {
   3.694964462130747*^9, 3.6949644669959793`*^9}, 3.694964509762739*^9, {
   3.694964705045744*^9, 3.694964723148797*^9}, {3.694964992988984*^9, 
-   3.6949649968504257`*^9}},ExpressionUUID->"c7103bd6-b539-4495-b98c-\
+   3.6949649968504257`*^9}, {3.758291687176977*^9, 
+   3.758291694181189*^9}},ExpressionUUID->"c7103bd6-b539-4495-b98c-\
 d4d12ac6cad8"],

 Cell["Gamma enum generation:", "Text",
@@ -1745,8 +1519,17 @@ namespace QCD {\>\""}]}], ";", "\[IndentingNewLine]",
       "\"\<\n\nconst std::array<const Gamma, 4> Gamma::gmu = {{\n  \
 Gamma(Gamma::Algebra::GammaX),\n  Gamma(Gamma::Algebra::GammaY),\n  \
 Gamma(Gamma::Algebra::GammaZ),\n  Gamma(Gamma::Algebra::GammaT)}};\n\nconst \
-std::array<const char *, Gamma::nGamma> Gamma::name = {{\n\>\""}]}], ";", 
-     "\[IndentingNewLine]", 
+std::array<const Gamma, 16> Gamma::gall = {{\n  \
+Gamma(Gamma::Algebra::Identity),\n  Gamma(Gamma::Algebra::Gamma5),\n  \
+Gamma(Gamma::Algebra::GammaX),\n  Gamma(Gamma::Algebra::GammaY),\n  \
+Gamma(Gamma::Algebra::GammaZ),\n  Gamma(Gamma::Algebra::GammaT),\n  \
+Gamma(Gamma::Algebra::GammaXGamma5),\n  Gamma(Gamma::Algebra::GammaYGamma5),\n\
+  Gamma(Gamma::Algebra::GammaZGamma5),\n  \
+Gamma(Gamma::Algebra::GammaTGamma5),\n  Gamma(Gamma::Algebra::SigmaXT),      \
+\n  Gamma(Gamma::Algebra::SigmaXY),      \n  Gamma(Gamma::Algebra::SigmaXZ),  \
+    \n  Gamma(Gamma::Algebra::SigmaYT),\n  Gamma(Gamma::Algebra::SigmaYZ),\n  \
+Gamma(Gamma::Algebra::SigmaZT)}};\n\nconst std::array<const char *, \
+Gamma::nGamma> Gamma::name = {{\n\>\""}]}], ";", "\[IndentingNewLine]", 
     RowBox[{"Do", "[", "\[IndentingNewLine]", 
      RowBox[{
       RowBox[{"out", " ", "=", " ", 
@@ -1847,7 +1630,9 @@ Gamma::nGamma> Gamma::mul = {{\\n\>\""}]}], ";", "\[IndentingNewLine]",
   3.694963031525289*^9}, {3.694963065828494*^9, 3.694963098327538*^9}, {
   3.6949632020836153`*^9, 3.6949632715940027`*^9}, {3.694963440035037*^9, 
   3.6949634418966017`*^9}, {3.6949651447067547`*^9, 3.694965161228381*^9}, {
-   3.694967957845581*^9, 3.694967958364184*^9}}],
+   3.694967957845581*^9, 3.694967958364184*^9}, {3.758291673792514*^9, 
+   3.758291676983432*^9}},ExpressionUUID->"b1b309f8-a3a7-4081-a781-\
+c3845e3cd372"],

 Cell[BoxData[
 RowBox[{
@@ -1867,8 +1652,8 @@ Cell[BoxData[""], "Input",
 },
 WindowSize->{1246, 1005},
 WindowMargins->{{282, Automatic}, {Automatic, 14}},
-FrontEndVersion->"11.2 for Mac OS X x86 (32-bit, 64-bit Kernel) (September \
-10, 2017)",
+FrontEndVersion->"11.3 for Mac OS X x86 (32-bit, 64-bit Kernel) (March 5, \
+2018)",
 StyleDefinitions->"Default.nb"
 ]
 (* End of Notebook Content *)
@@ -1888,75 +1673,48 @@ Cell[1948, 43, 570, 11, 73, "Input",ExpressionUUID->"5c937a3e-adfd-4d7e-8fde-afb
 Cell[2521, 56, 1172, 17, 34, "Output",ExpressionUUID->"72817ba6-2f6a-4a4d-8212-6f0970f49e7c"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[3730, 78, 248, 5, 30, "Input",ExpressionUUID->"d1a0fd03-85e1-43af-ba80-3ca4235675d8"],
-Cell[3981, 85, 299, 9, 34, "Output",ExpressionUUID->"16d3f953-4b24-4ed2-ae62-306dcab66ca7"]
+Cell[3730, 78, 174, 3, 67, "Section",ExpressionUUID->"a5b064b3-3011-4922-8559-ead857cad102"],
+Cell[3907, 83, 535, 16, 52, "Input",ExpressionUUID->"aa28f02b-31e1-4df2-9b5d-482177464b59"],
+Cell[4445, 101, 250, 4, 35, "Text",ExpressionUUID->"c8896b88-f1db-4ce4-b7a6-0c9838bdb8f1"],
+Cell[4698, 107, 5511, 169, 425, "Input",ExpressionUUID->"52a96ff6-047e-4043-86d0-e303866e5f8e"],
+Cell[CellGroupData[{
+Cell[10234, 280, 2183, 58, 135, "Input",ExpressionUUID->"8b0f4955-2c3f-418c-9226-9be8f87621e8"],
+Cell[12420, 340, 1027, 27, 56, "Output",ExpressionUUID->"edd0619f-6f12-4070-a1d2-6b547877fadc"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[4317, 99, 469, 14, 33, "Input",ExpressionUUID->"f0fa2a5c-3d81-4d75-a447-50c7ca3459ff"],
-Cell[4789, 115, 2423, 77, 56, "Output",ExpressionUUID->"d9825c95-24bb-442a-8734-4c0f47e99dfc"]
+Cell[13484, 372, 1543, 46, 114, "Input",ExpressionUUID->"fb45123c-c610-4075-99b0-7cd71c728ae7"],
+Cell[15030, 420, 1311, 32, 87, "Output",ExpressionUUID->"2ae14565-b412-4dc0-9dce-bd6c1ba5ef27"]
 }, Open  ]],
-Cell[7227, 195, 751, 18, 30, "Input",ExpressionUUID->"aea76313-c89e-45e8-b429-3f454091666d"],
+Cell[16356, 455, 179, 3, 35, "Text",ExpressionUUID->"af247231-a58d-417b-987a-26908dafffdb"],
+Cell[16538, 460, 2175, 65, 94, "Input",ExpressionUUID->"7c44cadd-e488-4f51-87d8-c64eef11f40c"],
+Cell[18716, 527, 193, 3, 35, "Text",ExpressionUUID->"856f1746-1107-4509-a5ce-ac9c7f56cdb1"],
 Cell[CellGroupData[{
-Cell[8003, 217, 323, 10, 30, "Input",ExpressionUUID->"07da3998-8eab-40ba-8c0b-ac6b130cb4fb"],
-Cell[CellGroupData[{
-Cell[8351, 231, 156, 3, 24, "Print",ExpressionUUID->"c577ba06-b67a-405a-9ff5-2bf7dc898d03"],
-Cell[8510, 236, 156, 3, 24, "Print",ExpressionUUID->"d041aa36-0cea-457c-9d4b-1fe9be66e2ab"],
-Cell[8669, 241, 155, 3, 24, "Print",ExpressionUUID->"bf141b55-86b2-4430-a994-5c03d5a19441"],
-Cell[8827, 246, 155, 3, 24, "Print",ExpressionUUID->"4968a660-4ecf-4b66-9071-8bd798c18d21"],
-Cell[8985, 251, 156, 3, 24, "Print",ExpressionUUID->"4e22d943-2680-416b-a1d7-a16ca20b781f"],
-Cell[9144, 256, 157, 3, 24, "Print",ExpressionUUID->"6dd38385-08b3-4dd9-932f-98a00c6db1b2"],
-Cell[9304, 261, 155, 3, 24, "Print",ExpressionUUID->"ef3baad3-91d1-4735-9a22-53495a624c15"],
-Cell[9462, 266, 154, 3, 24, "Print",ExpressionUUID->"413fbb68-5017-4272-a62a-fa234e6daaea"],
-Cell[9619, 271, 155, 3, 24, "Print",ExpressionUUID->"3a832a60-ae00-414b-a9ac-f5e86e67e917"],
-Cell[9777, 276, 155, 3, 24, "Print",ExpressionUUID->"bfc79ef6-f6c7-4f1e-88e8-005ac314be9c"],
-Cell[9935, 281, 154, 3, 24, "Print",ExpressionUUID->"0f892891-f885-489c-9925-ddef4d698410"],
-Cell[10092, 286, 156, 3, 24, "Print",ExpressionUUID->"2906f190-e673-4f33-9c34-e8e56efe7a27"]
-}, Open  ]],
-Cell[10263, 292, 376, 9, 34, "Output",ExpressionUUID->"500ca3c1-88d8-46e5-a1a1-86a7878e5638"]
+Cell[18934, 534, 536, 16, 30, "Input",ExpressionUUID->"8674484a-8543-434f-b177-3b27f9353212"],
+Cell[19473, 552, 1705, 35, 87, "Output",ExpressionUUID->"c3b3f84d-91f6-41af-af6b-a394ca020511"]
 }, Open  ]],
+Cell[21193, 590, 170, 3, 35, "Text",ExpressionUUID->"518a3040-54b1-4d43-8947-5c7d12efa94d"],
 Cell[CellGroupData[{
-Cell[10676, 306, 174, 3, 67, "Section",ExpressionUUID->"a5b064b3-3011-4922-8559-ead857cad102"],
-Cell[10853, 311, 535, 16, 52, "Input",ExpressionUUID->"aa28f02b-31e1-4df2-9b5d-482177464b59"],
-Cell[11391, 329, 250, 4, 35, "Text",ExpressionUUID->"c8896b88-f1db-4ce4-b7a6-0c9838bdb8f1"],
-Cell[11644, 335, 5511, 169, 425, "Input",ExpressionUUID->"52a96ff6-047e-4043-86d0-e303866e5f8e"],
-Cell[CellGroupData[{
-Cell[17180, 508, 2183, 58, 135, "Input",ExpressionUUID->"8b0f4955-2c3f-418c-9226-9be8f87621e8"],
-Cell[19366, 568, 1027, 27, 67, "Output",ExpressionUUID->"edd0619f-6f12-4070-a1d2-6b547877fadc"]
-}, Open  ]],
-Cell[CellGroupData[{
-Cell[20430, 600, 1543, 46, 114, "Input",ExpressionUUID->"fb45123c-c610-4075-99b0-7cd71c728ae7"],
-Cell[21976, 648, 1311, 32, 98, "Output",ExpressionUUID->"2ae14565-b412-4dc0-9dce-bd6c1ba5ef27"]
-}, Open  ]],
-Cell[23302, 683, 179, 3, 35, "Text",ExpressionUUID->"af247231-a58d-417b-987a-26908dafffdb"],
-Cell[23484, 688, 2175, 65, 94, "Input",ExpressionUUID->"7c44cadd-e488-4f51-87d8-c64eef11f40c"],
-Cell[25662, 755, 193, 3, 35, "Text",ExpressionUUID->"856f1746-1107-4509-a5ce-ac9c7f56cdb1"],
-Cell[CellGroupData[{
-Cell[25880, 762, 536, 16, 30, "Input",ExpressionUUID->"8674484a-8543-434f-b177-3b27f9353212"],
-Cell[26419, 780, 1705, 35, 87, "Output",ExpressionUUID->"c3b3f84d-91f6-41af-af6b-a394ca020511"]
-}, Open  ]],
-Cell[28139, 818, 170, 3, 35, "Text",ExpressionUUID->"518a3040-54b1-4d43-8947-5c7d12efa94d"],
-Cell[CellGroupData[{
-Cell[28334, 825, 536, 14, 30, "Input",ExpressionUUID->"61a2e974-2b39-4a07-8043-2dfd39a70569"],
-Cell[28873, 841, 6754, 167, 303, "Output",ExpressionUUID->"73480ac0-3043-4077-80cc-b952a94c822a"]
+Cell[21388, 597, 536, 14, 30, "Input",ExpressionUUID->"61a2e974-2b39-4a07-8043-2dfd39a70569"],
+Cell[21927, 613, 6754, 167, 303, "Output",ExpressionUUID->"73480ac0-3043-4077-80cc-b952a94c822a"]
 }, Open  ]]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[35676, 1014, 226, 4, 67, "Section",ExpressionUUID->"4e833cd6-9f0e-4aa3-a873-3d579e874720"],
-Cell[35905, 1020, 188, 4, 44, "Text",ExpressionUUID->"6d27fc04-3a60-4e03-8df7-3dd3aeee35b4"],
-Cell[36096, 1026, 2980, 53, 703, "Input",ExpressionUUID->"c7103bd6-b539-4495-b98c-d4d12ac6cad8"],
-Cell[39079, 1081, 221, 4, 44, "Text",ExpressionUUID->"0625593d-290f-4a39-9d80-8e2c6fdbc94e"],
-Cell[39303, 1087, 4936, 150, 682, "Input",ExpressionUUID->"1ad4904c-352f-4b1d-a7c7-91e1b0549409"],
-Cell[44242, 1239, 2645, 56, 199, "Input",ExpressionUUID->"0221674f-9b63-4662-91bc-ccc8c6ae9589"],
-Cell[46890, 1297, 209, 4, 44, "Text",ExpressionUUID->"d2d2257a-487b-416f-bc40-abd4482225f7"],
-Cell[47102, 1303, 15306, 397, 2131, "Input",ExpressionUUID->"daea68a9-c9e8-46ab-9bc8-5186e2cf477c"],
-Cell[62411, 1702, 137, 2, 44, "Text",ExpressionUUID->"76ba9d5a-7ee3-4888-be7e-6377003275e8"],
-Cell[62551, 1706, 521, 12, 30, "Input",ExpressionUUID->"4ec61f4c-3fd3-49ea-b5ef-6f7f04a16b34"]
+Cell[28730, 786, 226, 4, 67, "Section",ExpressionUUID->"4e833cd6-9f0e-4aa3-a873-3d579e874720"],
+Cell[28959, 792, 188, 4, 44, "Text",ExpressionUUID->"6d27fc04-3a60-4e03-8df7-3dd3aeee35b4"],
+Cell[29150, 798, 3104, 55, 724, "Input",ExpressionUUID->"c7103bd6-b539-4495-b98c-d4d12ac6cad8"],
+Cell[32257, 855, 221, 4, 44, "Text",ExpressionUUID->"0625593d-290f-4a39-9d80-8e2c6fdbc94e"],
+Cell[32481, 861, 4936, 150, 682, "Input",ExpressionUUID->"1ad4904c-352f-4b1d-a7c7-91e1b0549409"],
+Cell[37420, 1013, 2645, 56, 199, "Input",ExpressionUUID->"0221674f-9b63-4662-91bc-ccc8c6ae9589"],
+Cell[40068, 1071, 209, 4, 44, "Text",ExpressionUUID->"d2d2257a-487b-416f-bc40-abd4482225f7"],
+Cell[40280, 1077, 15306, 397, 2131, "Input",ExpressionUUID->"daea68a9-c9e8-46ab-9bc8-5186e2cf477c"],
+Cell[55589, 1476, 137, 2, 44, "Text",ExpressionUUID->"76ba9d5a-7ee3-4888-be7e-6377003275e8"],
+Cell[55729, 1480, 521, 12, 30, "Input",ExpressionUUID->"4ec61f4c-3fd3-49ea-b5ef-6f7f04a16b34"]
 }, Open  ]],
 Cell[CellGroupData[{
-Cell[63109, 1723, 167, 2, 67, "Section",ExpressionUUID->"a4458b3a-09b5-4e36-a1fc-781d6702b2dc"],
-Cell[63279, 1727, 5693, 122, 829, "Input",ExpressionUUID->"b1b309f8-a3a7-4081-a781-c3845e3cd372"],
-Cell[68975, 1851, 448, 10, 30, "Input",ExpressionUUID->"cba42949-b0f2-42ce-aebd-ffadfd83ef88"],
-Cell[69426, 1863, 94, 1, 30, "Input",ExpressionUUID->"6175b72c-af9f-43c2-b4ca-bd84c48a456d"]
+Cell[56287, 1497, 167, 2, 67, "Section",ExpressionUUID->"a4458b3a-09b5-4e36-a1fc-781d6702b2dc"],
+Cell[56457, 1501, 6464, 133, 1207, "Input",ExpressionUUID->"b1b309f8-a3a7-4081-a781-c3845e3cd372"],
+Cell[62924, 1636, 448, 10, 30, "Input",ExpressionUUID->"cba42949-b0f2-42ce-aebd-ffadfd83ef88"],
+Cell[63375, 1648, 94, 1, 30, "Input",ExpressionUUID->"6175b72c-af9f-43c2-b4ca-bd84c48a456d"]
 }, Open  ]]
 }
 ]
--- a/Grid/qcd/utils/CovariantSmearing.h
+++ b/Grid/qcd/utils/CovariantSmearing.h
@@ -0,0 +1,87 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid
+
+Source file: ./lib/qcd/action/scalar/CovariantLaplacian.h
+
+Copyright (C) 2016
+
+Author: Azusa Yamaguchi
+
+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
+*************************************************************************************/
+#pragma once
+
+namespace Grid {
+namespace QCD {
+
+template <class Gimpl> class CovariantSmearing : public Gimpl 
+{
+public:
+  INHERIT_GIMPL_TYPES(Gimpl);
+
+  typedef typename Gimpl::GaugeLinkField GaugeMat;
+  typedef typename Gimpl::GaugeField GaugeLorentz;
+
+  template<typename T>
+  static void GaussianSmear(const std::vector<LatticeColourMatrix>& U, 
+			    T& chi, 
+			    const Real& width, int Iterations, int orthog)
+  {
+    GridBase *grid = chi._grid;
+    T psi(grid);
+
+    ////////////////////////////////////////////////////////////////////////////////////
+    // Follow Chroma conventions for width to keep compatibility with previous data
+    // Free field iterates 
+    //   chi = (1 - w^2/4N p^2)^N chi
+    //
+    //       ~ (e^(-w^2/4N p^2)^N chi
+    //       ~ (e^(-w^2/4 p^2) chi
+    //       ~ (e^(-w'^2/2 p^2) chi          [ w' = w/sqrt(2) ]
+    //
+    // Which in coordinate space is proportional to
+    //
+    //   e^(-x^2/w^2) = e^(-x^2/2w'^2) 
+    //
+    // The 4 is a bit unconventional from Gaussian width perspective, but... it's Chroma convention.
+    // 2nd derivative approx d^2/dx^2  =  x+mu + x-mu - 2x
+    //
+    // d^2/dx^2 = - p^2
+    //
+    // chi = ( 1 + w^2/4N d^2/dx^2 )^N chi
+    //
+    ////////////////////////////////////////////////////////////////////////////////////
+    Real coeff = (width*width) / Real(4*Iterations);
+  
+    int dims = Nd;
+    if( orthog < Nd ) dims=Nd-1;
+
+    for(int n = 0; n < Iterations; ++n) {
+      psi = (-2.0*dims)*chi;
+      for(int mu=0;mu<Nd;mu++) {
+	if ( mu != orthog ) { 
+	  psi = psi + Gimpl::CovShiftForward(U[mu],mu,chi);    
+	  psi = psi + Gimpl::CovShiftBackward(U[mu],mu,chi);    
+	}
+      }
+      chi = chi + coeff*psi;
+    }
+  }
+};
+}}
--- a/Grid/qcd/utils/GaugeFix.h
+++ b/Grid/qcd/utils/GaugeFix.h
@@ -53,21 +53,32 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
  }  
  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) {
    GridBase *grid = Umu._grid;
+    GaugeMat xform(grid);
+    SteepestDescentGaugeFix(Umu,xform,alpha,maxiter,Omega_tol,Phi_tol,Fourier);
+  }
+  static void SteepestDescentGaugeFix(GaugeLorentz &Umu,GaugeMat &xform,Real & alpha,int maxiter,Real Omega_tol, Real Phi_tol,bool Fourier=false) {
+
+    GridBase *grid = Umu._grid;

    Real org_plaq      =WilsonLoops<Gimpl>::avgPlaquette(Umu);
    Real org_link_trace=WilsonLoops<Gimpl>::linkTrace(Umu); 
    Real old_trace = org_link_trace;
    Real trG;
    
+    xform=1.0;
+
    std::vector<GaugeMat> U(Nd,grid);
+
    GaugeMat dmuAmu(grid);

    for(int i=0;i<maxiter;i++){
+
      for(int mu=0;mu<Nd;mu++) U[mu]= PeekIndex<LorentzIndex>(Umu,mu);
+
      if ( Fourier==false ) { 
-	trG = SteepestDescentStep(U,alpha,dmuAmu);
+	trG = SteepestDescentStep(U,xform,alpha,dmuAmu);
      } else { 
-	trG = FourierAccelSteepestDescentStep(U,alpha,dmuAmu);
+	trG = FourierAccelSteepestDescentStep(U,xform,alpha,dmuAmu);
      }
      for(int mu=0;mu<Nd;mu++) PokeIndex<LorentzIndex>(Umu,U[mu],mu);
      // Monitor progress and convergence test 
@@ -84,7 +95,6 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
 	Real Phi  = 1.0 - old_trace / link_trace ;
 	Real Omega= 1.0 - trG;

-
 	std::cout << GridLogMessage << " Iteration "<<i<< " Phi= "<<Phi<< " Omega= " << Omega<< " trG " << trG <<std::endl;
 	if ( (Omega < Omega_tol) && ( ::fabs(Phi) < Phi_tol) ) {
 	  std::cout << GridLogMessage << "Converged ! "<<std::endl;
@@ -96,7 +106,7 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
      }
    }
  };
-  static Real SteepestDescentStep(std::vector<GaugeMat> &U,Real & alpha, GaugeMat & dmuAmu) {
+  static Real SteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu) {
    GridBase *grid = U[0]._grid;

    std::vector<GaugeMat> A(Nd,grid);
@@ -109,12 +119,13 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
    Real vol = grid->gSites();
    Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;

+    xform = g*xform ;
    SU<Nc>::GaugeTransform(U,g);

    return trG;
  }

-  static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,Real & alpha, GaugeMat & dmuAmu) {
+  static Real FourierAccelSteepestDescentStep(std::vector<GaugeMat> &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu) {

    GridBase *grid = U[0]._grid;

@@ -153,13 +164,6 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {
    Complex psqMax(16.0);
    Fp =  psqMax*one/psq;

-    /*
-    static int once;
-    if ( once == 0 ) { 
-      std::cout << " Fp " << Fp <<std::endl;
-      once ++;
-      }*/
-
    pokeSite(TComplex(1.0),Fp,coor);

    dmuAmu_p  = dmuAmu_p * Fp; 
@@ -173,6 +177,7 @@ class FourierAcceleratedGaugeFixer  : public Gimpl {

    Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc;

+    xform = g*xform ;
    SU<Nc>::GaugeTransform(U,g);

    return trG;
--- a/Grid/qcd/utils/SUn.h
+++ b/Grid/qcd/utils/SUn.h
@@ -676,9 +676,17 @@ class SU {
    }
  }
 /*
- add GaugeTrans
+ * Fundamental rep gauge xform
+ */
+  template<typename Fundamental,typename GaugeMat>
+  static void GaugeTransformFundamental( Fundamental &ferm, GaugeMat &g){
+    GridBase *grid = ferm._grid;
+    conformable(grid,g._grid);
+    ferm = g*ferm;
+  }
+/*
+ * Adjoint rep gauge xform
 */
-
  template<typename GaugeField,typename GaugeMat>
  static void GaugeTransform( GaugeField &Umu, GaugeMat &g){
    GridBase *grid = Umu._grid;
--- a/Grid/util/Init.cc
+++ b/Grid/util/Init.cc
@@ -289,6 +289,11 @@ void Grid_init(int *argc,char ***argv)
    std::cout << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the"<<std::endl;
    std::cout << "GNU General Public License for more details."<<std::endl;
    printHash();
+  #ifdef GRID_BUILD_REF
+  #define _GRID_BUILD_STR(x) #x
+  #define GRID_BUILD_STR(x) _GRID_BUILD_STR(x)
+    std::cout << "Build " << GRID_BUILD_STR(GRID_BUILD_REF) << std::endl;
+  #endif
    std::cout << std::endl;
  }

--- a/Hadrons/A2AMatrix.hpp
+++ b/Hadrons/A2AMatrix.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/A2AMatrix.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
--- a/Hadrons/A2AVectors.hpp
+++ b/Hadrons/A2AVectors.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/A2AVectors.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: fionnoh <fionnoh@gmail.com>
--- a/Hadrons/Application.cc
+++ b/Hadrons/Application.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Application.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -48,6 +48,9 @@ Application::Application(void)
 {
    initLogger();
    auto dim = GridDefaultLatt(), mpi = GridDefaultMpi(), loc(dim);
+
+    if (dim.size())
+    {
        locVol_ = 1;
        for (unsigned int d = 0; d < dim.size(); ++d)
        {
@@ -71,6 +74,7 @@ Application::Application(void)
        LOG(Message) << "Schur decomposition     : " << MACOUTS(HADRONS_DEFAULT_SCHUR) << std::endl;
        LOG(Message) << std::endl;
    }
+}

 Application::Application(const Application::GlobalPar &par)
 : Application()
--- a/Hadrons/Application.hpp
+++ b/Hadrons/Application.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Application.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Archive/Modules/ScalarVP.cc
+++ b/Hadrons/Archive/Modules/ScalarVP.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MScalar/ScalarVP.cc
+Source file: Hadrons/Archive/Modules/ScalarVP.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
--- a/Hadrons/Archive/Modules/ScalarVP.hpp
+++ b/Hadrons/Archive/Modules/ScalarVP.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MScalar/ScalarVP.hpp
+Source file: Hadrons/Archive/Modules/ScalarVP.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
--- a/Hadrons/Modules/MUtilities/TestSeqConserved.cc
+++ b/Hadrons/Modules/MUtilities/TestSeqConserved.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MUtilities/TestSeqConserved.cc
+Source file: Hadrons/Archive/Modules/TestSeqConserved.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MUtilities/TestSeqConserved.hpp
+++ b/Hadrons/Modules/MUtilities/TestSeqConserved.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MUtilities/TestSeqConserved.hpp
+Source file: Hadrons/Archive/Modules/TestSeqConserved.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MUtilities/TestSeqGamma.cc
+++ b/Hadrons/Modules/MUtilities/TestSeqGamma.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MUtilities/TestSeqGamma.cc
+Source file: Hadrons/Archive/Modules/TestSeqGamma.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MUtilities/TestSeqGamma.hpp
+++ b/Hadrons/Modules/MUtilities/TestSeqGamma.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MUtilities/TestSeqGamma.hpp
+Source file: Hadrons/Archive/Modules/TestSeqGamma.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Archive/Modules/VPCounterTerms.cc
+++ b/Hadrons/Archive/Modules/VPCounterTerms.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MScalar/VPCounterTerms.cc
+Source file: Hadrons/Archive/Modules/VPCounterTerms.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
--- a/Hadrons/Archive/Modules/VPCounterTerms.hpp
+++ b/Hadrons/Archive/Modules/VPCounterTerms.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MScalar/VPCounterTerms.hpp
+Source file: Hadrons/Archive/Modules/VPCounterTerms.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
--- a/Hadrons/Modules/MContraction/WardIdentity.cc
+++ b/Hadrons/Modules/MContraction/WardIdentity.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WardIdentity.cc
+Source file: Hadrons/Archive/Modules/WardIdentity.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MContraction/WardIdentity.hpp
+++ b/Hadrons/Modules/MContraction/WardIdentity.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WardIdentity.hpp
+Source file: Hadrons/Archive/Modules/WardIdentity.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/WeakHamiltonian.hpp
+++ b/Hadrons/Modules/MContraction/WeakHamiltonian.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WeakHamiltonian.hpp
+Source file: Hadrons/Archive/Modules/WeakHamiltonian.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
+++ b/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
+Source file: Hadrons/Archive/Modules/WeakHamiltonianEye.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp
+++ b/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp
+Source file: Hadrons/Archive/Modules/WeakHamiltonianEye.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
+++ b/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
+Source file: Hadrons/Archive/Modules/WeakHamiltonianNonEye.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp
+++ b/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp
+Source file: Hadrons/Archive/Modules/WeakHamiltonianNonEye.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
+++ b/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
+Source file: Hadrons/Archive/Modules/WeakNeutral4ptDisc.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp
+++ b/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp
+Source file: Hadrons/Archive/Modules/WeakNeutral4ptDisc.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/DilutedNoise.hpp
+++ b/Hadrons/DilutedNoise.hpp
@@ -4,10 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/DilutedNoise.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <Vera.Guelpers@ed.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
--- a/Hadrons/DiskVector.hpp
+++ b/Hadrons/DiskVector.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/DiskVector.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/EigenPack.hpp
+++ b/Hadrons/EigenPack.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/EigenPack.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -309,6 +309,8 @@ class CoarseEigenPack: public EigenPack<FineF, FineFIo>
 {
 public:
    typedef CoarseF   CoarseField;
+    typedef CoarseFIo CoarseFieldIo;
+public:      
    std::vector<CoarseF> evecCoarse;
    std::vector<RealD>   evalCoarse;
 public:
--- a/Hadrons/Environment.cc
+++ b/Hadrons/Environment.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Environment.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -45,13 +45,11 @@ Environment::Environment(void)
 {
    dim_ = GridDefaultLatt();
    nd_  = dim_.size();
-    createGrid<vComplex>(1);
    vol_ = 1.;
    for (auto d: dim_)
    {
        vol_ *= d;
    }
-    rng4d_.reset(new GridParallelRNG(getGrid()));
 }

 // grids ///////////////////////////////////////////////////////////////////////
@@ -76,8 +74,13 @@ double Environment::getVolume(void) const
 }

 // random number generator /////////////////////////////////////////////////////
-GridParallelRNG * Environment::get4dRng(void) const
+GridParallelRNG * Environment::get4dRng(void)
 {
+    if (rng4d_ == nullptr)
+    {
+        rng4d_.reset(new GridParallelRNG(getGrid()));
+    }
+
    return rng4d_.get();
 }

--- a/Hadrons/Environment.hpp
+++ b/Hadrons/Environment.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Environment.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -113,7 +113,7 @@ public:
    unsigned int            getNd(void) const;
    double                  getVolume(void) const;
    // random number generator
-    GridParallelRNG *       get4dRng(void) const;
+    GridParallelRNG *       get4dRng(void);
    // general memory management
    void                    addObject(const std::string name,
                                      const int moduleAddress = -1);
@@ -182,7 +182,7 @@ private:
    std::map<CoarseGridKey, GridPt>     gridCoarse5d_;
    unsigned int                        nd_;
    // random number generator
-    RngPt                               rng4d_;
+    RngPt                               rng4d_{nullptr};
    // object store
    std::vector<ObjInfo>                object_;
    std::map<std::string, unsigned int> objectAddress_;
--- a/Hadrons/Exceptions.cc
+++ b/Hadrons/Exceptions.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Exceptions.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Exceptions.hpp
+++ b/Hadrons/Exceptions.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Exceptions.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Factory.hpp
+++ b/Hadrons/Factory.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Factory.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/GeneticScheduler.hpp
+++ b/Hadrons/GeneticScheduler.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/GeneticScheduler.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Global.cc
+++ b/Hadrons/Global.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Global.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Global.hpp
+++ b/Hadrons/Global.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Global.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
@@ -44,6 +44,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #define DEFAULT_ASCII_PREC 16
 #endif

+#define ARG(...) __VA_ARGS__
+
 /* the 'using Grid::operator<<;' statement prevents a very nasty compilation
 * error with GCC 5 (clang & GCC 6 compile fine without it).
 */
@@ -101,15 +103,16 @@ BEGIN_HADRONS_NAMESPACE
 typedef typename Impl::Field                         ScalarField##suffix;\
 typedef typename Impl::PropagatorField               PropagatorField##suffix;\
 typedef typename Impl::SitePropagator::scalar_object SitePropagator##suffix;\
-typedef std::vector<SitePropagator##suffix>          SlicedPropagator##suffix;
+typedef typename Impl::ComplexField                  ComplexField##suffix;\
+typedef std::vector<SitePropagator##suffix>          SlicedPropagator##suffix;\
+typedef std::vector<typename ComplexField##suffix::vector_object::scalar_object> SlicedComplex##suffix;

 #define FERM_TYPE_ALIASES(FImpl, suffix)\
 BASIC_TYPE_ALIASES(FImpl, suffix);\
 typedef FermionOperator<FImpl>            FMat##suffix;\
 typedef typename FImpl::FermionField      FermionField##suffix;\
 typedef typename FImpl::GaugeField        GaugeField##suffix;\
-typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;\
-typedef typename FImpl::ComplexField      ComplexField##suffix;
+typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;

 #define GAUGE_TYPE_ALIASES(GImpl, suffix)\
 typedef typename GImpl::GaugeField GaugeField##suffix;
@@ -263,6 +266,15 @@ void tokenReplace(std::string &str, const std::string token,
    }
 }

+// generic correlator class
+template <typename Metadata, typename Scalar = Complex>
+struct Correlator: Serializable
+{
+    GRID_SERIALIZABLE_CLASS_MEMBERS(ARG(Correlator<Metadata, Scalar>),
+                                    Metadata,             info,
+                                    std::vector<Complex>, corr);
+};
+
 END_HADRONS_NAMESPACE

 #include <Hadrons/Exceptions.hpp>
--- a/Hadrons/Graph.hpp
+++ b/Hadrons/Graph.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Graph.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Module.cc
+++ b/Hadrons/Module.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Module.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Module.hpp
+++ b/Hadrons/Module.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Module.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -65,7 +65,6 @@ static ns##mod##ModuleRegistrar ns##mod##ModuleRegistrarInstance;
 extern template class base;\
 MODULE_REGISTER(mod, ARG(base), ns);

-#define ARG(...) __VA_ARGS__
 #define HADRONS_MACRO_REDIRECT_12(arg1, arg2, macro, ...) macro
 #define HADRONS_MACRO_REDIRECT_23(arg1, arg2, arg3, macro, ...) macro

@@ -78,6 +77,15 @@ env().template getGrid<typename latticeType::vector_type>(Ls)
 #define envGetGrid(...)\
 HADRONS_MACRO_REDIRECT_12(__VA_ARGS__, envGetGrid5, envGetGrid4)(__VA_ARGS__)

+#define envGetCoarseGrid4(latticeType, blockSize)\
+env().template getCoarseGrid<typename latticeType::vector_type>(blockSize)
+
+#define envGetCoarseGrid5(latticeType, blockSize, Ls)\
+env().template getCoarseGrid<typename latticeType::vector_type>(blockSize, Ls)
+
+#define envGetCoarseGrid(...)\
+HADRONS_MACRO_REDIRECT_23(__VA_ARGS__, envGetCoarseGrid5, envGetCoarseGrid4)(__VA_ARGS__)
+
 #define envGetRbGrid4(latticeType)\
 env().template getRbGrid<typename latticeType::vector_type>()

--- a/Hadrons/ModuleFactory.hpp
+++ b/Hadrons/ModuleFactory.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/ModuleFactory.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules.hpp
+++ b/Hadrons/Modules.hpp
@@ -1,15 +1,13 @@
+#include <Hadrons/Modules/MContraction/WeakEye3pt.hpp>
 #include <Hadrons/Modules/MContraction/Baryon.hpp>
 #include <Hadrons/Modules/MContraction/A2AAslashField.hpp>
+#include <Hadrons/Modules/MContraction/A2ALoop.hpp>
 #include <Hadrons/Modules/MContraction/A2AMesonField.hpp>
 #include <Hadrons/Modules/MContraction/Meson.hpp>
-#include <Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
-#include <Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
 #include <Hadrons/Modules/MContraction/DiscLoop.hpp>
-#include <Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
 #include <Hadrons/Modules/MContraction/Gamma3pt.hpp>
-#include <Hadrons/Modules/MContraction/WardIdentity.hpp>
-#include <Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
 #include <Hadrons/Modules/MContraction/WeakMesonDecayKl2.hpp>
+#include <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
 #include <Hadrons/Modules/MFermion/FreeProp.hpp>
 #include <Hadrons/Modules/MFermion/GaugeProp.hpp>
 #include <Hadrons/Modules/MSource/SeqGamma.hpp>
@@ -23,28 +21,23 @@
 #include <Hadrons/Modules/MSink/Point.hpp>
 #include <Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp>
 #include <Hadrons/Modules/MSolver/LocalCoherenceLanczos.hpp>
+#include <Hadrons/Modules/MSolver/A2AAslashVectors.hpp>
 #include <Hadrons/Modules/MSolver/Guesser.hpp>
 #include <Hadrons/Modules/MSolver/RBPrecCG.hpp>
 #include <Hadrons/Modules/MSolver/A2AVectors.hpp>
-#include <Hadrons/Modules/MSolver/A2AAslashVectors.hpp>
 #include <Hadrons/Modules/MGauge/UnitEm.hpp>
 #include <Hadrons/Modules/MGauge/StoutSmearing.hpp>
 #include <Hadrons/Modules/MGauge/Unit.hpp>
-#include <Hadrons/Modules/MGauge/Electrify.hpp>
 #include <Hadrons/Modules/MGauge/Random.hpp>
 #include <Hadrons/Modules/MGauge/GaugeFix.hpp>
 #include <Hadrons/Modules/MGauge/FundtoHirep.hpp>
 #include <Hadrons/Modules/MGauge/StochEm.hpp>
+#include <Hadrons/Modules/MGauge/Electrify.hpp>
 #include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
 #include <Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp>
 #include <Hadrons/Modules/MUtilities/PrecisionCast.hpp>
 #include <Hadrons/Modules/MUtilities/RandomVectors.hpp>
-#include <Hadrons/Modules/MUtilities/TestSeqGamma.hpp>
-#include <Hadrons/Modules/MUtilities/TestSeqConserved.hpp>
-#include <Hadrons/Modules/MLoop/NoiseLoop.hpp>
 #include <Hadrons/Modules/MScalar/FreeProp.hpp>
-#include <Hadrons/Modules/MScalar/VPCounterTerms.hpp>
-#include <Hadrons/Modules/MScalar/ScalarVP.hpp>
 #include <Hadrons/Modules/MScalar/Scalar.hpp>
 #include <Hadrons/Modules/MScalar/ChargedProp.hpp>
 #include <Hadrons/Modules/MNPR/Bilinear.hpp>
@@ -58,7 +51,6 @@
 #include <Hadrons/Modules/MAction/ScaledDWF.hpp>
 #include <Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
 #include <Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
-#include <Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
 #include <Hadrons/Modules/MScalarSUN/Div.hpp>
 #include <Hadrons/Modules/MScalarSUN/TrMag.hpp>
 #include <Hadrons/Modules/MScalarSUN/EMT.hpp>
--- a/Hadrons/Modules/MAction/DWF.cc
+++ b/Hadrons/Modules/MAction/DWF.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/DWF.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MAction/DWF.hpp
+++ b/Hadrons/Modules/MAction/DWF.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/DWF.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
@@ -112,8 +112,6 @@ void TDWF<FImpl>::setup(void)
                 << par().mass << ", M5= " << par().M5 << " and Ls= "
                 << par().Ls << " using gauge field '" << par().gauge << "'"
                 << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                 
    auto &U    = envGet(GaugeField, par().gauge);
    auto &g4   = *envGetGrid(FermionField);
@@ -121,8 +119,26 @@ void TDWF<FImpl>::setup(void)
    auto &g5   = *envGetGrid(FermionField, par().Ls);
    auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
    typename DomainWallFermion<FImpl>::ImplParams implParams;
+    if (!par().boundary.empty())
+    {
        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
    envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
                     grb5, g4, grb4, par().mass, par().M5, implParams);
 }
--- a/Hadrons/Modules/MAction/MobiusDWF.cc
+++ b/Hadrons/Modules/MAction/MobiusDWF.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/MobiusDWF.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MAction/MobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/MobiusDWF.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/MobiusDWF.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -112,8 +112,6 @@ void TMobiusDWF<FImpl>::setup(void)
                 << ", b= " << par().b << ", c= " << par().c
                 << " using gauge field '" << par().gauge << "'"
                 << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                 
    auto &U    = envGet(GaugeField, par().gauge);
    auto &g4   = *envGetGrid(FermionField);
@@ -121,8 +119,26 @@ void TMobiusDWF<FImpl>::setup(void)
    auto &g5   = *envGetGrid(FermionField, par().Ls);
    auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
    typename MobiusFermion<FImpl>::ImplParams implParams;
+    if (!par().boundary.empty())
+    {
        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
    envCreateDerived(FMat, MobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                     grb5, g4, grb4, par().mass, par().M5, par().b, par().c,
                     implParams);
--- a/Hadrons/Modules/MAction/ScaledDWF.cc
+++ b/Hadrons/Modules/MAction/ScaledDWF.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/ScaledDWF.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MAction/ScaledDWF.hpp
+++ b/Hadrons/Modules/MAction/ScaledDWF.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/ScaledDWF.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -111,8 +111,6 @@ void TScaledDWF<FImpl>::setup(void)
                 << ", scale= " << par().scale
                 << " using gauge field '" << par().gauge << "'"
                 << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;

    auto &U    = envGet(GaugeField, par().gauge);
    auto &g4   = *envGetGrid(FermionField);
@@ -120,8 +118,26 @@ void TScaledDWF<FImpl>::setup(void)
    auto &g5   = *envGetGrid(FermionField, par().Ls);
    auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
    typename ScaledShamirFermion<FImpl>::ImplParams implParams;
+    if (!par().boundary.empty())
+    {
        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
    envCreateDerived(FMat, ScaledShamirFermion<FImpl>, getName(), par().Ls, U, g5,
                     grb5, g4, grb4, par().mass, par().M5, par().scale,
                     implParams);
--- a/Hadrons/Modules/MAction/Wilson.cc
+++ b/Hadrons/Modules/MAction/Wilson.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/Wilson.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MAction/Wilson.hpp
+++ b/Hadrons/Modules/MAction/Wilson.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/Wilson.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
@@ -109,15 +109,31 @@ void TWilson<FImpl>::setup(void)
 {
    LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass
                 << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                 
    auto &U      = envGet(GaugeField, par().gauge);
    auto &grid   = *envGetGrid(FermionField);
    auto &gridRb = *envGetRbGrid(FermionField);
    typename WilsonFermion<FImpl>::ImplParams implParams;
+    if (!par().boundary.empty())
+    {
        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
    envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
                     par().mass, implParams);
 }
--- a/Hadrons/Modules/MAction/WilsonClover.cc
+++ b/Hadrons/Modules/MAction/WilsonClover.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/WilsonClover.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MAction/WilsonClover.hpp
+++ b/Hadrons/Modules/MAction/WilsonClover.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/WilsonClover.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
@@ -112,17 +112,34 @@ void TWilsonClover<FImpl>::setup(void)
 {
    LOG(Message) << "Setting up Wilson clover fermion matrix with m= " << par().mass
                 << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary 
-                 << std::endl;
    LOG(Message) << "Clover term csw_r: " << par().csw_r
                 << " csw_t: " << par().csw_t
                 << std::endl;
+                 
    auto &U      = envGet(GaugeField, par().gauge);
    auto &grid   = *envGetGrid(FermionField);
    auto &gridRb = *envGetRbGrid(FermionField);
    typename WilsonCloverFermion<FImpl>::ImplParams implParams;
+    if (!par().boundary.empty())
+    {
        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
    envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid,
                     gridRb, par().mass, par().csw_r, par().csw_t, 
                     par().clover_anisotropy, implParams); 
--- a/Hadrons/Modules/MAction/ZMobiusDWF.cc
+++ b/Hadrons/Modules/MAction/ZMobiusDWF.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/ZMobiusDWF.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MAction/ZMobiusDWF.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -118,10 +118,7 @@ void TZMobiusDWF<FImpl>::setup(void)
    {
        LOG(Message) << "  omega[" << i << "]= " << par().omega[i] << std::endl;
    }
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;

-    env().createGrid(par().Ls);
    auto &U    = envGet(GaugeField, par().gauge);
    auto &g4   = *envGetGrid(FermionField);
    auto &grb4 = *envGetRbGrid(FermionField);
@@ -129,8 +126,26 @@ void TZMobiusDWF<FImpl>::setup(void)
    auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
    auto omega = par().omega;
    typename ZMobiusFermion<FImpl>::ImplParams implParams;
+    if (!par().boundary.empty())
+    {
        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
    envCreateDerived(FMat, ZMobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                     grb5, g4, grb4, par().mass, par().M5, omega,
                     par().b, par().c, implParams);
--- a/Hadrons/Modules/MContraction/A2AAslashField.cc
+++ b/Hadrons/Modules/MContraction/A2AAslashField.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/A2AAslashField.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MContraction/A2AAslashField.hpp
+++ b/Hadrons/Modules/MContraction/A2AAslashField.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/A2AAslashField.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MContraction/A2ALoop.cc
+++ b/Hadrons/Modules/MContraction/A2ALoop.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MLoop/NoiseLoop.cc
+Source file: Hadrons/Modules/MContraction/A2ALoop.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -25,11 +25,10 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Hadrons/Modules/MLoop/NoiseLoop.hpp>
+#include <Hadrons/Modules/MContraction/A2ALoop.hpp>

 using namespace Grid;
 using namespace Hadrons;
-using namespace MLoop;
-
-template class Grid::Hadrons::MLoop::TNoiseLoop<FIMPL>;
+using namespace MContraction;

+template class Grid::Hadrons::MContraction::TA2ALoop<FIMPL>;
--- a/Hadrons/Modules/MContraction/A2ALoop.hpp
+++ b/Hadrons/Modules/MContraction/A2ALoop.hpp
@@ -2,12 +2,11 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MLoop/NoiseLoop.hpp
+Source file: Hadrons/Modules/MContraction/A2ALoop.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
-Author: Lanny91 <andrew.lawson@gmail.com>

 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -26,9 +25,8 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-
-#ifndef Hadrons_MLoop_NoiseLoop_hpp_
-#define Hadrons_MLoop_NoiseLoop_hpp_
+#ifndef Hadrons_MContraction_A2ALoop_hpp_
+#define Hadrons_MContraction_A2ALoop_hpp_

 #include <Hadrons/Global.hpp>
 #include <Hadrons/Module.hpp>
@@ -36,99 +34,90 @@ See the full license in the file "LICENSE" in the top level distribution directo

 BEGIN_HADRONS_NAMESPACE

-/*
- 
- Noise loop propagator
- -----------------------------
- * loop_x = q_x * adj(eta_x)
- 
- * options:
- - q = Result of inversion on noise source.
- - eta = noise source.
-
- */
-
-
 /******************************************************************************
- *                         NoiseLoop                                          *
+ *                  From closed loop from all-to-all vectors                  *
 ******************************************************************************/
-BEGIN_MODULE_NAMESPACE(MLoop)
+BEGIN_MODULE_NAMESPACE(MContraction)

-class NoiseLoopPar: Serializable
+class A2ALoopPar: Serializable
 {
 public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(NoiseLoopPar,
-                                    std::string, q,
-                                    std::string, eta);
+    GRID_SERIALIZABLE_CLASS_MEMBERS(A2ALoopPar,
+                                    std::string, left,
+                                    std::string, right);
 };

 template <typename FImpl>
-class TNoiseLoop: public Module<NoiseLoopPar>
+class TA2ALoop: public Module<A2ALoopPar>
 {
 public:
    FERM_TYPE_ALIASES(FImpl,);
 public:
    // constructor
-    TNoiseLoop(const std::string name);
+    TA2ALoop(const std::string name);
    // destructor
-    virtual ~TNoiseLoop(void) {};
+    virtual ~TA2ALoop(void) {};
    // dependency relation
    virtual std::vector<std::string> getInput(void);
    virtual std::vector<std::string> getOutput(void);
-protected:
    // setup
    virtual void setup(void);
    // execution
    virtual void execute(void);
 };

-MODULE_REGISTER_TMP(NoiseLoop, TNoiseLoop<FIMPL>, MLoop);
+MODULE_REGISTER_TMP(A2ALoop, TA2ALoop<FIMPL>, MContraction);

 /******************************************************************************
- *                 TNoiseLoop implementation                                  *
+ *                        TA2ALoop implementation                             *
 ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
 template <typename FImpl>
-TNoiseLoop<FImpl>::TNoiseLoop(const std::string name)
-: Module<NoiseLoopPar>(name)
+TA2ALoop<FImpl>::TA2ALoop(const std::string name)
+: Module<A2ALoopPar>(name)
 {}

 // dependencies/products ///////////////////////////////////////////////////////
 template <typename FImpl>
-std::vector<std::string> TNoiseLoop<FImpl>::getInput(void)
+std::vector<std::string> TA2ALoop<FImpl>::getInput(void)
 {
-    std::vector<std::string> in = {par().q, par().eta};
+    std::vector<std::string> in = {par().left, par().right};
    
    return in;
 }

 template <typename FImpl>
-std::vector<std::string> TNoiseLoop<FImpl>::getOutput(void)
+std::vector<std::string> TA2ALoop<FImpl>::getOutput(void)
 {
-    std::vector<std::string> out = {getName()};
+    std::vector<std::string> out = {};
    
    return out;
 }

 // setup ///////////////////////////////////////////////////////////////////////
 template <typename FImpl>
-void TNoiseLoop<FImpl>::setup(void)
+void TA2ALoop<FImpl>::setup(void)
 {
    envCreateLat(PropagatorField, getName());
 }

 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
-void TNoiseLoop<FImpl>::execute(void)
+void TA2ALoop<FImpl>::execute(void)
 {
    auto &loop  = envGet(PropagatorField, getName());
-    auto &q    = envGet(PropagatorField, par().q);
-    auto &eta  = envGet(PropagatorField, par().eta);
-    loop = q*adj(eta);
+    auto &left  = envGet(std::vector<FermionField>, par().left);
+    auto &right = envGet(std::vector<FermionField>, par().right);
+
+    loop = zero;
+    for (unsigned int i = 0; i < left.size(); ++i)
+    {
+        loop += outerProduct(adj(left[i]), right[i]);
+    }
 }

 END_MODULE_NAMESPACE

 END_HADRONS_NAMESPACE

-#endif // Hadrons_MLoop_NoiseLoop_hpp_
+#endif // Hadrons_MContraction_A2ALoop_hpp_
--- a/Hadrons/Modules/MContraction/A2AMesonField.cc
+++ b/Hadrons/Modules/MContraction/A2AMesonField.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/A2AMesonField.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: paboyle <paboyle@ph.ed.ac.uk>
--- a/Hadrons/Modules/MContraction/A2AMesonField.hpp
+++ b/Hadrons/Modules/MContraction/A2AMesonField.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/A2AMesonField.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
--- a/Hadrons/Modules/MContraction/Baryon.cc
+++ b/Hadrons/Modules/MContraction/Baryon.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/Baryon.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MContraction/Baryon.hpp
+++ b/Hadrons/Modules/MContraction/Baryon.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/Baryon.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/DiscLoop.cc
+++ b/Hadrons/Modules/MContraction/DiscLoop.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/DiscLoop.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MContraction/DiscLoop.hpp
+++ b/Hadrons/Modules/MContraction/DiscLoop.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/DiscLoop.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/Gamma3pt.cc
+++ b/Hadrons/Modules/MContraction/Gamma3pt.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/Gamma3pt.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MContraction/Gamma3pt.hpp
+++ b/Hadrons/Modules/MContraction/Gamma3pt.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/Gamma3pt.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/Meson.cc
+++ b/Hadrons/Modules/MContraction/Meson.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/Meson.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MContraction/Meson.hpp
+++ b/Hadrons/Modules/MContraction/Meson.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MContraction/Meson.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Lanny91 <andrew.lawson@gmail.com>
--- a/Hadrons/Modules/MContraction/WeakEye3pt.cc
+++ b/Hadrons/Modules/MContraction/WeakEye3pt.cc
@@ -2,9 +2,9 @@

 Grid physics library, www.github.com/paboyle/Grid 

-Source file: Hadrons/Modules/MScalarSUN/ShiftProbe.cc
+Source file: Hadrons/Modules/MContraction/WeakEye3pt.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

@@ -25,15 +25,10 @@ with this program; if not, write to the Free Software Foundation, Inc.,
 See the full license in the file "LICENSE" in the top level distribution directory
 *************************************************************************************/
 /*  END LEGAL */
-#include <Hadrons/Modules/MScalarSUN/ShiftProbe.hpp>
+#include <Hadrons/Modules/MContraction/WeakEye3pt.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>>;
+using namespace MContraction;

+template class Grid::Hadrons::MContraction::TWeakEye3pt<FIMPL>;
--- a/Hadrons/Modules/MContraction/WeakEye3pt.hpp
+++ b/Hadrons/Modules/MContraction/WeakEye3pt.hpp
@@ -0,0 +1,200 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MContraction/WeakEye3pt.hpp
+
+Copyright (C) 2015-2019
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Lanny91 <andrew.lawson@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef Hadrons_MContraction_WeakEye3pt_hpp_
+#define Hadrons_MContraction_WeakEye3pt_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ * Weak Hamiltonian meson 3-pt diagrams, eye topologies.
+ * 
+ * Schematics:       loop                 |                  
+ *                  /-<-¬                 |                             
+ *                 /     \                |            qbl     G     qbr
+ *                 \     /                |        /----<------*------<----¬         
+ *            qbl   \   /    qbr          |       /          /-*-¬          \
+ *       /-----<-----* *-----<----¬       |      /          /  G  \          \
+ *  gIn *            G G           * gOut | gIn *           \     /  loop    * gOut
+ *       \                        /       |      \           \->-/          /   
+ *        \                      /        |       \                        /       
+ *         \---------->---------/         |        \----------->----------/        
+ *                   qs                   |                   qs                  
+ *                                        |
+ *                one trace               |                two traces
+ * 
+ * one trace : tr(qbr*gOut*qs*adj(gIn)*g5*adj(qbl)*g5*G*loop*G*qbr*gOut)
+ * two traces: tr(qbr*gOut*qs*adj(gIn)*g5*adj(qbl)*g5*G)*tr(loop*G)
+ * 
+ */
+
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class WeakEye3ptPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WeakEye3ptPar,
+                                    std::string,    qBarLeft,
+                                    std::string,    qBarRight,
+                                    std::string,    qSpectator,
+                                    std::string,    loop,
+                                    unsigned int,   tOut,
+                                    Gamma::Algebra, gammaIn,
+                                    Gamma::Algebra, gammaOut,
+                                    std::string,    output);
+};
+
+template <typename FImpl>
+class TWeakEye3pt: public Module<WeakEye3ptPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    class Metadata: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata,
+                                        Gamma::Algebra, in,
+                                        Gamma::Algebra, out,
+                                        Gamma::Algebra, op,
+                                        unsigned int,   trace);
+    };
+    typedef Correlator<Metadata> Result;
+public:
+    // constructor
+    TWeakEye3pt(const std::string name);
+    // destructor
+    virtual ~TWeakEye3pt(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(WeakEye3pt, TWeakEye3pt<FIMPL>, MContraction);
+
+/******************************************************************************
+ *                        TWeakEye3pt implementation                          *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWeakEye3pt<FImpl>::TWeakEye3pt(const std::string name)
+: Module<WeakEye3ptPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWeakEye3pt<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().qBarLeft, par().qBarRight, 
+                                   par().qSpectator};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWeakEye3pt<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWeakEye3pt<FImpl>::setup(void)
+{
+    envTmpLat(ComplexField, "corr");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWeakEye3pt<FImpl>::execute(void)
+{
+    LOG(Message) << "Computing mesonic weak 3pt contractions, eye topologies" << std::endl;
+    LOG(Message) << "gIn : " << par().gammaIn << std::endl;
+    LOG(Message) << "gOut: " << par().gammaIn << std::endl;
+    LOG(Message) << "tOut: " << par().tOut << std::endl;
+    LOG(Message) << "qbl : " << par().qBarLeft << std::endl;
+    LOG(Message) << "qbr : " << par().qBarRight << std::endl;
+    LOG(Message) << "qs  : " << par().qSpectator << std::endl;
+    LOG(Message) << "loop: " << par().loop << std::endl;
+
+    std::vector<Result> result;
+    Result              r;
+    auto                &qbl  = envGet(PropagatorField, par().qBarLeft);
+    auto                &qbr  = envGet(PropagatorField, par().qBarRight);
+    auto                &loop = envGet(PropagatorField, par().loop);
+    auto                &qs   = envGet(SlicedPropagator, par().qSpectator);
+    auto                qst   = qs[par().tOut];
+    Gamma               gIn(par().gammaIn), gOut(par().gammaOut);
+    Gamma               g5(Gamma::Algebra::Gamma5);
+
+    envGetTmp(ComplexField, corr);
+    r.info.in  = par().gammaIn;
+    r.info.out = par().gammaOut;
+    for (auto &G: Gamma::gall)
+    {
+        SlicedComplex buf;
+
+        r.info.op = G.g;
+        // one trace
+        corr = trace(qbr*gOut*qst*adj(gIn)*g5*adj(qbl)*g5*G*loop*G*qbr*gOut);
+        sliceSum(corr, buf, Tp);
+        r.corr.clear();
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            r.corr.push_back(TensorRemove(buf[t]));
+        }
+        r.info.trace = 1;
+        result.push_back(r);
+        // two traces
+        corr = trace(qbr*gOut*qst*adj(gIn)*g5*adj(qbl)*g5*G)*trace(loop*G);
+        sliceSum(corr, buf, Tp);
+        r.corr.clear();
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            r.corr.push_back(TensorRemove(buf[t]));
+        }
+        r.info.trace = 2;
+        result.push_back(r);
+    }
+    saveResult(par().output, "weakEye3pt", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MContraction_WeakEye3pt_hpp_
--- a/Hadrons/Modules/MContraction/WeakNonEye3pt.cc
+++ b/Hadrons/Modules/MContraction/WeakNonEye3pt.cc
@@ -0,0 +1,34 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MContraction/WeakNonEye3pt.cc
+
+Copyright (C) 2015-2019
+
+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 <Hadrons/Modules/MContraction/WeakNonEye3pt.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TWeakNonEye3pt<FIMPL>;
--- a/Hadrons/Modules/MContraction/WeakNonEye3pt.hpp
+++ b/Hadrons/Modules/MContraction/WeakNonEye3pt.hpp
@@ -0,0 +1,198 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: Hadrons/Modules/MContraction/WeakNonEye3pt.hpp
+
+Copyright (C) 2015-2019
+
+Author: Antonin Portelli <antonin.portelli@me.com>
+Author: Lanny91 <andrew.lawson@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+#ifndef Hadrons_MContraction_WeakNonEye3pt_hpp_
+#define Hadrons_MContraction_WeakNonEye3pt_hpp_
+
+#include <Hadrons/Global.hpp>
+#include <Hadrons/Module.hpp>
+#include <Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ * Weak Hamiltonian meson 3-pt diagrams, non-eye topologies.
+ * 
+ * Schematic:     
+ *            qbl           qbr            |            qbl             qbr
+ *          /--<--¬       /--<--¬          |          /--<--¬         /--<--¬       
+ *         /       \     /       \         |         /       \       /       \      
+ *        /         \   /         \        |        /         \     /         \     
+ *       /           \ /           \       |       /           \   /           \    
+ *  gIn *             * G           * gOut |  gIn *           G * * G           * gOut
+ *      \             * G           |      |       \           /   \           /
+ *       \           / \           /       |        \         /     \         /    
+ *        \         /   \         /        |         \       /       \       /  
+ *         \       /     \       /         |          \-->--/         \-->--/      
+ *          \-->--/       \-->--/          |            ql              qr 
+ *            ql            qr             |
+ *               one trace                 |              two traces
+ *
+ * one trace : tr(ql*adj(gIn)*g5*adj(qbl)*g5*G*qbr*gOut*g5*adj(qr)*g5*G)
+ * two traces: tr(ql*adj(gIn)*g5*adj(qbl)*g5*G)*tr(qbr*gOut*g5*adj(qr)*g5*G)
+ * 
+ */
+
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class WeakNonEye3ptPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WeakNonEye3ptPar,
+                                    std::string,    qLeft,
+                                    std::string,    qBarLeft,
+                                    std::string,    qRight,
+                                    std::string,    qBarRight,
+                                    Gamma::Algebra, gammaIn,
+                                    Gamma::Algebra, gammaOut,
+                                    std::string,    output);
+};
+
+template <typename FImpl>
+class TWeakNonEye3pt: public Module<WeakNonEye3ptPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    class Metadata: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata,
+                                        Gamma::Algebra, in,
+                                        Gamma::Algebra, out,
+                                        Gamma::Algebra, op,
+                                        unsigned int,   trace);
+    };
+    typedef Correlator<Metadata> Result;
+public:
+    // constructor
+    TWeakNonEye3pt(const std::string name);
+    // destructor
+    virtual ~TWeakNonEye3pt(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(WeakNonEye3pt, TWeakNonEye3pt<FIMPL>, MContraction);
+
+/******************************************************************************
+ *                 TWeakNonEye3pt implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWeakNonEye3pt<FImpl>::TWeakNonEye3pt(const std::string name)
+: Module<WeakNonEye3ptPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWeakNonEye3pt<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().qLeft, par().qBarLeft, 
+                                   par().qRight, par().qBarRight};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWeakNonEye3pt<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWeakNonEye3pt<FImpl>::setup(void)
+{
+    envTmpLat(ComplexField, "corr");
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWeakNonEye3pt<FImpl>::execute(void)
+{
+    LOG(Message) << "Computing mesonic weak 3pt contractions, non-eye topologies" << std::endl;
+    LOG(Message) << "gIn : " << par().gammaIn << std::endl;
+    LOG(Message) << "gOut: " << par().gammaIn << std::endl;
+    LOG(Message) << "ql  : " << par().qLeft << std::endl;
+    LOG(Message) << "qbl : " << par().qBarLeft << std::endl;
+    LOG(Message) << "qr  : " << par().qRight << std::endl;
+    LOG(Message) << "qbr : " << par().qBarRight << std::endl;
+
+    std::vector<Result> result;
+    Result              r;
+    auto                &ql  = envGet(PropagatorField, par().qLeft);
+    auto                &qbl = envGet(PropagatorField, par().qBarLeft);
+    auto                &qr  = envGet(PropagatorField, par().qRight);
+    auto                &qbr = envGet(PropagatorField, par().qBarRight);
+    Gamma               gIn(par().gammaIn), gOut(par().gammaOut);
+    Gamma               g5(Gamma::Algebra::Gamma5);
+
+    envGetTmp(ComplexField, corr);
+    r.info.in  = par().gammaIn;
+    r.info.out = par().gammaOut;
+    for (auto &G: Gamma::gall)
+    {
+        SlicedComplex buf;
+
+        r.info.op = G.g;
+        // one trace
+        corr = trace(ql*adj(gIn)*g5*adj(qbl)*g5*G*qbr*gOut*g5*adj(qr)*g5*G);
+        sliceSum(corr, buf, Tp);
+        r.corr.clear();
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            r.corr.push_back(TensorRemove(buf[t]));
+        }
+        r.info.trace = 1;
+        result.push_back(r);
+        // two traces
+        corr = trace(ql*adj(gIn)*g5*adj(qbl)*g5*G)*trace(qbr*gOut*g5*adj(qr)*g5*G);
+        sliceSum(corr, buf, Tp);
+        r.corr.clear();
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            r.corr.push_back(TensorRemove(buf[t]));
+        }
+        r.info.trace = 2;
+        result.push_back(r);
+    }
+    saveResult(par().output, "weakNonEye3pt", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MContraction_WeakNonEye3pt_hpp_
--- a/Hadrons/Modules/MFermion/FreeProp.cc
+++ b/Hadrons/Modules/MFermion/FreeProp.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MFermion/FreeProp.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <vmg1n14@soton.ac.uk>
--- a/Hadrons/Modules/MFermion/FreeProp.hpp
+++ b/Hadrons/Modules/MFermion/FreeProp.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MFermion/FreeProp.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <vmg1n14@soton.ac.uk>
--- a/Hadrons/Modules/MFermion/GaugeProp.cc
+++ b/Hadrons/Modules/MFermion/GaugeProp.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MFermion/GaugeProp.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MFermion/GaugeProp.hpp
+++ b/Hadrons/Modules/MFermion/GaugeProp.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MFermion/GaugeProp.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
--- a/Hadrons/Modules/MGauge/Electrify.cc
+++ b/Hadrons/Modules/MGauge/Electrify.cc
@@ -4,9 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/Electrify.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <Vera.Guelpers@ed.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
--- a/Hadrons/Modules/MGauge/Electrify.hpp
+++ b/Hadrons/Modules/MGauge/Electrify.hpp
@@ -4,9 +4,11 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/Electrify.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

+Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Vera Guelpers <Vera.Guelpers@ed.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
--- a/Hadrons/Modules/MGauge/FundtoHirep.cc
+++ b/Hadrons/Modules/MGauge/FundtoHirep.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/FundtoHirep.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Guido Cossu <guido.cossu@ed.ac.uk>
--- a/Hadrons/Modules/MGauge/FundtoHirep.hpp
+++ b/Hadrons/Modules/MGauge/FundtoHirep.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/FundtoHirep.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: pretidav <david.preti@csic.es>
--- a/Hadrons/Modules/MGauge/GaugeFix.cc
+++ b/Hadrons/Modules/MGauge/GaugeFix.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/GaugeFix.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
--- a/Hadrons/Modules/MGauge/GaugeFix.hpp
+++ b/Hadrons/Modules/MGauge/GaugeFix.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/GaugeFix.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: Peter Boyle <paboyle@ph.ed.ac.uk>
--- a/Hadrons/Modules/MGauge/Random.cc
+++ b/Hadrons/Modules/MGauge/Random.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/Random.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MGauge/Random.hpp
+++ b/Hadrons/Modules/MGauge/Random.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/Random.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MGauge/StochEm.cc
+++ b/Hadrons/Modules/MGauge/StochEm.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/StochEm.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
--- a/Hadrons/Modules/MGauge/StochEm.hpp
+++ b/Hadrons/Modules/MGauge/StochEm.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/StochEm.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>
 Author: James Harrison <jch1g10@soton.ac.uk>
--- a/Hadrons/Modules/MGauge/StoutSmearing.cc
+++ b/Hadrons/Modules/MGauge/StoutSmearing.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/StoutSmearing.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MGauge/StoutSmearing.hpp
+++ b/Hadrons/Modules/MGauge/StoutSmearing.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/StoutSmearing.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MGauge/Unit.cc
+++ b/Hadrons/Modules/MGauge/Unit.cc
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/Unit.cc

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Hadrons/Modules/MGauge/Unit.hpp
+++ b/Hadrons/Modules/MGauge/Unit.hpp
@@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid

 Source file: Hadrons/Modules/MGauge/Unit.hpp

-Copyright (C) 2015-2018
+Copyright (C) 2015-2019

 Author: Antonin Portelli <antonin.portelli@me.com>

--- a/Show More
+++ b/Show More