Merge branch 'develop' into feature/hmc_generalise

Makefile.am

@@ -3,6 +3,8 @@ SUBDIRS = lib benchmarks tests extras
 
 include $(top_srcdir)/doxygen.inc
 
+bin_SCRIPTS=grid-config
+
 tests: all
 	$(MAKE) -C tests tests
 

benchmarks/Benchmark_dwf.cc

@@ -145,6 +145,7 @@ int main (int argc, char ** argv)
   RealD M5  =1.8;
 
   RealD NP = UGrid->_Nprocessors;
+  RealD NN = UGrid->NodeCount();
 
   std::cout << GridLogMessage<< "*****************************************************************" <<std::endl;
   std::cout << GridLogMessage<< "* Kernel options --dslash-generic, --dslash-unroll, --dslash-asm" <<std::endl;
@@ -154,6 +155,10 @@ int main (int argc, char ** argv)
   std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplex::Nsimd()<<std::endl;
   if ( sizeof(Real)==4 )   std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
   if ( sizeof(Real)==8 )   std::cout << GridLogMessage<< "* DOUBLE precision "<<std::endl;
+#ifdef GRID_OMP
+  if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsAndCompute ) std::cout << GridLogMessage<< "* Using Overlapped Comms/Compute" <<std::endl;
+  if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsThenCompute) std::cout << GridLogMessage<< "* Using sequential comms compute" <<std::endl;
+#endif
   if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptGeneric   ) std::cout << GridLogMessage<< "* Using GENERIC Nc WilsonKernels" <<std::endl;
   if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptHandUnroll) std::cout << GridLogMessage<< "* Using Nc=3       WilsonKernels" <<std::endl;
   if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptInlineAsm ) std::cout << GridLogMessage<< "* Using Asm Nc=3   WilsonKernels" <<std::endl;
@@ -183,6 +188,7 @@ int main (int argc, char ** argv)
     //    std::cout<<GridLogMessage << "norm ref    "<< norm2(ref)<<std::endl;
     std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
     std::cout<<GridLogMessage << "mflop/s per rank =  "<< flops/(t1-t0)/NP<<std::endl;
+    std::cout<<GridLogMessage << "mflop/s per node =  "<< flops/(t1-t0)/NN<<std::endl;
     err = ref-result; 
     std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<<std::endl;
 
@@ -219,6 +225,7 @@ int main (int argc, char ** argv)
     std::cout<<GridLogMessage << "Called half prec comms Dw "<<ncall<<" times in "<<t1-t0<<" us"<<std::endl;
     std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
     std::cout<<GridLogMessage << "mflop/s per rank =  "<< flops/(t1-t0)/NP<<std::endl;
+    std::cout<<GridLogMessage << "mflop/s per node =  "<< flops/(t1-t0)/NN<<std::endl;
     err = ref-result; 
     std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<<std::endl;
 
@@ -234,6 +241,10 @@ int main (int argc, char ** argv)
     std::cout << GridLogMessage<< "* Vectorising fifth dimension by "<<vComplex::Nsimd()<<std::endl;
     if ( sizeof(Real)==4 )   std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
     if ( sizeof(Real)==8 )   std::cout << GridLogMessage<< "* DOUBLE precision "<<std::endl;
+#ifdef GRID_OMP
+  if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsAndCompute ) std::cout << GridLogMessage<< "* Using Overlapped Comms/Compute" <<std::endl;
+  if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsThenCompute) std::cout << GridLogMessage<< "* Using sequential comms compute" <<std::endl;
+#endif
     if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptGeneric   ) std::cout << GridLogMessage<< "* Using GENERIC Nc WilsonKernels" <<std::endl;
     if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptHandUnroll) std::cout << GridLogMessage<< "* Using Nc=3       WilsonKernels" <<std::endl;
     if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptInlineAsm ) std::cout << GridLogMessage<< "* Using Asm Nc=3   WilsonKernels" <<std::endl;
@@ -265,6 +276,7 @@ int main (int argc, char ** argv)
     std::cout<<GridLogMessage << "Called Dw s_inner "<<ncall<<" times in "<<t1-t0<<" us"<<std::endl;
     std::cout<<GridLogMessage << "mflop/s =   "<< flops/(t1-t0)<<std::endl;
     std::cout<<GridLogMessage << "mflop/s per rank =  "<< flops/(t1-t0)/NP<<std::endl;
+    std::cout<<GridLogMessage << "mflop/s per node =  "<< flops/(t1-t0)/NN<<std::endl;
     //    std::cout<<GridLogMessage<< "res norms "<< norm2(result)<<" " <<norm2(sresult)<<std::endl;
     sDw.Report();
     RealD sum=0;
@@ -297,6 +309,10 @@ int main (int argc, char ** argv)
       std::cout << GridLogMessage<< "* Vectorising fifth dimension by "<<vComplex::Nsimd()<<std::endl;
       if ( sizeof(Real)==4 )   std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
       if ( sizeof(Real)==8 )   std::cout << GridLogMessage<< "* DOUBLE precision "<<std::endl;
+#ifdef GRID_OMP
+  if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsAndCompute ) std::cout << GridLogMessage<< "* Using Overlapped Comms/Compute" <<std::endl;
+  if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsThenCompute) std::cout << GridLogMessage<< "* Using sequential comms compute" <<std::endl;
+#endif
       if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptGeneric   ) 
 	std::cout << GridLogMessage<< "* Using GENERIC Nc WilsonKernels" <<std::endl;
       if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptHandUnroll) 
@@ -336,6 +352,7 @@ int main (int argc, char ** argv)
 
       std::cout<<GridLogMessage << "sDeo mflop/s =   "<< flops/(t1-t0)<<std::endl;
       std::cout<<GridLogMessage << "sDeo mflop/s per rank   "<< flops/(t1-t0)/NP<<std::endl;
+      std::cout<<GridLogMessage << "sDeo mflop/s per node   "<< flops/(t1-t0)/NN<<std::endl;
       sDw.Report();
 
       sDw.DhopEO(ssrc_o,sr_e,DaggerNo);
@@ -414,14 +431,15 @@ int main (int argc, char ** argv)
 
 
   // S-direction is INNERMOST and takes no part in the parity.
-  static int Opt;  // these are a temporary hack
-  static int Comms;  // these are a temporary hack
-
   std::cout << GridLogMessage<< "*********************************************************" <<std::endl;
   std::cout << GridLogMessage<< "* Benchmarking DomainWallFermionR::DhopEO                "<<std::endl;
   std::cout << GridLogMessage<< "* Vectorising space-time by "<<vComplex::Nsimd()<<std::endl;
   if ( sizeof(Real)==4 )   std::cout << GridLogMessage<< "* SINGLE precision "<<std::endl;
   if ( sizeof(Real)==8 )   std::cout << GridLogMessage<< "* DOUBLE precision "<<std::endl;
+#ifdef GRID_OMP
+  if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsAndCompute ) std::cout << GridLogMessage<< "* Using Overlapped Comms/Compute" <<std::endl;
+  if ( WilsonKernelsStatic::Comms == WilsonKernelsStatic::CommsThenCompute) std::cout << GridLogMessage<< "* Using sequential comms compute" <<std::endl;
+#endif
   if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptGeneric   ) std::cout << GridLogMessage<< "* Using GENERIC Nc WilsonKernels" <<std::endl;
   if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptHandUnroll) std::cout << GridLogMessage<< "* Using Nc=3       WilsonKernels" <<std::endl;
   if ( WilsonKernelsStatic::Opt == WilsonKernelsStatic::OptInlineAsm ) std::cout << GridLogMessage<< "* Using Asm Nc=3   WilsonKernels" <<std::endl;
@@ -442,6 +460,7 @@ int main (int argc, char ** argv)
 
     std::cout<<GridLogMessage << "Deo mflop/s =   "<< flops/(t1-t0)<<std::endl;
     std::cout<<GridLogMessage << "Deo mflop/s per rank   "<< flops/(t1-t0)/NP<<std::endl;
+    std::cout<<GridLogMessage << "Deo mflop/s per node   "<< flops/(t1-t0)/NN<<std::endl;
     Dw.Report();
   }
   Dw.DhopEO(src_o,r_e,DaggerNo);

configure.ac

@@ -1,5 +1,5 @@
 AC_PREREQ([2.63])
-AC_INIT([Grid], [0.6.0], [https://github.com/paboyle/Grid], [Grid])
+AC_INIT([Grid], [0.6.0-dev], [https://github.com/paboyle/Grid], [Grid])
 AC_CANONICAL_BUILD
 AC_CANONICAL_HOST
 AC_CANONICAL_TARGET
@@ -32,7 +32,7 @@ AC_TYPE_SIZE_T
 AC_TYPE_UINT32_T
 AC_TYPE_UINT64_T
 
-############### OpenMP 
+############### OpenMP
 AC_OPENMP
 ac_openmp=no
 if test "${OPENMP_CXXFLAGS}X" != "X"; then
@@ -63,8 +63,8 @@ AC_ARG_WITH([mpfr],
     [AM_CXXFLAGS="-I$with_mpfr/include $AM_CXXFLAGS"]
     [AM_LDFLAGS="-L$with_mpfr/lib $AM_LDFLAGS"])
 
-############### FFTW3 
+############### FFTW3
-AC_ARG_WITH([fftw],    
+AC_ARG_WITH([fftw],
             [AS_HELP_STRING([--with-fftw=prefix],
             [try this for a non-standard install prefix of the FFTW3 library])],
             [AM_CXXFLAGS="-I$with_fftw/include $AM_CXXFLAGS"]
@@ -77,9 +77,9 @@ AC_ARG_WITH([lime],
             [AM_CXXFLAGS="-I$with_lime/include $AM_CXXFLAGS"]
             [AM_LDFLAGS="-L$with_lime/lib $AM_LDFLAGS"])
 
-############### lapack 
+############### lapack
 AC_ARG_ENABLE([lapack],
-    [AC_HELP_STRING([--enable-lapack=yes|no|prefix], [enable LAPACK])], 
+    [AC_HELP_STRING([--enable-lapack=yes|no|prefix], [enable LAPACK])],
     [ac_LAPACK=${enable_lapack}], [ac_LAPACK=no])
 
 case ${ac_LAPACK} in
@@ -95,7 +95,7 @@ esac
 
 ############### FP16 conversions
 AC_ARG_ENABLE([sfw-fp16],
-    [AC_HELP_STRING([--enable-sfw-fp16=yes|no], [enable software fp16 comms])], 
+    [AC_HELP_STRING([--enable-sfw-fp16=yes|no], [enable software fp16 comms])],
     [ac_SFW_FP16=${enable_sfw_fp16}], [ac_SFW_FP16=yes])
 case ${ac_SFW_FP16} in
     yes)
@@ -130,7 +130,7 @@ AC_ARG_WITH([hdf5],
 
 ############### first-touch
 AC_ARG_ENABLE([numa],
-    [AC_HELP_STRING([--enable-numa=yes|no|prefix], [enable first touch numa opt])], 
+    [AC_HELP_STRING([--enable-numa=yes|no|prefix], [enable first touch numa opt])],
     [ac_NUMA=${enable_NUMA}],[ac_NUMA=no])
 
 case ${ac_NUMA} in
@@ -156,8 +156,8 @@ if test "${ac_MKL}x" != "nox"; then
 fi
 
 AC_SEARCH_LIBS([__gmpf_init], [gmp],
-               [AC_SEARCH_LIBS([mpfr_init], [mpfr], 
+               [AC_SEARCH_LIBS([mpfr_init], [mpfr],
-                               [AC_DEFINE([HAVE_LIBMPFR], [1], 
+                               [AC_DEFINE([HAVE_LIBMPFR], [1],
                                           [Define to 1 if you have the `MPFR' library])]
                                [have_mpfr=true], [AC_MSG_ERROR([MPFR library not found])])]
                [AC_DEFINE([HAVE_LIBGMP], [1], [Define to 1 if you have the `GMP' library])]
@@ -166,7 +166,7 @@ AC_SEARCH_LIBS([__gmpf_init], [gmp],
 if test "${ac_LAPACK}x" != "nox"; then
     AC_SEARCH_LIBS([LAPACKE_sbdsdc], [lapack], [],
                    [AC_MSG_ERROR("LAPACK enabled but library not found")])
-fi   
+fi
 
 AC_SEARCH_LIBS([fftw_execute], [fftw3],
                [AC_SEARCH_LIBS([fftwf_execute], [fftw3f], [],
@@ -334,7 +334,7 @@ case ${ac_COMMS} in
         comms_type='shmem'
      ;;
      *)
-        AC_MSG_ERROR([${ac_COMMS} unsupported --enable-comms option]); 
+        AC_MSG_ERROR([${ac_COMMS} unsupported --enable-comms option]);
      ;;
 esac
 case ${ac_COMMS} in
@@ -371,7 +371,7 @@ case ${ac_RNG} in
       AC_DEFINE([RNG_SITMO],[1],[RNG_SITMO] )
      ;;
      *)
-      AC_MSG_ERROR([${ac_RNG} unsupported --enable-rng option]); 
+      AC_MSG_ERROR([${ac_RNG} unsupported --enable-rng option]);
      ;;
 esac
 
@@ -388,7 +388,7 @@ case ${ac_TIMERS} in
       AC_DEFINE([TIMERS_OFF],[1],[TIMERS_OFF] )
      ;;
      *)
-      AC_MSG_ERROR([${ac_TIMERS} unsupported --enable-timers option]); 
+      AC_MSG_ERROR([${ac_TIMERS} unsupported --enable-timers option]);
      ;;
 esac
 
@@ -400,7 +400,7 @@ case ${ac_CHROMA} in
      yes|no)
      ;;
      *)
-       AC_MSG_ERROR([${ac_CHROMA} unsupported --enable-chroma option]); 
+       AC_MSG_ERROR([${ac_CHROMA} unsupported --enable-chroma option]);
      ;;
 esac
 
@@ -421,29 +421,23 @@ DX_INIT_DOXYGEN([$PACKAGE_NAME], [doxygen.cfg])
 
 ############### Ouput
 cwd=`pwd -P`; cd ${srcdir}; abs_srcdir=`pwd -P`; cd ${cwd}
+GRID_CXXFLAGS="$AM_CXXFLAGS $CXXFLAGS"
+GRID_LDFLAGS="$AM_LDFLAGS $LDFLAGS"
+GRID_LIBS=$LIBS
+GRID_SHORT_SHA=`git rev-parse --short HEAD`
+GRID_SHA=`git rev-parse HEAD`
+GRID_BRANCH=`git rev-parse --abbrev-ref HEAD`
 AM_CXXFLAGS="-I${abs_srcdir}/include $AM_CXXFLAGS"
 AM_CFLAGS="-I${abs_srcdir}/include $AM_CFLAGS"
 AM_LDFLAGS="-L${cwd}/lib $AM_LDFLAGS"
 AC_SUBST([AM_CFLAGS])
 AC_SUBST([AM_CXXFLAGS])
 AC_SUBST([AM_LDFLAGS])
-AC_CONFIG_FILES(Makefile)
+AC_SUBST([GRID_CXXFLAGS])
-AC_CONFIG_FILES(lib/Makefile)
+AC_SUBST([GRID_LDFLAGS])
-AC_CONFIG_FILES(tests/Makefile)
+AC_SUBST([GRID_LIBS])
-AC_CONFIG_FILES(tests/IO/Makefile)
+AC_SUBST([GRID_SHA])
-AC_CONFIG_FILES(tests/core/Makefile)
+AC_SUBST([GRID_BRANCH])
-AC_CONFIG_FILES(tests/debug/Makefile)
-AC_CONFIG_FILES(tests/forces/Makefile)
-AC_CONFIG_FILES(tests/hadrons/Makefile)
-AC_CONFIG_FILES(tests/hmc/Makefile)
-AC_CONFIG_FILES(tests/solver/Makefile)
-AC_CONFIG_FILES(tests/qdpxx/Makefile)
-AC_CONFIG_FILES(tests/smearing/Makefile)
-AC_CONFIG_FILES(tests/testu01/Makefile)
-AC_CONFIG_FILES(benchmarks/Makefile)
-AC_CONFIG_FILES(extras/Makefile)
-AC_CONFIG_FILES(extras/Hadrons/Makefile)
-AC_OUTPUT
 
 git_commit=`cd $srcdir && ./scripts/configure.commit`
 
@@ -452,7 +446,6 @@ Summary of configuration for $PACKAGE v$VERSION
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 ----- GIT VERSION -------------------------------------
 $git_commit
-
 ----- PLATFORM ----------------------------------------
 architecture (build)        : $build_cpu
 os (build)                  : $build_os
@@ -462,11 +455,11 @@ compiler vendor             : ${ax_cv_cxx_compiler_vendor}
 compiler version            : ${ax_cv_gxx_version}
 ----- BUILD OPTIONS -----------------------------------
 SIMD                        : ${ac_SIMD}${SIMD_GEN_WIDTH_MSG}
-Threading                   : ${ac_openmp} 
+Threading                   : ${ac_openmp}
 Communications type         : ${comms_type}
 Default precision           : ${ac_PRECISION}
 Software FP16 conversion    : ${ac_SFW_FP16}
-RNG choice                  : ${ac_RNG} 
+RNG choice                  : ${ac_RNG}
 GMP                         : `if test "x$have_gmp" = xtrue; then echo yes; else echo no; fi`
 LAPACK                      : ${ac_LAPACK}
 FFTW                        : `if test "x$have_fftw" = xtrue; then echo yes; else echo no; fi`
@@ -481,6 +474,29 @@ LDFLAGS:
 LIBS:
 `echo ${LIBS} | tr ' ' '\n' | sed 's/^-/    -/g'`
 -------------------------------------------------------" > grid.configure.summary
+
+GRID_SUMMARY="`cat grid.configure.summary`"
+AM_SUBST_NOTMAKE([GRID_SUMMARY])
+AC_SUBST([GRID_SUMMARY])
+
+AC_CONFIG_FILES([grid-config], [chmod +x grid-config])
+AC_CONFIG_FILES(Makefile)
+AC_CONFIG_FILES(lib/Makefile)
+AC_CONFIG_FILES(tests/Makefile)
+AC_CONFIG_FILES(tests/IO/Makefile)
+AC_CONFIG_FILES(tests/core/Makefile)
+AC_CONFIG_FILES(tests/debug/Makefile)
+AC_CONFIG_FILES(tests/forces/Makefile)
+AC_CONFIG_FILES(tests/hadrons/Makefile)
+AC_CONFIG_FILES(tests/hmc/Makefile)
+AC_CONFIG_FILES(tests/solver/Makefile)
+AC_CONFIG_FILES(tests/qdpxx/Makefile)
+AC_CONFIG_FILES(tests/testu01/Makefile)
+AC_CONFIG_FILES(benchmarks/Makefile)
+AC_CONFIG_FILES(extras/Makefile)
+AC_CONFIG_FILES(extras/Hadrons/Makefile)
+AC_OUTPUT
+
 echo ""
 cat grid.configure.summary
 echo ""

extras/Hadrons/Application.cc

@@ -162,7 +162,8 @@ void Application::saveParameterFile(const std::string parameterFileName)
 sizeString((size)*locVol_) << " (" << sizeString(size)  << "/site)"
 
 #define DEFINE_MEMPEAK \
-auto memPeak = [this](const std::vector<unsigned int> &program)\
+GeneticScheduler<unsigned int>::ObjFunc memPeak = \
+[this](const std::vector<unsigned int> &program)\
 {\
     unsigned int memPeak;\
     bool         msg;\

extras/Hadrons/Global.hpp

@@ -145,6 +145,15 @@ std::string typeName(void)
     return typeName(typeIdPt<T>());
 }
 
+// default writers/readers
+#ifdef HAVE_HDF5
+typedef Hdf5Reader CorrReader;
+typedef Hdf5Writer CorrWriter;
+#else
+typedef XmlReader CorrReader;
+typedef XmlWriter CorrWriter;
+#endif
+
 END_HADRONS_NAMESPACE
 
 #endif // Hadrons_Global_hpp_

extras/Hadrons/Modules.hpp

@@ -29,12 +29,20 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MAction/DWF.hpp>
 #include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
 #include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
+#include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
+#include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
 #include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
 #include <Grid/Hadrons/Modules/MGauge/Load.hpp>
 #include <Grid/Hadrons/Modules/MGauge/Random.hpp>
 #include <Grid/Hadrons/Modules/MGauge/Unit.hpp>
+#include <Grid/Hadrons/Modules/MLoop/NoiseLoop.hpp>
 #include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
 #include <Grid/Hadrons/Modules/MSource/Point.hpp>
 #include <Grid/Hadrons/Modules/MSource/SeqGamma.hpp>
+#include <Grid/Hadrons/Modules/MSource/Wall.hpp>
 #include <Grid/Hadrons/Modules/MSource/Z2.hpp>
 #include <Grid/Hadrons/Modules/Quark.hpp>

extras/Hadrons/Modules/MContraction/Baryon.hpp

@@ -112,7 +112,7 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
                  << " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
                  << par().q3 << "'" << std::endl;
     
-    XmlWriter             writer(par().output);
+    CorrWriter             writer(par().output);
     PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
     PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
     PropagatorField3      &q3 = *env().template getObject<PropagatorField3>(par().q2);

extras/Hadrons/Modules/MContraction/DiscLoop.hpp

@@ -0,0 +1,144 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/DiscLoop.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_DiscLoop_hpp_
+#define Hadrons_DiscLoop_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                DiscLoop                                    *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class DiscLoopPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(DiscLoopPar,
+                                    std::string,    q_loop,
+                                    Gamma::Algebra, gamma,
+                                    std::string,    output);
+};
+
+template <typename FImpl>
+class TDiscLoop: public Module<DiscLoopPar>
+{
+    TYPE_ALIASES(FImpl,);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        Gamma::Algebra, gamma,
+                                        std::vector<Complex>, corr);
+    };
+public:
+    // constructor
+    TDiscLoop(const std::string name);
+    // destructor
+    virtual ~TDiscLoop(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(DiscLoop, TDiscLoop<FIMPL>, MContraction);
+
+/******************************************************************************
+ *                       TDiscLoop implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TDiscLoop<FImpl>::TDiscLoop(const std::string name)
+: Module<DiscLoopPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TDiscLoop<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q_loop};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TDiscLoop<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TDiscLoop<FImpl>::setup(void)
+{
+    
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TDiscLoop<FImpl>::execute(void)
+{
+    LOG(Message) << "Computing disconnected loop contraction '" << getName() 
+                 << "' using '" << par().q_loop << "' with " << par().gamma 
+                 << " insertion." << std::endl;
+
+    CorrWriter            writer(par().output);
+    PropagatorField       &q_loop = *env().template getObject<PropagatorField>(par().q_loop);
+    LatticeComplex        c(env().getGrid());
+    Gamma                 gamma(par().gamma);
+    std::vector<TComplex> buf;
+    Result                result;
+
+    c = trace(gamma*q_loop);
+    sliceSum(c, buf, Tp);
+
+    result.gamma = par().gamma;
+    result.corr.resize(buf.size());
+    for (unsigned int t = 0; t < buf.size(); ++t)
+    {
+        result.corr[t] = TensorRemove(buf[t]);
+    }
+
+    write(writer, "disc", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_DiscLoop_hpp_

extras/Hadrons/Modules/MContraction/Gamma3pt.hpp

@@ -0,0 +1,170 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/Gamma3pt.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_Gamma3pt_hpp_
+#define Hadrons_Gamma3pt_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ * 3pt contraction with gamma matrix insertion.
+ *
+ * Schematic:
+ *
+ *             q2           q3
+ *        /----<------*------<----¬
+ *       /          gamma          \
+ *      /                           \
+ *   i *                            * f
+ *      \                          /
+ *       \                        /
+ *        \----------->----------/
+ *                   q1
+ *
+ *      trace(g5*q1*adj(q2)*g5*gamma*q3)
+ */
+
+/******************************************************************************
+ *                               Gamma3pt                                     *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+class Gamma3ptPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(Gamma3ptPar,
+                                    std::string,    q1,
+                                    std::string,    q2,
+                                    std::string,    q3,
+                                    Gamma::Algebra, gamma,
+                                    std::string,    output);
+};
+
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+class TGamma3pt: public Module<Gamma3ptPar>
+{
+    TYPE_ALIASES(FImpl1, 1);
+    TYPE_ALIASES(FImpl2, 2);
+    TYPE_ALIASES(FImpl3, 3);
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        Gamma::Algebra, gamma,
+                                        std::vector<Complex>, corr);
+    };
+public:
+    // constructor
+    TGamma3pt(const std::string name);
+    // destructor
+    virtual ~TGamma3pt(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Gamma3pt, ARG(TGamma3pt<FIMPL, FIMPL, FIMPL>), MContraction);
+
+/******************************************************************************
+ *                       TGamma3pt implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+TGamma3pt<FImpl1, FImpl2, FImpl3>::TGamma3pt(const std::string name)
+: Module<Gamma3ptPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getInput(void)
+{
+    std::vector<std::string> in = {par().q1, par().q2, par().q3};
+    
+    return in;
+}
+
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+std::vector<std::string> TGamma3pt<FImpl1, FImpl2, FImpl3>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+void TGamma3pt<FImpl1, FImpl2, FImpl3>::setup(void)
+{
+    
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl1, typename FImpl2, typename FImpl3>
+void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
+{
+    LOG(Message) << "Computing 3pt contractions '" << getName() << "' using"
+                 << " quarks '" << par().q1 << "', '" << par().q2 << "' and '"
+                 << par().q3 << "', with " << par().gamma << " insertion." 
+                 << std::endl;
+
+    CorrWriter            writer(par().output);
+    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
+    PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
+    PropagatorField3      &q3 = *env().template getObject<PropagatorField3>(par().q3);
+    LatticeComplex        c(env().getGrid());
+    Gamma                 g5(Gamma::Algebra::Gamma5);
+    Gamma                 gamma(par().gamma);
+    std::vector<TComplex> buf;
+    Result                result;
+
+    c = trace(g5*q1*adj(q2)*(g5*gamma)*q3);
+    sliceSum(c, buf, Tp);
+
+    result.gamma = par().gamma;
+    result.corr.resize(buf.size());
+    for (unsigned int t = 0; t < buf.size(); ++t)
+    {
+        result.corr[t] = TensorRemove(buf[t]);
+    }
+
+    write(writer, "gamma3pt", result);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_Gamma3pt_hpp_

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -6,8 +6,10 @@ Source file: extras/Hadrons/Modules/MContraction/Meson.hpp
 
 Copyright (C) 2015
 Copyright (C) 2016
+Copyright (C) 2017
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+        Andrew Lawson    <andrew.lawson1991@gmail.com>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -36,20 +38,39 @@ See the full license in the file "LICENSE" in the top level distribution directo
 
 BEGIN_HADRONS_NAMESPACE
 
+/*
+ 
+ Meson contractions
+ -----------------------------
+ 
+ * options:
+ - q1: input propagator 1 (string)
+ - q2: input propagator 2 (string)
+ - gammas: gamma products to insert at sink & source, pairs of gamma matrices 
+           (space-separated strings) in angled brackets (i.e. <g_sink g_src>),
+           in a sequence (e.g. "<Gamma5 Gamma5><Gamma5 GammaT>").
+
+           Special values: "all" - perform all possible contractions.
+ - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0."),
+        given as multiples of (2*pi) / L.
+*/
+
 /******************************************************************************
  *                                TMeson                                       *
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MContraction)
 
+typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
+
 class MesonPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar,
-                                    std::string,    q1,
+                                    std::string, q1,
-                                    std::string,    q2,
+                                    std::string, q2,
-                                    std::string,    output,
+                                    std::string, gammas,
-                                    Gamma::Algebra, gammaSource,
+                                    std::string, mom,
-                                    Gamma::Algebra, gammaSink);
+                                    std::string, output);
 };
 
 template <typename FImpl1, typename FImpl2>
@@ -61,7 +82,10 @@ public:
     class Result: Serializable
     {
     public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result, std::vector<Complex>, corr);
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                                        Gamma::Algebra, gamma_snk,
+                                        Gamma::Algebra, gamma_src,
+                                        std::vector<Complex>, corr);
     };
 public:
     // constructor
@@ -71,6 +95,7 @@ public:
     // dependencies/products
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+    virtual void parseGammaString(std::vector<GammaPair> &gammaList);
     // execution
     virtual void execute(void);
 };
@@ -103,6 +128,32 @@ std::vector<std::string> TMeson<FImpl1, FImpl2>::getOutput(void)
     return output;
 }
 
+template <typename FImpl1, typename FImpl2>
+void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
+{
+    // Determine gamma matrices to insert at source/sink.
+    if (par().gammas.compare("all") == 0)
+    {
+        // Do all contractions.
+        unsigned int n_gam = Ns * Ns;
+        gammaList.resize(n_gam*n_gam);
+        for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
+        {
+            for (unsigned int j = 1; j < Gamma::nGamma; j += 2)
+            {
+                gammaList.push_back(std::make_pair((Gamma::Algebra)i, 
+                                                   (Gamma::Algebra)j));
+            }
+        }
+    }
+    else
+    {
+        // Parse individual contractions from input string.
+        gammaList = strToVec<GammaPair>(par().gammas);
+    }
+}
+
+
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl1, typename FImpl2>
 void TMeson<FImpl1, FImpl2>::execute(void)
@@ -111,21 +162,44 @@ void TMeson<FImpl1, FImpl2>::execute(void)
                  << " quarks '" << par().q1 << "' and '" << par().q2 << "'"
                  << std::endl;
     
-    XmlWriter             writer(par().output);
+    CorrWriter              writer(par().output);
-    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
+    PropagatorField1       &q1 = *env().template getObject<PropagatorField1>(par().q1);
-    PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
+    PropagatorField2       &q2 = *env().template getObject<PropagatorField2>(par().q2);
-    LatticeComplex        c(env().getGrid());
+    LatticeComplex         c(env().getGrid());
-    Gamma                 gSrc(par().gammaSource), gSnk(par().gammaSink);
+    Gamma                  g5(Gamma::Algebra::Gamma5);
-    Gamma                 g5(Gamma::Algebra::Gamma5);
+    std::vector<GammaPair> gammaList;
-    std::vector<TComplex> buf;
+    std::vector<TComplex>  buf;
-    Result                result;
+    std::vector<Result>    result;
-    
+    std::vector<Real>      p;
-    c = trace(gSnk*q1*adj(gSrc)*g5*adj(q2)*g5);
+
-    sliceSum(c, buf, Tp);
+    p  = strToVec<Real>(par().mom);
-    result.corr.resize(buf.size());
+    LatticeComplex         ph(env().getGrid()), coor(env().getGrid());
-    for (unsigned int t = 0; t < buf.size(); ++t)
+    Complex                i(0.0,1.0);
+    ph = zero;
+    for(unsigned int mu = 0; mu < env().getNd(); mu++)
     {
-        result.corr[t] = TensorRemove(buf[t]);
+        LatticeCoordinate(coor, mu);
+        ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
+    }
+    ph = exp((Real)(2*M_PI)*i*ph);
+    
+    parseGammaString(gammaList);
+
+    result.resize(gammaList.size());
+    for (unsigned int i = 0; i < result.size(); ++i)
+    {
+        Gamma gSnk(gammaList[i].first);
+        Gamma gSrc(gammaList[i].second);
+        c = trace((g5*gSnk)*q1*(adj(gSrc)*g5)*adj(q2))*ph;
+        sliceSum(c, buf, Tp);
+
+        result[i].gamma_snk = gammaList[i].first;
+        result[i].gamma_src = gammaList[i].second;
+        result[i].corr.resize(buf.size());
+        for (unsigned int t = 0; t < buf.size(); ++t)
+        {
+            result[i].corr[t] = TensorRemove(buf[t]);
+        }
     }
     write(writer, "meson", result);
 }

extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp

@@ -0,0 +1,114 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_WeakHamiltonian_hpp_
+#define Hadrons_WeakHamiltonian_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         WeakHamiltonian                                    *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+/*******************************************************************************
+ * Utilities for contractions involving the Weak Hamiltonian.
+ ******************************************************************************/
+//// Sum and store correlator.
+#define MAKE_DIAG(exp, buf, res, n)\
+sliceSum(exp, buf, Tp);\
+res.name = (n);\
+res.corr.resize(buf.size());\
+for (unsigned int t = 0; t < buf.size(); ++t)\
+{\
+    res.corr[t] = TensorRemove(buf[t]);\
+}
+
+//// Contraction of mu index: use 'mu' variable in exp.
+#define SUM_MU(buf,exp)\
+buf = zero;\
+for (unsigned int mu = 0; mu < ndim; ++mu)\
+{\
+    buf += exp;\
+}
+
+enum 
+{
+  i_V = 0,
+  i_A = 1,
+  n_i = 2
+};
+
+class WeakHamiltonianPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WeakHamiltonianPar,
+                                    std::string, q1,
+                                    std::string, q2,
+                                    std::string, q3,
+                                    std::string, q4,
+                                    std::string, output);
+};
+
+#define MAKE_WEAK_MODULE(modname)\
+class T##modname: public Module<WeakHamiltonianPar>\
+{\
+public:\
+    TYPE_ALIASES(FIMPL,)\
+    class Result: Serializable\
+    {\
+    public:\
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,\
+                                        std::string, name,\
+                                        std::vector<Complex>, corr);\
+    };\
+public:\
+    /* constructor */ \
+    T##modname(const std::string name);\
+    /* destructor */ \
+    virtual ~T##modname(void) = default;\
+    /* dependency relation */ \
+    virtual std::vector<std::string> getInput(void);\
+    virtual std::vector<std::string> getOutput(void);\
+    /* setup */ \
+    virtual void setup(void);\
+    /* execution */ \
+    virtual void execute(void);\
+    std::vector<std::string> VA_label = {"V", "A"};\
+};\
+MODULE_REGISTER_NS(modname, T##modname, MContraction);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WeakHamiltonian_hpp_

extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc

@@ -0,0 +1,137 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+/*
+ * Weak Hamiltonian current-current contractions, Eye-type.
+ * 
+ * These contractions are generated by the Q1 and Q2 operators in the physical
+ * basis (see e.g. Fig 3 of arXiv:1507.03094).
+ * 
+ * Schematics:        q4                 |                  
+ *                  /-<-¬                |                             
+ *                 /     \               |             q2           q3
+ *                 \     /               |        /----<------*------<----¬                        
+ *            q2    \   /    q3          |       /          /-*-¬          \
+ *       /-----<-----* *-----<----¬      |      /          /     \          \
+ *    i *            H_W           * f   |   i *           \     /  q4      * f
+ *       \                        /      |      \           \->-/          /   
+ *        \                      /       |       \                        /       
+ *         \---------->---------/        |        \----------->----------/        
+ *                   q1                  |                   q1                  
+ *                                       |
+ *                Saucer (S)             |                  Eye (E)
+ * 
+ * S: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1]*q4*gL[mu][p_2])
+ * E: trace(q3*g5*q1*adj(q2)*g5*gL[mu][p_1])*trace(q4*gL[mu][p_2])
+ */
+
+/******************************************************************************
+ *                  TWeakHamiltonianEye implementation                        *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+TWeakHamiltonianEye::TWeakHamiltonianEye(const std::string name)
+: Module<WeakHamiltonianPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> TWeakHamiltonianEye::getInput(void)
+{
+    std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
+    
+    return in;
+}
+
+std::vector<std::string> TWeakHamiltonianEye::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void TWeakHamiltonianEye::setup(void)
+{
+
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void TWeakHamiltonianEye::execute(void)
+{
+    LOG(Message) << "Computing Weak Hamiltonian (Eye type) contractions '" 
+                 << getName() << "' using quarks '" << par().q1 << "', '" 
+                 << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
+                 << "'." << std::endl;
+
+    CorrWriter             writer(par().output);
+    PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
+    PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
+    PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
+    PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
+    Gamma g5            = Gamma(Gamma::Algebra::Gamma5);
+    LatticeComplex        expbuf(env().getGrid());
+    std::vector<TComplex> corrbuf;
+    std::vector<Result>   result(n_eye_diag);
+    unsigned int ndim   = env().getNd();
+
+    PropagatorField              tmp1(env().getGrid());
+    LatticeComplex               tmp2(env().getGrid());
+    std::vector<PropagatorField> S_body(ndim, tmp1);
+    std::vector<PropagatorField> S_loop(ndim, tmp1);
+    std::vector<LatticeComplex>  E_body(ndim, tmp2);
+    std::vector<LatticeComplex>  E_loop(ndim, tmp2);
+
+    // Setup for S-type contractions.
+    for (int mu = 0; mu < ndim; ++mu)
+    {
+        S_body[mu] = MAKE_SE_BODY(q1, q2, q3, GammaL(Gamma::gmu[mu]));
+        S_loop[mu] = MAKE_SE_LOOP(q4, GammaL(Gamma::gmu[mu]));
+    }
+
+    // Perform S-type contractions.    
+    SUM_MU(expbuf, trace(S_body[mu]*S_loop[mu]))
+    MAKE_DIAG(expbuf, corrbuf, result[S_diag], "HW_S")
+
+    // Recycle sub-expressions for E-type contractions.
+    for (unsigned int mu = 0; mu < ndim; ++mu)
+    {
+        E_body[mu] = trace(S_body[mu]);
+        E_loop[mu] = trace(S_loop[mu]);
+    }
+
+    // Perform E-type contractions.
+    SUM_MU(expbuf, E_body[mu]*E_loop[mu])
+    MAKE_DIAG(expbuf, corrbuf, result[E_diag], "HW_E")
+
+    write(writer, "HW_Eye", result);
+}

extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp

@@ -0,0 +1,58 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_WeakHamiltonianEye_hpp_
+#define Hadrons_WeakHamiltonianEye_hpp_
+
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         WeakHamiltonianEye                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+enum
+{
+    S_diag = 0,
+    E_diag = 1,
+    n_eye_diag = 2
+};
+
+// Saucer and Eye subdiagram contractions.
+#define MAKE_SE_BODY(Q_1, Q_2, Q_3, gamma) (Q_3*g5*Q_1*adj(Q_2)*g5*gamma)
+#define MAKE_SE_LOOP(Q_loop, gamma) (Q_loop*gamma)
+
+MAKE_WEAK_MODULE(WeakHamiltonianEye)
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WeakHamiltonianEye_hpp_

extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc

@@ -0,0 +1,139 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+/*
+ * Weak Hamiltonian current-current contractions, Non-Eye-type.
+ * 
+ * These contractions are generated by the Q1 and Q2 operators in the physical
+ * basis (see e.g. Fig 3 of arXiv:1507.03094).
+ * 
+ * Schematic:     
+ *            q2             q3          |           q2              q3
+ *          /--<--¬       /--<--¬        |        /--<--¬         /--<--¬       
+ *         /       \     /       \       |       /       \       /       \      
+ *        /         \   /         \      |      /         \     /         \     
+ *       /           \ /           \     |     /           \   /           \    
+ *    i *             * H_W         *  f |  i *             * * H_W         * f 
+ *      \             *             |    |     \           /   \           /
+ *       \           / \           /     |      \         /     \         /    
+ *        \         /   \         /      |       \       /       \       /  
+ *         \       /     \       /       |        \-->--/         \-->--/      
+ *          \-->--/       \-->--/        |          q1               q4 
+ *            q1             q4          |
+ *                Connected (C)          |                 Wing (W)
+ *
+ * C: trace(q1*adj(q2)*g5*gL[mu]*q3*adj(q4)*g5*gL[mu])
+ * W: trace(q1*adj(q2)*g5*gL[mu])*trace(q3*adj(q4)*g5*gL[mu])
+ * 
+ */
+
+/******************************************************************************
+ *                  TWeakHamiltonianNonEye implementation                     *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+TWeakHamiltonianNonEye::TWeakHamiltonianNonEye(const std::string name)
+: Module<WeakHamiltonianPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> TWeakHamiltonianNonEye::getInput(void)
+{
+    std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
+    
+    return in;
+}
+
+std::vector<std::string> TWeakHamiltonianNonEye::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void TWeakHamiltonianNonEye::setup(void)
+{
+
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void TWeakHamiltonianNonEye::execute(void)
+{
+    LOG(Message) << "Computing Weak Hamiltonian (Non-Eye type) contractions '" 
+                 << getName() << "' using quarks '" << par().q1 << "', '" 
+                 << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
+                 << "'." << std::endl;
+    
+    CorrWriter             writer(par().output);
+    PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
+    PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
+    PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
+    PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
+    Gamma g5            = Gamma(Gamma::Algebra::Gamma5);
+    LatticeComplex        expbuf(env().getGrid());
+    std::vector<TComplex> corrbuf;
+    std::vector<Result>   result(n_noneye_diag); 
+    unsigned int ndim   = env().getNd();
+
+    PropagatorField              tmp1(env().getGrid());
+    LatticeComplex               tmp2(env().getGrid());
+    std::vector<PropagatorField> C_i_side_loop(ndim, tmp1);
+    std::vector<PropagatorField> C_f_side_loop(ndim, tmp1);
+    std::vector<LatticeComplex>  W_i_side_loop(ndim, tmp2);
+    std::vector<LatticeComplex>  W_f_side_loop(ndim, tmp2);
+
+    // Setup for C-type contractions.
+    for (int mu = 0; mu < ndim; ++mu)
+    {
+        C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, GammaL(Gamma::gmu[mu]));
+        C_f_side_loop[mu] = MAKE_CW_SUBDIAG(q3, q4, GammaL(Gamma::gmu[mu]));
+    }
+
+    // Perform C-type contractions.    
+    SUM_MU(expbuf, trace(C_i_side_loop[mu]*C_f_side_loop[mu]))
+    MAKE_DIAG(expbuf, corrbuf, result[C_diag], "HW_C")
+
+    // Recycle sub-expressions for W-type contractions.
+    for (unsigned int mu = 0; mu < ndim; ++mu)
+    {
+        W_i_side_loop[mu] = trace(C_i_side_loop[mu]);
+        W_f_side_loop[mu] = trace(C_f_side_loop[mu]);
+    }
+
+    // Perform W-type contractions.
+    SUM_MU(expbuf, W_i_side_loop[mu]*W_f_side_loop[mu])
+    MAKE_DIAG(expbuf, corrbuf, result[W_diag], "HW_W")
+
+    write(writer, "HW_NonEye", result);
+}

extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp

@@ -0,0 +1,57 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_WeakHamiltonianNonEye_hpp_
+#define Hadrons_WeakHamiltonianNonEye_hpp_
+
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         WeakHamiltonianNonEye                              *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+enum
+{
+    W_diag = 0,
+    C_diag = 1,
+    n_noneye_diag = 2
+};
+
+// Wing and Connected subdiagram contractions
+#define MAKE_CW_SUBDIAG(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma)
+
+MAKE_WEAK_MODULE(WeakHamiltonianNonEye)
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WeakHamiltonianNonEye_hpp_

extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc

@@ -0,0 +1,135 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+/*
+ * Weak Hamiltonian + current contractions, disconnected topology for neutral 
+ * mesons.
+ * 
+ * These contractions are generated by operators Q_1,...,10 of the dS=1 Weak
+ * Hamiltonian in the physical basis and an additional current J (see e.g. 
+ * Fig 11 of arXiv:1507.03094).
+ * 
+ * Schematic:
+ *                        
+ *           q2          q4             q3
+ *       /--<--¬     /---<--¬       /---<--¬
+ *     /         \ /         \     /        \
+ *  i *           * H_W      |  J *          * f
+ *     \         / \         /     \        /
+ *      \--->---/   \-------/       \------/
+ *          q1 
+ * 
+ * options
+ * - q1: input propagator 1 (string)
+ * - q2: input propagator 2 (string)
+ * - q3: input propagator 3 (string), assumed to be sequential propagator 
+ * - q4: input propagator 4 (string), assumed to be a loop
+ * 
+ * type 1: trace(q1*adj(q2)*g5*gL[mu])*trace(loop*gL[mu])*trace(q3*g5)
+ * type 2: trace(q1*adj(q2)*g5*gL[mu]*loop*gL[mu])*trace(q3*g5)
+ */
+
+/*******************************************************************************
+ *                  TWeakNeutral4ptDisc implementation                         *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+TWeakNeutral4ptDisc::TWeakNeutral4ptDisc(const std::string name)
+: Module<WeakHamiltonianPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+std::vector<std::string> TWeakNeutral4ptDisc::getInput(void)
+{
+    std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
+    
+    return in;
+}
+
+std::vector<std::string> TWeakNeutral4ptDisc::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+void TWeakNeutral4ptDisc::setup(void)
+{
+
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+void TWeakNeutral4ptDisc::execute(void)
+{
+    LOG(Message) << "Computing Weak Hamiltonian neutral disconnected contractions '" 
+                 << getName() << "' using quarks '" << par().q1 << "', '" 
+                 << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
+                 << "'." << std::endl;
+
+    CorrWriter             writer(par().output);
+    PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
+    PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
+    PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
+    PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
+    Gamma g5            = Gamma(Gamma::Algebra::Gamma5);
+    LatticeComplex        expbuf(env().getGrid());
+    std::vector<TComplex> corrbuf;
+    std::vector<Result>   result(n_neut_disc_diag);
+    unsigned int ndim   = env().getNd();
+
+    PropagatorField              tmp(env().getGrid());
+    std::vector<PropagatorField> meson(ndim, tmp);
+    std::vector<PropagatorField> loop(ndim, tmp);
+    LatticeComplex               curr(env().getGrid());
+
+    // Setup for type 1 contractions.
+    for (int mu = 0; mu < ndim; ++mu)
+    {
+        meson[mu] = MAKE_DISC_MESON(q1, q2, GammaL(Gamma::gmu[mu]));
+        loop[mu] = MAKE_DISC_LOOP(q4, GammaL(Gamma::gmu[mu]));
+    }
+    curr = MAKE_DISC_CURR(q3, GammaL(Gamma::Algebra::Gamma5));
+
+    // Perform type 1 contractions.    
+    SUM_MU(expbuf, trace(meson[mu]*loop[mu]))
+    expbuf *= curr;
+    MAKE_DIAG(expbuf, corrbuf, result[neut_disc_1_diag], "HW_disc0_1")
+
+    // Perform type 2 contractions.
+    SUM_MU(expbuf, trace(meson[mu])*trace(loop[mu]))
+    expbuf *= curr;
+    MAKE_DIAG(expbuf, corrbuf, result[neut_disc_2_diag], "HW_disc0_2")
+
+    write(writer, "HW_disc0", result);
+}

extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp

@@ -0,0 +1,59 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson    <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_WeakNeutral4ptDisc_hpp_
+#define Hadrons_WeakNeutral4ptDisc_hpp_
+
+#include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                         WeakNeutral4ptDisc                                 *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MContraction)
+
+enum
+{
+    neut_disc_1_diag = 0,
+    neut_disc_2_diag = 1,
+    n_neut_disc_diag = 2
+};
+
+// Neutral 4pt disconnected subdiagram contractions.
+#define MAKE_DISC_MESON(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma)
+#define MAKE_DISC_LOOP(Q_LOOP, gamma) (Q_LOOP*gamma)
+#define MAKE_DISC_CURR(Q_c, gamma) (trace(Q_c*gamma))
+
+MAKE_WEAK_MODULE(WeakNeutral4ptDisc)
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WeakNeutral4ptDisc_hpp_

extras/Hadrons/Modules/MLoop/NoiseLoop.hpp

@@ -0,0 +1,132 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MLoop/NoiseLoop.hpp
+
+Copyright (C) 2016
+
+Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_NoiseLoop_hpp_
+#define Hadrons_NoiseLoop_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+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                                          *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MLoop)
+
+class NoiseLoopPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(NoiseLoopPar,
+                                    std::string, q,
+                                    std::string, eta);
+};
+
+template <typename FImpl>
+class TNoiseLoop: public Module<NoiseLoopPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TNoiseLoop(const std::string name);
+    // destructor
+    virtual ~TNoiseLoop(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(NoiseLoop, TNoiseLoop<FIMPL>, MLoop);
+
+/******************************************************************************
+ *                 TNoiseLoop implementation                                  *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TNoiseLoop<FImpl>::TNoiseLoop(const std::string name)
+: Module<NoiseLoopPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TNoiseLoop<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q, par().eta};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TNoiseLoop<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TNoiseLoop<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TNoiseLoop<FImpl>::execute(void)
+{
+    PropagatorField &loop = *env().template createLattice<PropagatorField>(getName());
+    PropagatorField &q    = *env().template getObject<PropagatorField>(par().q);
+    PropagatorField &eta  = *env().template getObject<PropagatorField>(par().eta);
+    loop = q*adj(eta);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_NoiseLoop_hpp_

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MSource/SeqGamma.hpp
 
 Copyright (C) 2015
 Copyright (C) 2016
+Copyright (C) 2017
 
 Author: Antonin Portelli <antonin.portelli@me.com>
 
@@ -149,9 +150,9 @@ void TSeqGamma<FImpl>::execute(void)
     for(unsigned int mu = 0; mu < env().getNd(); mu++)
     {
         LatticeCoordinate(coor, mu);
-        ph = ph + p[mu]*coor;
+        ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
     }
-    ph = exp(i*ph);
+    ph = exp((Real)(2*M_PI)*i*ph);
     LatticeCoordinate(t, Tp);
     src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q);
 }

extras/Hadrons/Modules/MSource/Wall.hpp

@@ -0,0 +1,147 @@
+/*************************************************************************************
+
+Grid physics library, www.github.com/paboyle/Grid 
+
+Source file: extras/Hadrons/Modules/MSource/Wall.hpp
+
+Copyright (C) 2017
+
+Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+This program is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 2 of the License, or
+(at your option) any later version.
+
+This program is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along
+with this program; if not, write to the Free Software Foundation, Inc.,
+51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+See the full license in the file "LICENSE" in the top level distribution directory
+*************************************************************************************/
+/*  END LEGAL */
+
+#ifndef Hadrons_WallSource_hpp_
+#define Hadrons_WallSource_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/*
+ 
+ Wall source
+ -----------------------------
+ * src_x = delta(x_3 - tW) * exp(i x.mom)
+ 
+ * options:
+ - tW: source timeslice (integer)
+ - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.")
+ 
+ */
+
+/******************************************************************************
+ *                         Wall                                               *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSource)
+
+class WallPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(WallPar,
+                                    unsigned int, tW,
+                                    std::string, mom);
+};
+
+template <typename FImpl>
+class TWall: public Module<WallPar>
+{
+public:
+    TYPE_ALIASES(FImpl,);
+public:
+    // constructor
+    TWall(const std::string name);
+    // destructor
+    virtual ~TWall(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Wall, TWall<FIMPL>, MSource);
+
+/******************************************************************************
+ *                 TWall implementation                                       *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TWall<FImpl>::TWall(const std::string name)
+: Module<WallPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TWall<FImpl>::getInput(void)
+{
+    std::vector<std::string> in;
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TWall<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWall<FImpl>::setup(void)
+{
+    env().template registerLattice<PropagatorField>(getName());
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TWall<FImpl>::execute(void)
+{    
+    LOG(Message) << "Generating wall source at t = " << par().tW 
+                 << " with momentum " << par().mom << std::endl;
+    
+    PropagatorField &src = *env().template createLattice<PropagatorField>(getName());
+    Lattice<iScalar<vInteger>> t(env().getGrid());
+    LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
+    std::vector<Real>          p;
+    Complex                    i(0.0,1.0);
+    
+    p  = strToVec<Real>(par().mom);
+    ph = zero;
+    for(unsigned int mu = 0; mu < Nd; mu++)
+    {
+        LatticeCoordinate(coor, mu);
+        ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
+    }
+    ph = exp((Real)(2*M_PI)*i*ph);
+    LatticeCoordinate(t, Tp);
+    src = 1.;
+    src = where((t == par().tW), src*ph, 0.*src);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WallSource_hpp_

extras/Hadrons/Modules/Quark.hpp

@@ -173,7 +173,7 @@ void TQuark<FImpl>::execute(void)
                 *env().template getObject<PropagatorField>(getName());
             
             axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0);
-            axpby_ssp_pplus(sol, 0., sol, 1., sol, 0, Ls_-1);
+            axpby_ssp_pplus(sol, 1., sol, 1., sol, 0, Ls_-1);
             ExtractSlice(tmp, sol, 0, 0);
             FermToProp(p4d, tmp, s, c);
         }

extras/Hadrons/modules.inc

@@ -1,4 +1,7 @@
 modules_cc =\
+  Modules/MContraction/WeakHamiltonianEye.cc \
+  Modules/MContraction/WeakHamiltonianNonEye.cc \
+  Modules/MContraction/WeakNeutral4ptDisc.cc \
   Modules/MGauge/Load.cc \
   Modules/MGauge/Random.cc \
   Modules/MGauge/Unit.cc
@@ -7,13 +10,21 @@ modules_hpp =\
   Modules/MAction/DWF.hpp \
   Modules/MAction/Wilson.hpp \
   Modules/MContraction/Baryon.hpp \
+  Modules/MContraction/DiscLoop.hpp \
+  Modules/MContraction/Gamma3pt.hpp \
   Modules/MContraction/Meson.hpp \
+  Modules/MContraction/WeakHamiltonian.hpp \
+  Modules/MContraction/WeakHamiltonianEye.hpp \
+  Modules/MContraction/WeakHamiltonianNonEye.hpp \
+  Modules/MContraction/WeakNeutral4ptDisc.hpp \
   Modules/MGauge/Load.hpp \
   Modules/MGauge/Random.hpp \
   Modules/MGauge/Unit.hpp \
+  Modules/MLoop/NoiseLoop.hpp \
   Modules/MSolver/RBPrecCG.hpp \
   Modules/MSource/Point.hpp \
   Modules/MSource/SeqGamma.hpp \
+  Modules/MSource/Wall.hpp \
   Modules/MSource/Z2.hpp \
   Modules/Quark.hpp
 

grid-config.in

@@ -0,0 +1,86 @@
+#! /bin/sh
+
+prefix=@prefix@
+exec_prefix=@exec_prefix@
+includedir=@includedir@
+
+usage()
+{
+  cat <<EOF
+Usage: grid-config [OPTION]
+
+Known values for OPTION are:
+
+  --prefix     show Grid installation prefix
+  --cxxflags   print pre-processor and compiler flags
+  --ldflags    print library linking flags
+  --libs       print library linking information
+  --summary    print full build summary
+  --help       display this help and exit
+  --version    output version information
+  --git        print git revision
+
+EOF
+  
+  exit $1
+}
+
+if test $# -eq 0; then
+  usage 1
+fi
+
+cflags=false
+libs=false
+
+while test $# -gt 0; do
+  case "$1" in
+    -*=*) optarg=`echo "$1" | sed 's/[-_a-zA-Z0-9]*=//'` ;;
+    *) optarg= ;;
+  esac
+  
+  case "$1" in
+    --prefix)
+      echo $prefix
+    ;;
+    
+    --version)
+      echo @VERSION@
+      exit 0
+    ;;
+    
+    --git)
+      echo "@GRID_BRANCH@ @GRID_SHA@"
+      exit 0
+    ;;
+    
+    --help)
+      usage 0
+    ;;
+    
+    --cxxflags)
+      echo @GRID_CXXFLAGS@
+    ;;
+    
+    --ldflags)
+      echo @GRID_LDFLAGS@
+    ;;
+    
+    --libs)
+      echo @GRID_LIBS@
+    ;;
+    
+    --summary)
+      echo ""
+      echo "@GRID_SUMMARY@"
+      echo ""
+    ;;
+    
+    *)
+      usage
+      exit 1
+    ;;
+  esac
+  shift
+done
+
+exit 0

lib/GridCore.h

@@ -38,28 +38,7 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
 #ifndef GRID_BASE_H
 #define GRID_BASE_H
 
-///////////////////
+#include <Grid/GridStd.h>
-// Std C++ dependencies
-///////////////////
-#include <cassert>
-#include <complex>
-#include <vector>
-#include <iostream>
-#include <iomanip>
-#include <random>
-#include <functional>
-#include <stdio.h>
-#include <stdlib.h>
-#include <stdio.h>
-#include <signal.h>
-#include <ctime>
-#include <sys/time.h>
-#include <chrono>
-
-///////////////////
-// Grid headers
-///////////////////
-#include "Config.h"
 
 #include <Grid/perfmon/Timer.h>
 #include <Grid/perfmon/PerfCount.h>

lib/GridStd.h

@@ -0,0 +1,27 @@
+#ifndef GRID_STD_H
+#define GRID_STD_H
+
+///////////////////
+// Std C++ dependencies
+///////////////////
+#include <cassert>
+#include <complex>
+#include <vector>
+#include <iostream>
+#include <iomanip>
+#include <random>
+#include <functional>
+#include <stdio.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <signal.h>
+#include <ctime>
+#include <sys/time.h>
+#include <chrono>
+
+///////////////////
+// Grid config
+///////////////////
+#include "Config.h"
+
+#endif /* GRID_STD_H */

lib/algorithms/approx/Remez.h

@@ -16,7 +16,7 @@
 #define INCLUDED_ALG_REMEZ_H
 
 #include <stddef.h>
-#include <Config.h>
+#include <Grid/GridStd.h>
 
 #ifdef HAVE_LIBGMP
 #include "bigfloat.h"

lib/qcd/action/fermion/CayleyFermion5D.cc

@@ -320,7 +320,7 @@ void CayleyFermion5D<Impl>::MDeriv  (GaugeField &mat,const FermionField &U,const
     this->DhopDeriv(mat,U,Din,dag);
   } else {
     //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
-    MeooeDag5D(U,Din);
+    Meooe5D(U,Din);
     this->DhopDeriv(mat,Din,V,dag);
   }
 };
@@ -335,8 +335,8 @@ void CayleyFermion5D<Impl>::MoeDeriv(GaugeField &mat,const FermionField &U,const
     this->DhopDerivOE(mat,U,Din,dag);
   } else {
     //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
-      MeooeDag5D(U,Din);
+    Meooe5D(U,Din);
-      this->DhopDerivOE(mat,Din,V,dag);
+    this->DhopDerivOE(mat,Din,V,dag);
   }
 };
 template<class Impl>
@@ -350,7 +350,7 @@ void CayleyFermion5D<Impl>::MeoDeriv(GaugeField &mat,const FermionField &U,const
     this->DhopDerivEO(mat,U,Din,dag);
   } else {
     //      U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call
-    MeooeDag5D(U,Din);
+    Meooe5D(U,Din);
     this->DhopDerivEO(mat,Din,V,dag);
   }
 };
@@ -380,6 +380,8 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Co
   ///////////////////////////////////////////////////////////
   // The Cayley coeffs (unprec)
   ///////////////////////////////////////////////////////////
+  assert(gamma.size()==Ls);
+
   omega.resize(Ls);
   bs.resize(Ls);
   cs.resize(Ls);
@@ -412,10 +414,11 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Co
   for(int i=0; i < Ls; i++){
     as[i] = 1.0;
     omega[i] = gamma[i]*zolo_hi; //NB reciprocal relative to Chroma NEF code
+    //    assert(fabs(omega[i])>0.0);
     bs[i] = 0.5*(bpc/omega[i] + bmc);
     cs[i] = 0.5*(bpc/omega[i] - bmc);
   }
-  
+
   ////////////////////////////////////////////////////////
   // Constants for the preconditioned matrix Cayley form
   ////////////////////////////////////////////////////////
@@ -425,12 +428,12 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Co
   ceo.resize(Ls);
   
   for(int i=0;i<Ls;i++){
-    bee[i]=as[i]*(bs[i]*(4.0-this->M5) +1.0);
+    bee[i]=as[i]*(bs[i]*(4.0-this->M5) +1.0);     
+    //    assert(fabs(bee[i])>0.0);
     cee[i]=as[i]*(1.0-cs[i]*(4.0-this->M5));
     beo[i]=as[i]*bs[i];
     ceo[i]=-as[i]*cs[i];
   }
-  
   aee.resize(Ls);
   aeo.resize(Ls);
   for(int i=0;i<Ls;i++){
@@ -474,14 +477,16 @@ void CayleyFermion5D<Impl>::SetCoefficientsInternal(RealD zolo_hi,std::vector<Co
 	
   { 
     Coeff_t delta_d=mass*cee[Ls-1];
-    for(int j=0;j<Ls-1;j++) delta_d *= cee[j]/bee[j];
+    for(int j=0;j<Ls-1;j++) {
+      //      assert(fabs(bee[j])>0.0);
+      delta_d *= cee[j]/bee[j];
+    }
     dee[Ls-1] += delta_d;
   }  
 
   int inv=1;
   this->MooeeInternalCompute(0,inv,MatpInv,MatmInv);
   this->MooeeInternalCompute(1,inv,MatpInvDag,MatmInvDag);
-
 }
 
 
@@ -495,7 +500,9 @@ void CayleyFermion5D<Impl>::MooeeInternalCompute(int dag, int inv,
   GridBase *grid = this->FermionRedBlackGrid();
   int LLs = grid->_rdimensions[0];
 
-  if ( LLs == Ls ) return; // Not vectorised in 5th direction
+  if ( LLs == Ls ) {
+    return; // Not vectorised in 5th direction
+  }
 
   Eigen::MatrixXcd Pplus  = Eigen::MatrixXcd::Zero(Ls,Ls);
   Eigen::MatrixXcd Pminus = Eigen::MatrixXcd::Zero(Ls,Ls);

lib/qcd/action/fermion/DomainWallFermion.h

@@ -68,7 +68,7 @@ namespace Grid {
 	Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham
 	assert(zdata->n==this->Ls);
 	
-	//	std::cout<<GridLogMessage << "DomainWallFermion with Ls="<<this->Ls<<std::endl;
+	std::cout<<GridLogMessage << "DomainWallFermion with Ls="<<this->Ls<<std::endl;
 	// Call base setter
 	this->SetCoefficientsTanh(zdata,1.0,0.0);
 

lib/serialisation/BaseIO.h

@@ -138,6 +138,54 @@ namespace Grid {
     unsigned int               dimInd_{0};
   };
   
+  // Pair IO utilities /////////////////////////////////////////////////////////
+  // helper function to parse input in the format "<obj1 obj2>"
+  template <typename T1, typename T2>
+  inline std::istream & operator>>(std::istream &is, std::pair<T1, T2> &buf)
+  {
+    T1 buf1;
+    T2 buf2;
+    char c;
+
+    // Search for "pair" delimiters.
+    do
+    {
+      is.get(c);
+    } while (c != '<' && !is.eof());
+    if (c == '<')
+    {
+      int start = is.tellg();
+      do
+      {
+        is.get(c);
+      } while (c != '>' && !is.eof());
+      if (c == '>')
+      {
+        int end = is.tellg();
+        int psize = end - start - 1;
+
+        // Only read data between pair limiters.
+        is.seekg(start);
+        std::string tmpstr(psize, ' ');
+        is.read(&tmpstr[0], psize);
+        std::istringstream temp(tmpstr);
+        temp >> buf1 >> buf2;
+        buf = std::make_pair(buf1, buf2);
+        is.seekg(end);
+      }
+    }
+    is.peek();
+    return is;
+  }
+  
+  // output to streams for pairs
+  template <class T1, class T2>
+  inline std::ostream & operator<<(std::ostream &os, const std::pair<T1, T2> &p)
+  {
+    os << "<" << p.first << " " << p.second << ">";
+    return os;
+  }
+
   // Abstract writer/reader classes ////////////////////////////////////////////
   // static polymorphism implemented using CRTP idiom
   class Serializable;

scripts/configure.commit

@@ -16,4 +16,4 @@ else # configure.ac is modified
   printf 'commit: %s-dirty\n' `git rev-parse --short HEAD`
   printf 'branch: %s\n' `git rev-parse --abbrev-ref HEAD`
   printf 'date  : %s\n' `git log -1 --date=short --pretty=format:%cd`
- fi
+fi

tests/IO/Test_serialisation.cc

@@ -115,6 +115,7 @@ int main(int argc,char **argv)
   // test serializable class writing
   myclass              obj(1234); // non-trivial constructor
   std::vector<myclass> vec;
+  std::pair<myenum, myenum> pair;
   
   std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl;
   write(WR,"obj",obj);
@@ -122,6 +123,8 @@ int main(int argc,char **argv)
   vec.push_back(myclass(1234));
   vec.push_back(myclass(5678));
   vec.push_back(myclass(3838));
+  pair = std::make_pair(myenum::red, myenum::blue);
+
   write(WR, "objvec", vec);
   std::cout << "-- serialisable class writing to std::cout:" << std::endl;
   std::cout << obj << std::endl;
@@ -129,21 +132,30 @@ int main(int argc,char **argv)
   std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl;
   std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl;
   
+  write(WR, "objpair", pair);
+  std::cout << "-- pair writing to std::cout:" << std::endl;
+  std::cout << pair << std::endl;
+  
   // read tests
   std::cout << "\n==== IO self-consistency tests" << std::endl;
   //// XML
   ioTest<XmlWriter, XmlReader>("iotest.xml", obj, "XML    (object)           ");
   ioTest<XmlWriter, XmlReader>("iotest.xml", vec, "XML    (vector of objects)");
+  ioTest<XmlWriter, XmlReader>("iotest.xml", pair, "XML    (pair of objects)");
   //// binary
   ioTest<BinaryWriter, BinaryReader>("iotest.bin", obj, "binary (object)           ");
   ioTest<BinaryWriter, BinaryReader>("iotest.bin", vec, "binary (vector of objects)");
+  ioTest<BinaryWriter, BinaryReader>("iotest.bin", pair, "binary (pair of objects)");
   //// text
   ioTest<TextWriter, TextReader>("iotest.dat", obj, "text   (object)           ");
   ioTest<TextWriter, TextReader>("iotest.dat", vec, "text   (vector of objects)");
+  ioTest<TextWriter, TextReader>("iotest.dat", pair, "text   (pair of objects)");
   //// HDF5
+#undef HAVE_HDF5
 #ifdef HAVE_HDF5
   ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", obj, "HDF5   (object)           ");
   ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", vec, "HDF5   (vector of objects)");
+  ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", pair, "HDF5   (pair of objects)");
 #endif
   
   std::cout << "\n==== vector flattening/reconstruction" << std::endl;

tests/core/Test_mobius_even_odd.cc

@@ -0,0 +1,287 @@
+    /*************************************************************************************
+
+    Grid physics library, www.github.com/paboyle/Grid 
+
+    Source file: ./tests/Test_dwf_even_odd.cc
+
+    Copyright (C) 2015
+
+Author: Peter Boyle <paboyle@ph.ed.ac.uk>
+Author: paboyle <paboyle@ph.ed.ac.uk>
+
+    This program is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License along
+    with this program; if not, write to the Free Software Foundation, Inc.,
+    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+    See the full license in the file "LICENSE" in the top level distribution directory
+    *************************************************************************************/
+    /*  END LEGAL */
+#include <Grid/Grid.h>
+
+using namespace std;
+using namespace Grid;
+using namespace Grid::QCD;
+
+template<class d>
+struct scal {
+  d internal;
+};
+
+  Gamma::Algebra Gmu [] = {
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+  };
+
+
+int main (int argc, char ** argv)
+{
+  Grid_init(&argc,&argv);
+
+  int threads = GridThread::GetThreads();
+  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
+
+
+  const int Ls=10;
+  GridCartesian         * UGrid   = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+  GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
+  GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
+  GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
+
+  std::vector<int> seeds4({1,2,3,4});
+  std::vector<int> seeds5({5,6,7,8});
+
+  GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
+  GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
+
+  LatticeFermion src   (FGrid); random(RNG5,src);
+  LatticeFermion phi   (FGrid); random(RNG5,phi);
+  LatticeFermion chi   (FGrid); random(RNG5,chi);
+  LatticeFermion result(FGrid); result=zero;
+  LatticeFermion    ref(FGrid);    ref=zero;
+  LatticeFermion    tmp(FGrid);    tmp=zero;
+  LatticeFermion    err(FGrid);    tmp=zero;
+  LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
+  std::vector<LatticeColourMatrix> U(4,UGrid);
+
+  // Only one non-zero (y)
+  Umu=zero;
+  for(int nn=0;nn<Nd;nn++){
+    random(RNG4,U[nn]);
+    if ( nn>0 ) 
+      U[nn]=zero;
+    PokeIndex<LorentzIndex>(Umu,U[nn],nn);
+  }
+
+  RealD mass=0.1;
+  RealD M5  =1.8;
+  std::vector < std::complex<double>  > omegas;
+#if 0
+  for(int i=0;i<Ls;i++){
+    double imag = 0.;
+    if (i==0) imag=1.;
+    if (i==Ls-1) imag=-1.;
+    std::complex<double> temp (0.25+0.01*i, imag*0.01);
+    omegas.push_back(temp);
+  }
+#else
+  omegas.push_back( std::complex<double>(1.45806438985048,-0) );
+  omegas.push_back( std::complex<double>(1.18231318389348,-0) );
+  omegas.push_back( std::complex<double>(0.830951166685955,-0) );
+  omegas.push_back( std::complex<double>(0.542352409156791,-0) );
+  omegas.push_back( std::complex<double>(0.341985020453729,-0) );
+  omegas.push_back( std::complex<double>(0.21137902619029,-0) );
+  omegas.push_back( std::complex<double>(0.126074299502912,-0) );
+  omegas.push_back( std::complex<double>(0.0990136651962626,-0) );
+  omegas.push_back( std::complex<double>(0.0686324988446592,0.0550658530827402) );
+  omegas.push_back( std::complex<double>(0.0686324988446592,-0.0550658530827402) );
+#endif
+
+  MobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, 0.5,0.5);
+  //  DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
+
+  LatticeFermion src_e (FrbGrid);
+  LatticeFermion src_o (FrbGrid);
+  LatticeFermion r_e   (FrbGrid);
+  LatticeFermion r_o   (FrbGrid);
+  LatticeFermion r_eo  (FGrid);
+  LatticeFermion r_eeoo(FGrid);
+
+  std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
+  std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
+
+  pickCheckerboard(Even,src_e,src);
+  pickCheckerboard(Odd,src_o,src);
+
+  Ddwf.Meooe(src_e,r_o);  std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
+  Ddwf.Meooe(src_o,r_e);  std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
+  setCheckerboard(r_eo,r_o);
+  setCheckerboard(r_eo,r_e);
+
+  Ddwf.Mooee(src_e,r_e);  std::cout<<GridLogMessage<<"Applied Mee"<<std::endl;
+  Ddwf.Mooee(src_o,r_o);  std::cout<<GridLogMessage<<"Applied Moo"<<std::endl;
+  setCheckerboard(r_eeoo,r_e);
+  setCheckerboard(r_eeoo,r_o);
+
+  r_eo=r_eo+r_eeoo;
+  Ddwf.M(src,ref);  
+
+  //  std::cout<<GridLogMessage << r_eo<<std::endl;
+  //  std::cout<<GridLogMessage << ref <<std::endl;
+
+  err= ref - r_eo;
+  std::cout<<GridLogMessage << "EO norm diff   "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
+    
+  LatticeComplex cerr(FGrid);
+  cerr = localInnerProduct(err,err);
+  //  std::cout<<GridLogMessage << cerr<<std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MooeeDagger is the dagger of Mooee by requiring                "<<std::endl;
+  std::cout<<GridLogMessage<<"=  < phi | Deo | chi > * = < chi | Deo^dag| phi>  "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  LatticeFermion chi_e   (FrbGrid);
+  LatticeFermion chi_o   (FrbGrid);
+
+  LatticeFermion dchi_e  (FrbGrid);
+  LatticeFermion dchi_o  (FrbGrid);
+
+  LatticeFermion phi_e   (FrbGrid);
+  LatticeFermion phi_o   (FrbGrid);
+
+  LatticeFermion dphi_e  (FrbGrid);
+  LatticeFermion dphi_o  (FrbGrid);
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+  pickCheckerboard(Even,phi_e,phi);
+  pickCheckerboard(Odd ,phi_o,phi);
+
+  Ddwf.Mooee(chi_e,dchi_o);
+  Ddwf.Mooee(chi_o,dchi_e);
+  Ddwf.MooeeDag(phi_e,dphi_o);
+  Ddwf.MooeeDag(phi_o,dphi_e);
+
+  ComplexD pDce = innerProduct(phi_e,dchi_e);
+  ComplexD pDco = innerProduct(phi_o,dchi_o);
+  ComplexD cDpe = innerProduct(chi_e,dphi_e);
+  ComplexD cDpo = innerProduct(chi_o,dphi_o);
+
+
+  std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
+  std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
+
+  std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
+  std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring                "<<std::endl;
+  std::cout<<GridLogMessage<<"=  < phi | Deo | chi > * = < chi | Deo^dag| phi>  "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+  pickCheckerboard(Even,phi_e,phi);
+  pickCheckerboard(Odd ,phi_o,phi);
+
+  Ddwf.Meooe(chi_e,dchi_o);
+  Ddwf.Meooe(chi_o,dchi_e);
+  Ddwf.MeooeDag(phi_e,dphi_o);
+  Ddwf.MeooeDag(phi_o,dphi_e);
+
+  pDce = innerProduct(phi_e,dchi_e);
+  pDco = innerProduct(phi_o,dchi_o);
+  cDpe = innerProduct(chi_e,dphi_e);
+  cDpo = innerProduct(chi_o,dphi_o);
+
+  std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
+  std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
+
+  std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
+  std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1                                         "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+
+  Ddwf.Mooee(chi_e,src_e);
+  Ddwf.MooeeInv(src_e,phi_e);
+
+  Ddwf.Mooee(chi_o,src_o);
+  Ddwf.MooeeInv(src_o,phi_o);
+  
+  setCheckerboard(phi,phi_e);
+  setCheckerboard(phi,phi_o);
+
+  err = phi-chi;
+  std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<< std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1                                   "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+
+  Ddwf.MooeeDag(chi_e,src_e);
+  Ddwf.MooeeInvDag(src_e,phi_e);
+
+  Ddwf.MooeeDag(chi_o,src_o);
+  Ddwf.MooeeInvDag(src_o,phi_o);
+  
+  setCheckerboard(phi,phi_e);
+  setCheckerboard(phi,phi_o);
+
+  err = phi-chi;
+  std::cout<<GridLogMessage << "norm diff   "<< norm2(err)<< std::endl;
+
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian              "<<std::endl;
+  std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
+  
+  random(RNG5,phi);
+  random(RNG5,chi);
+  pickCheckerboard(Even,chi_e,chi);
+  pickCheckerboard(Odd ,chi_o,chi);
+  pickCheckerboard(Even,phi_e,phi);
+  pickCheckerboard(Odd ,phi_o,phi);
+  RealD t1,t2;
+
+
+  SchurDiagMooeeOperator<MobiusFermionR,LatticeFermion> HermOpEO(Ddwf);
+  HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
+  HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
+
+  HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
+  HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
+
+  pDce = innerProduct(phi_e,dchi_e);
+  pDco = innerProduct(phi_o,dchi_o);
+  cDpe = innerProduct(chi_e,dphi_e);
+  cDpo = innerProduct(chi_o,dphi_o);
+
+  std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
+  std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
+
+  std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
+  std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
+  
+  Grid_finalize();
+}

tests/forces/Test_mobius_force.cc

@@ -46,19 +46,8 @@ int main (int argc, char ** argv)
   GridCartesian         * FGrid   = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
   GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
 
-  // Want a different conf at every run
+  std::vector<int> seeds4({1,2,3,4});
-  // First create an instance of an engine.
+  std::vector<int> seeds5({5,6,7,8});
-  std::random_device rnd_device;
-  // Specify the engine and distribution.
-  std::mt19937 mersenne_engine(rnd_device());
-  std::uniform_int_distribution<int> dist(1, 100);
-  
-  auto gen = std::bind(dist, mersenne_engine);
-  std::vector<int> seeds4(4);
-  generate(begin(seeds4), end(seeds4), gen);
-  
-  //std::vector<int> seeds4({1,2,3,5});
-  std::vector<int> seeds5({5,6,7,1});
   GridParallelRNG          RNG5(FGrid);  RNG5.SeedFixedIntegers(seeds5);
   GridParallelRNG          RNG4(UGrid);  RNG4.SeedFixedIntegers(seeds4);
 
@@ -76,15 +65,13 @@ int main (int argc, char ** argv)
   ////////////////////////////////////
   // Unmodified matrix element
   ////////////////////////////////////
-  RealD mass = 0.01; 
+  RealD mass=0.01; 
-  RealD M5   = 1.8;
+  RealD M5=1.8; 
-  RealD b    = 2.0;
+  RealD b=0.5;
-  RealD c    = 1.0;
+  RealD c=0.5;
   MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
   Ddwf.M   (phi,Mphi);
 
-
-  
   ComplexD S    = innerProduct(Mphi,Mphi); // pdag MdagM p
 
   // get the deriv of phidag MdagM phi with respect to "U"

tests/forces/Test_wilson_force_phiMdagMphi.cc

@@ -1,167 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_wilson_force_phiMdagMphi.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
-#include <Grid/Grid.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
- 
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
-
-  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
-
-  int threads = GridThread::GetThreads();
-  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
-
-  std::vector<int> seeds({1,2,3,4});
-
-  GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
-
-  LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
-  LatticeFermion Mphi       (&Grid); 
-  LatticeFermion Mdagphi       (&Grid); 
-  LatticeFermion MphiPrime  (&Grid); 
-  LatticeFermion MdagphiPrime  (&Grid); 
-  LatticeFermion dMphi      (&Grid); 
-
-  LatticeGaugeField U(&Grid);
-
-
-  SU3::HotConfiguration(pRNG,U);
-  //  SU3::ColdConfiguration(pRNG,U);
-  
-  ////////////////////////////////////
-  // Unmodified matrix element
-  ////////////////////////////////////
-  RealD mass=-4.0; //kills the diagonal term
-  WilsonFermionR Dw     (U,     Grid,RBGrid,mass);
-  Dw.M   (phi,Mphi);
-  Dw.Mdag(phi,Mdagphi);
-
-  ComplexD S    = innerProduct(Mphi,Mphi); // pdag MdagM p
-  ComplexD Sdag = innerProduct(Mdagphi,Mdagphi); // pdag MMdag p
-
-  // get the deriv of phidag MdagM phi with respect to "U"
-  LatticeGaugeField UdSdU(&Grid);
-  LatticeGaugeField UdSdUdag(&Grid);
-  LatticeGaugeField tmp(&Grid);
-
-  Dw.MDeriv(tmp , Mphi,  phi,DaggerNo );  UdSdU=tmp;
-
-  Dw.MDeriv(tmp , Mdagphi,  phi,DaggerYes );  UdSdUdag=tmp;
-
-
-  LatticeFermion dMdagphi      (&Grid);  dMdagphi=zero;
-  LatticeFermion Ftmp      (&Grid);
-
-
-  //  Dw.MDeriv(UdSdU,Mdagphi,  phi,DaggerYes );// UdSdU =UdSdU +tmp;
-
-  ////////////////////////////////////
-  // Modify the gauge field a little in one dir
-  ////////////////////////////////////
-  RealD dt = 1.0e-3;
-
-  LatticeColourMatrix mommu(&Grid); 
-  LatticeGaugeField mom(&Grid); 
-  LatticeGaugeField Uprime(&Grid); 
-
-  for(int mu=0;mu<Nd;mu++){
-
-    SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
-
-    //    Dw.DoubleStore(Dw.Umu,Uprime); // update U _and_ Udag
-    Dw.DhopDirDisp(phi,Ftmp,mu,mu+4,DaggerYes); 
-    dMdagphi=dMdagphi+mommu*Ftmp*dt;
-    
-    PokeIndex<LorentzIndex>(mom,mommu,mu);
-
-    parallel_for(auto i=mom.begin();i<mom.end();i++){
-      Uprime[i](mu) =U[i](mu)+ mom[i](mu)*U[i](mu)*dt;
-      Dw.Umu[i](mu) =Uprime[i](mu); // update U but _not_ Udag
-    }
-
-  }
-
-  Dw.Mdag(phi,MdagphiPrime);
-  Dw.M   (phi,MphiPrime);
-
-  std::cout << GridLogMessage << "deltaMdag phi    "<< norm2(dMdagphi) <<std::endl;
-  Ftmp=MdagphiPrime - Mdagphi;
-  std::cout << GridLogMessage << "diff Mdag phi    "<< norm2(Ftmp) <<std::endl;
-  Ftmp = Ftmp - dMdagphi;
-  std::cout << GridLogMessage << "err  Mdag phi    "<< norm2(Ftmp) <<std::endl;
-  std::cout << dMdagphi<<std::endl;
-  Ftmp=MdagphiPrime - Mdagphi;
-  std::cout << Ftmp<<std::endl;
-
-
-  ComplexD Sprime    = innerProduct(Mphi   ,MphiPrime);
-  ComplexD Sprimedag = innerProduct(Mdagphi,MdagphiPrime);
-
-  ComplexD deltaSdag = innerProduct(Mdagphi,dMdagphi);
-  std::cout << GridLogMessage << "deltaSdag from inner prod of mom* M[u]     "<<deltaSdag<<std::endl;
-
-  //////////////////////////////////////////////
-  // Use derivative to estimate dS
-  //////////////////////////////////////////////
-  LatticeComplex dS(&Grid); dS = zero;
-  LatticeComplex dSdag(&Grid); dSdag = zero;
-  parallel_for(auto i=mom.begin();i<mom.end();i++){
-    for(int mu=0;mu<Nd;mu++){
-      //      dS[i]() = dS[i]()+trace(mom[i](mu) * UdSdU[i](mu) - mom[i](mu)* adj( UdSdU[i](mu)) )*dt;
-      dS[i]()    =    dS[i]()+trace(mom[i](mu) * UdSdU[i](mu) )*dt;
-      dSdag[i]() = dSdag[i]()+trace(mom[i](mu) * UdSdUdag[i](mu) )*dt;
-    }
-  }
-  Complex dSpred    = sum(dS);
-  Complex dSdagpred = sum(dSdag);
-
-  std::cout << GridLogMessage << " S      "<<S<<std::endl;
-  std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
-  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
-  std::cout << GridLogMessage << "predict dS    "<< dSpred <<std::endl;
-  std::cout << "\n\n"<<std::endl;
-  std::cout << GridLogMessage << " Sdag      "<<Sdag<<std::endl;
-  std::cout << GridLogMessage << " Sprimedag "<<Sprimedag<<std::endl;
-  std::cout << GridLogMessage << "dSdag      "<<Sprimedag-Sdag<<std::endl;
-  std::cout << GridLogMessage << "predict dSdag    "<< dSdagpred <<std::endl;
-
-  std::cout<< GridLogMessage << "Done" <<std::endl;
-  Grid_finalize();
-}

tests/forces/Test_wilson_force_phiMphi.cc

@@ -1,189 +0,0 @@
-    /*************************************************************************************
-
-    Grid physics library, www.github.com/paboyle/Grid 
-
-    Source file: ./tests/Test_wilson_force_phiMphi.cc
-
-    Copyright (C) 2015
-
-Author: Peter Boyle <paboyle@ph.ed.ac.uk>
-
-    This program is free software; you can redistribute it and/or modify
-    it under the terms of the GNU General Public License as published by
-    the Free Software Foundation; either version 2 of the License, or
-    (at your option) any later version.
-
-    This program is distributed in the hope that it will be useful,
-    but WITHOUT ANY WARRANTY; without even the implied warranty of
-    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-    GNU General Public License for more details.
-
-    You should have received a copy of the GNU General Public License along
-    with this program; if not, write to the Free Software Foundation, Inc.,
-    51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-
-    See the full license in the file "LICENSE" in the top level distribution directory
-    *************************************************************************************/
-    /*  END LEGAL */
-#include <Grid/Grid.h>
-
-using namespace std;
-using namespace Grid;
-using namespace Grid::QCD;
-
- 
-
-int main (int argc, char ** argv)
-{
-  Grid_init(&argc,&argv);
-
-  std::vector<int> latt_size   = GridDefaultLatt();
-  std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
-  std::vector<int> mpi_layout  = GridDefaultMpi();
-
-  GridCartesian               Grid(latt_size,simd_layout,mpi_layout);
-  GridRedBlackCartesian     RBGrid(latt_size,simd_layout,mpi_layout);
-
-  int threads = GridThread::GetThreads();
-  std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
-
-  std::vector<int> seeds({1,2,3,4});
-
-  GridParallelRNG          pRNG(&Grid);
-  pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
-
-  LatticeFermion phi        (&Grid); gaussian(pRNG,phi);
-  LatticeFermion Mphi       (&Grid); 
-  LatticeFermion MphiPrime  (&Grid); 
-  LatticeFermion dMphi      (&Grid); 
-
-  LatticeGaugeField U(&Grid);
-
-
-  SU3::HotConfiguration(pRNG,U);
-
-  ////////////////////////////////////
-  // Unmodified matrix element
-  ////////////////////////////////////
-  RealD mass=-4.0; //kills the diagonal term
-  WilsonFermionR Dw     (U,     Grid,RBGrid,mass);
-  Dw.M(phi,Mphi);
-
-  ComplexD S = innerProduct(phi,Mphi);
-
-  // get the deriv
-  LatticeGaugeField UdSdU(&Grid);
-  Dw.MDeriv(UdSdU,phi, phi,DaggerNo ); 
-
-
-  ////////////////////////////////////
-  // Modify the gauge field a little in one dir
-  ////////////////////////////////////
-  RealD dt = 1.0e-3;
-  Complex Complex_i(0,1);
-
-  LatticeColourMatrix Umu(&Grid);
-  LatticeColourMatrix Umu_save(&Grid);
-  LatticeColourMatrix dU (&Grid);
-  LatticeColourMatrix mom(&Grid); 
-  SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mom); // Traceless antihermitian momentum; gaussian in lie alg
-
-
-  // check mom is as i expect
-  LatticeColourMatrix tmpmom(&Grid); 
-  tmpmom = mom+adj(mom);
-  std::cout << GridLogMessage << "mom anti-herm check "<< norm2(tmpmom)<<std::endl;
-  std::cout << GridLogMessage << "mom tr check "<< norm2(trace(mom))<<std::endl;
-  
-  const int mu=0;
-  Umu = PeekIndex<LorentzIndex>(U,mu);
-  Umu_save=Umu;
-  dU = mom * Umu * dt;
-  Umu= Umu+dU;
-  PokeIndex<LorentzIndex>(Dw.Umu,Umu,mu);
-
-  Dw.M(phi,MphiPrime);
-
-  ComplexD Sprime = innerProduct(phi,MphiPrime);
-
-  std::cout << GridLogMessage << " S      "<<S<<std::endl;
-  std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
-  std::cout << GridLogMessage << "dS      "<<Sprime-S<<std::endl;
-
-  Dw.Umu=zero;
-  PokeIndex<LorentzIndex>(Dw.Umu,dU,mu);
-  Dw.M(phi,dMphi);
-
-
-  ComplexD deltaS = innerProduct(phi,dMphi);
-  std::cout << GridLogMessage << "deltaS      "<<deltaS<<std::endl;
-
-  Dw.Umu=zero;
-  PokeIndex<LorentzIndex>(Dw.Umu,Umu_save,mu);
-  Dw.Mdir(phi,dMphi,mu,1);
-  dMphi = dt*mom*dMphi;
-
-  deltaS = innerProduct(phi,dMphi);
-  std::cout << GridLogMessage << "deltaS from inner prod of mom* M[u]     "<<deltaS<<std::endl;
-
-  deltaS = sum(trace(outerProduct(dMphi,phi)));
-  std::cout << GridLogMessage << "deltaS from trace outer prod of deltaM      "<<deltaS<<std::endl;
-
-/*
-  LatticeComplex lip(&Grid);
-  lip  = localInnerProduct(phi,dMphi);
-  
-  LatticeComplex trop(&Grid);
-  trop = trace(outerProduct(dMphi,phi));
-
-  LatticeSpinColourMatrix op(&Grid);
-  op = outerProduct(dMphi,phi);
-
-  LatticeSpinColourMatrix hop(&Grid);
-  LatticeComplex op_cpt(&Grid);
-  for(int s1=0;s1<Ns;s1++){
-  for(int s2=0;s2<Ns;s2++){
-  for(int c1=0;c1<Nc;c1++){
-  for(int c2=0;c2<Nc;c2++){
-
-    op_cpt = peekColour(peekSpin(dMphi,s1),c1) * adj(peekColour(peekSpin(phi,s2),c2));
-
-    parallel_for(auto i=hop.begin();i<hop.end();i++){
-      hop[i]()(s1,s2)(c1,c2) = op_cpt[i]()()();
-    }
-    
-  }}}}
-
-  LatticeSpinColourMatrix diffop(&Grid);
-
-  diffop = hop - op;
-  std::cout << GridLogMessage << "hand outer prod diff   "<<norm2(diffop)<<std::endl;
-
-  deltaS = sum(trace(hop));
-  std::cout << GridLogMessage << "deltaS hop   "<<deltaS<<std::endl;
-
-  std::cout << GridLogMessage<< "  phi[0] : "<<  phi._odata[0]<<std::endl;
-  std::cout << GridLogMessage<< "dMphi[0] : "<<dMphi._odata[0]<<std::endl;
-  std::cout << GridLogMessage<< "hop[0]   : "<<  hop._odata[0]<<std::endl;
-  std::cout << GridLogMessage<< " op[0]   : "<<   op._odata[0]<<std::endl;
-
-
-  std::cout << GridLogMessage << "lip      "<<lip<<std::endl;
-  std::cout << GridLogMessage << "trop     "<<trop<<std::endl;
-
-*/  
-  
-  //  std::cout << GridLogMessage << " UdSdU " << UdSdU << std::endl;
-
-  LatticeComplex dS(&Grid); dS = zero;
-  parallel_for(auto i=mom.begin();i<mom.end();i++){
-    dS[i]() = trace(mom[i]() * UdSdU[i](mu) )*dt;
-  }
-  Complex dSpred = sum(dS);
-
-  std::cout << GridLogMessage << "predict dS    "<< dSpred <<std::endl;
-
-
-  cout<< GridLogMessage << "Done" <<std::endl;
-  Grid_finalize();
-}

tests/hadrons/Test_hadrons.hpp

@@ -0,0 +1,368 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+
+ Source file: tests/hadrons/Test_hadrons.hpp
+
+ Copyright (C) 2017
+
+ Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include <Grid/Hadrons/Application.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+
+/*******************************************************************************
+ * Macros to reduce code duplication.
+ ******************************************************************************/
+// Useful definitions
+#define ZERO_MOM "0. 0. 0. 0."
+#define INIT_INDEX(s, n) (std::string(s) + "_" + std::to_string(n))
+#define ADD_INDEX(s, n) (s + "_" + std::to_string(n))
+#define LABEL_3PT(s, t1, t2) ADD_INDEX(INIT_INDEX(s, t1), t2)
+#define LABEL_4PT(s, t1, t2, t3) ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3)
+#define LABEL_4PT_NOISE(s, t1, t2, t3, nn) ADD_INDEX(ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3), nn)
+
+// Wall source/sink macros
+#define NAME_3MOM_WALL_SOURCE(t, mom) ("wall_" + std::to_string(t) + "_" + mom)
+#define NAME_WALL_SOURCE(t) NAME_3MOM_WALL_SOURCE(t, ZERO_MOM)
+#define NAME_POINT_SOURCE(pos) ("point_" + pos)
+
+#define MAKE_3MOM_WALL_PROP(tW, mom, propName, solver)\
+{\
+    std::string srcName = NAME_3MOM_WALL_SOURCE(tW, mom);\
+    makeWallSource(application, srcName, tW, mom);\
+    makePropagator(application, propName, srcName, solver);\
+}
+
+#define MAKE_WALL_PROP(tW, propName, solver)\
+        MAKE_3MOM_WALL_PROP(tW, ZERO_MOM, propName, solver)
+
+// Sequential source macros
+#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, propName, solver)\
+{\
+    std::string srcName = ADD_INDEX(qSrc + "_seq", tS);\
+    makeSequentialSource(application, srcName, qSrc, tS, mom);\
+    makePropagator(application, propName, srcName, solver);\
+}
+
+// Point source macros
+#define MAKE_POINT_PROP(pos, propName, solver)\
+{\
+    std::string srcName = NAME_POINT_SOURCE(pos);\
+    makePointSource(application, srcName, pos);\
+    makePropagator(application, propName, srcName, solver);\
+}
+
+/*******************************************************************************
+ * Functions for propagator construction.
+ ******************************************************************************/
+ 
+/*******************************************************************************
+ * Name: makePointSource
+ * Purpose: Construct point source and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             pos         - Position of point source.
+ * Returns: None.
+ ******************************************************************************/
+inline void makePointSource(Application &application, std::string srcName,
+                            std::string pos)
+{
+    // If the source already exists, don't make the module again.
+    if (!(Environment::getInstance().hasModule(srcName)))
+    {
+        MSource::Point::Par pointPar;
+        pointPar.position = pos;
+        application.createModule<MSource::Point>(srcName, pointPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeSequentialSource
+ * Purpose: Construct sequential source and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             qSrc        - Input quark for sequential inversion.
+ *             tS          - sequential source timeslice.
+ *             mom         - momentum insertion (default is zero).
+ * Returns: None.
+ ******************************************************************************/
+inline void makeSequentialSource(Application &application, std::string srcName,
+                                 std::string qSrc, unsigned int tS,
+                                 std::string mom = ZERO_MOM)
+{
+    // If the source already exists, don't make the module again.
+    if (!(Environment::getInstance().hasModule(srcName)))
+    {
+        MSource::SeqGamma::Par seqPar;
+        seqPar.q   = qSrc;
+        seqPar.tA  = tS;
+        seqPar.tB  = tS;
+        seqPar.mom = mom;
+        seqPar.gamma = Gamma::Algebra::GammaT;
+        application.createModule<MSource::SeqGamma>(srcName, seqPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeWallSource
+ * Purpose: Construct wall source and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             tW          - wall source timeslice.
+ *             mom         - momentum insertion (default is zero).
+ * Returns: None.
+ ******************************************************************************/
+inline void makeWallSource(Application &application, std::string srcName,
+                           unsigned int tW, std::string mom = ZERO_MOM)
+{
+    // If the source already exists, don't make the module again.
+    if (!(Environment::getInstance().hasModule(srcName)))
+    {
+        MSource::Wall::Par wallPar;
+        wallPar.tW  = tW;
+        wallPar.mom = mom;
+        application.createModule<MSource::Wall>(srcName, wallPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeWallSink
+ * Purpose: Wall sink smearing of a propagator.
+ * Parameters: application - main application that stores modules.
+ *             propName    - name of input propagator.
+ *             wallName    - name of smeared propagator.
+ *             mom         - momentum insertion (default is zero).
+ * Returns: None.
+ ******************************************************************************/
+inline void makeWallSink(Application &application, std::string propName,
+                         std::string wallName, std::string mom = ZERO_MOM)
+{
+    // If the propagator has already been smeared, don't smear it again.
+    // Temporarily removed, strategy for sink smearing likely to change.
+    /*if (!(Environment::getInstance().hasModule(wallName)))
+    {
+        MSink::Wall::Par wallPar;
+        wallPar.q   = propName;
+        wallPar.mom = mom;
+        application.createModule<MSink::Wall>(wallName, wallPar);
+    }*/
+}
+
+/*******************************************************************************
+ * Name: makePropagator
+ * Purpose: Construct source and propagator then add to application module.
+ * Parameters: application - main application that stores modules.
+ *             propName    - name of propagator module to create.
+ *             srcName     - name of source module to use.
+ *             solver      - solver to use (default is CG).
+ * Returns: None.
+ ******************************************************************************/
+inline void makePropagator(Application &application, std::string &propName,
+                           std::string &srcName, std::string &solver)
+{
+    // If the propagator already exists, don't make the module again.
+    if (!(Environment::getInstance().hasModule(propName)))
+    {
+        Quark::Par         quarkPar;
+        quarkPar.source = srcName;
+        quarkPar.solver = solver;
+        application.createModule<Quark>(propName, quarkPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: makeLoop
+ * Purpose: Use noise source and inversion result to make loop propagator, then 
+ *          add to application module.
+ * Parameters: application - main application that stores modules.
+ *             propName    - name of propagator module to create.
+ *             srcName     - name of noise source module to use.
+ *             resName     - name of inversion result on given noise source.
+ * Returns: None.
+ ******************************************************************************/
+inline void makeLoop(Application &application, std::string &propName,
+                     std::string &srcName, std::string &resName)
+{
+    // If the loop propagator already exists, don't make the module again.
+    if (!(Environment::getInstance().hasModule(propName)))
+    {
+        MLoop::NoiseLoop::Par loopPar;
+        loopPar.q   = resName;
+        loopPar.eta = srcName;
+        application.createModule<MLoop::NoiseLoop>(propName, loopPar);
+    }
+}
+
+/*******************************************************************************
+ * Contraction module creation.
+ ******************************************************************************/
+
+/*******************************************************************************
+ * Name: mesonContraction
+ * Purpose: Create meson contraction module and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             npt         - specify n-point correlator (for labelling).
+ *             q1          - quark propagator 1.
+ *             q2          - quark propagator 2.
+ *             label       - unique label to construct module name.
+ *             mom         - momentum to project (default is zero)
+ *             gammas      - gamma insertions at source and sink.
+ * Returns: None.
+ ******************************************************************************/
+inline void mesonContraction(Application &application, unsigned int npt, 
+                             std::string &q1, std::string &q2,
+                             std::string &label, 
+                             std::string mom = ZERO_MOM,
+                             std::string gammas = "<Gamma5 Gamma5>")
+{
+    std::string modName = std::to_string(npt) + "pt_" + label;
+    if (!(Environment::getInstance().hasModule(modName)))
+    {
+        MContraction::Meson::Par mesPar;
+        mesPar.output = std::to_string(npt) + "pt/" + label;
+        mesPar.q1 = q1;
+        mesPar.q2 = q2;
+        mesPar.mom = mom;
+        mesPar.gammas = gammas;
+        application.createModule<MContraction::Meson>(modName, mesPar);
+    }
+ }
+
+/*******************************************************************************
+ * Name: gamma3ptContraction
+ * Purpose: Create gamma3pt contraction module and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             npt         - specify n-point correlator (for labelling).
+ *             q1          - quark propagator 1.
+ *             q2          - quark propagator 2.
+ *             q3          - quark propagator 3.
+ *             label       - unique label to construct module name.
+ *             gamma       - gamma insertions between q2 and q3.
+ * Returns: None.
+ ******************************************************************************/
+inline void gamma3ptContraction(Application &application, unsigned int npt, 
+                                std::string &q1, std::string &q2,
+                                std::string &q3, std::string &label, 
+                                Gamma::Algebra gamma = Gamma::Algebra::Identity)
+{
+    std::string modName = std::to_string(npt) + "pt_" + label;
+    if (!(Environment::getInstance().hasModule(modName)))
+    {
+        MContraction::Gamma3pt::Par gamma3ptPar;
+        gamma3ptPar.output = std::to_string(npt) + "pt/" + label;
+        gamma3ptPar.q1 = q1;
+        gamma3ptPar.q2 = q2;
+        gamma3ptPar.q3 = q3;
+        gamma3ptPar.gamma = gamma;
+        application.createModule<MContraction::Gamma3pt>(modName, gamma3ptPar);
+    }
+ }
+
+/*******************************************************************************
+ * Name: weakContraction[Eye,NonEye]
+ * Purpose: Create Weak Hamiltonian contraction module for Eye/NonEye topology
+ *          and add to application module.
+ * Parameters: application - main application that stores modules.
+ *             npt         - specify n-point correlator (for labelling).
+ *             q1          - quark propagator 1.
+ *             q2          - quark propagator 2.
+ *             q3          - quark propagator 3.
+ *             q4          - quark propagator 4.
+ *             label       - unique label to construct module name.
+ * Returns: None.
+ ******************************************************************************/
+#define HW_CONTRACTION(top) \
+inline void weakContraction##top(Application &application, unsigned int npt,\
+                                 std::string &q1, std::string &q2, \
+                                 std::string &q3, std::string &q4, \
+                                 std::string &label)\
+{\
+    std::string modName = std::to_string(npt) + "pt_" + label;\
+    if (!(Environment::getInstance().hasModule(modName)))\
+    {\
+        MContraction::WeakHamiltonian##top::Par weakPar;\
+        weakPar.output = std::to_string(npt) + "pt/" + label;\
+        weakPar.q1 = q1;\
+        weakPar.q2 = q2;\
+        weakPar.q3 = q3;\
+        weakPar.q4 = q4;\
+        application.createModule<MContraction::WeakHamiltonian##top>(modName, weakPar);\
+    }\
+}
+HW_CONTRACTION(Eye)    // weakContractionEye
+HW_CONTRACTION(NonEye) // weakContractionNonEye
+
+/*******************************************************************************
+ * Name: disc0Contraction
+ * Purpose: Create contraction module for 4pt Weak Hamiltonian + current
+ *          disconnected topology for neutral mesons and add to application 
+ *          module.
+ * Parameters: application - main application that stores modules.
+ *             q1          - quark propagator 1.
+ *             q2          - quark propagator 2.
+ *             q3          - quark propagator 3.
+ *             q4          - quark propagator 4.
+ *             label       - unique label to construct module name.
+ * Returns: None.
+ ******************************************************************************/
+inline void disc0Contraction(Application &application, 
+                             std::string &q1, std::string &q2,
+                             std::string &q3, std::string &q4,
+                             std::string &label)
+{
+    std::string modName = "4pt_" + label;
+    if (!(Environment::getInstance().hasModule(modName)))
+    {
+        MContraction::WeakNeutral4ptDisc::Par disc0Par;
+        disc0Par.output = "4pt/" + label;
+        disc0Par.q1 = q1;
+        disc0Par.q2 = q2;
+        disc0Par.q3 = q3;
+        disc0Par.q4 = q4;
+        application.createModule<MContraction::WeakNeutral4ptDisc>(modName, disc0Par);
+    }
+ }
+
+/*******************************************************************************
+ * Name: discLoopContraction
+ * Purpose: Create contraction module for disconnected loop and add to
+ *          application module.
+ * Parameters: application - main application that stores modules.
+ *             q_loop      - loop quark propagator.
+ *             modName     - unique module name.
+ *             gamma       - gamma matrix to use in contraction.
+ * Returns: None.
+ ******************************************************************************/
+inline void discLoopContraction(Application &application,
+                                std::string &q_loop, std::string &modName,
+                                Gamma::Algebra gamma = Gamma::Algebra::Identity)
+{
+    if (!(Environment::getInstance().hasModule(modName)))
+    {
+        MContraction::DiscLoop::Par discPar;
+        discPar.output = "disc/" + modName;
+        discPar.q_loop = q_loop;
+        discPar.gamma  = gamma;
+        application.createModule<MContraction::DiscLoop>(modName, discPar);
+    }
+ }

tests/hadrons/Test_hadrons_meson_3pt.cc

@@ -30,14 +30,6 @@
 using namespace Grid;
 using namespace Hadrons;
 
-static Gamma::Algebra gmu[4] =
-{
-    Gamma::Algebra::GammaX,
-    Gamma::Algebra::GammaY,
-    Gamma::Algebra::GammaZ,
-    Gamma::Algebra::GammaT
-};
-
 int main(int argc, char *argv[])
 {
     // initialization //////////////////////////////////////////////////////////
@@ -110,7 +102,7 @@ int main(int argc, char *argv[])
             seqName.push_back(std::vector<std::string>(Nd));
             for (unsigned int mu = 0; mu < Nd; ++mu)
             {
-                seqPar.gamma   = gmu[mu];
+                seqPar.gamma   = 0x1 << mu;
                 seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
                                  + "_" + std::to_string(seqPar.tA) + "-"
                                  + qName[i];
@@ -135,11 +127,11 @@ int main(int argc, char *argv[])
         for (unsigned int i = 0; i < flavour.size(); ++i)
         for (unsigned int j = i; j < flavour.size(); ++j)
         {
-            mesPar.output      = "mesons/Z2_" + flavour[i] + flavour[j];
+            mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
-            mesPar.q1          = qName[i];
+            mesPar.q1     = qName[i];
-            mesPar.q2          = qName[j];
+            mesPar.q2     = qName[j];
-            mesPar.gammaSource = Gamma::Algebra::Gamma5;
+            mesPar.gammas = "all";
-            mesPar.gammaSink   = Gamma::Algebra::Gamma5;
+            mesPar.mom    = "0. 0. 0. 0.";
             application.createModule<MContraction::Meson>("meson_Z2_"
                                                           + std::to_string(t)
                                                           + "_"
@@ -157,6 +149,8 @@ int main(int argc, char *argv[])
                             + std::to_string(mu);
             mesPar.q1     = qName[i];
             mesPar.q2     = seqName[j][mu];
+            mesPar.gammas = "all";
+            mesPar.mom    = "0. 0. 0. 0.";
             application.createModule<MContraction::Meson>("3pt_Z2_"
                                                           + std::to_string(t)
                                                           + "_"

tests/hadrons/Test_hadrons_rarekaon.cc

@@ -0,0 +1,337 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+
+ Source file: tests/hadrons/Test_hadrons_rarekaon.cc
+
+ Copyright (C) 2017
+
+ Author: Andrew Lawson <andrew.lawson1991@gmail.com>
+
+ This program is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 2 of the License, or
+ (at your option) any later version.
+
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License along
+ with this program; if not, write to the Free Software Foundation, Inc.,
+ 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+
+ See the full license in the file "LICENSE" in the top level distribution
+ directory.
+ *******************************************************************************/
+
+#include "Test_hadrons.hpp"
+
+using namespace Grid;
+using namespace Hadrons;
+
+enum quarks
+{
+   light   = 0,
+   strange = 1,
+   charm   = 2  
+};
+
+int main(int argc, char *argv[])
+{
+    // parse command line //////////////////////////////////////////////////////
+    std::string configStem;
+    
+    if (argc < 2)
+    {
+        std::cerr << "usage: " << argv[0] << " <configuration filestem> [Grid options]";
+        std::cerr << std::endl;
+        std::exit(EXIT_FAILURE);
+    }
+    configStem = argv[1];
+    
+    // initialization //////////////////////////////////////////////////////////
+    Grid_init(&argc, &argv);
+    HadronsLogError.Active(GridLogError.isActive());
+    HadronsLogWarning.Active(GridLogWarning.isActive());
+    HadronsLogMessage.Active(GridLogMessage.isActive());
+    HadronsLogIterative.Active(GridLogIterative.isActive());
+    HadronsLogDebug.Active(GridLogDebug.isActive());
+    LOG(Message) << "Grid initialized" << std::endl;
+
+    // run setup ///////////////////////////////////////////////////////////////
+    Application              application;
+    std::vector<double>       mass    = {.01, .04, .2};
+    std::vector<std::string>  flavour = {"l", "s", "c"};
+    std::vector<std::string>  solvers = {"CG_l", "CG_s", "CG_c"};
+    std::string               kmom    = "0. 0. 0. 0.";
+    std::string               pmom    = "1. 0. 0. 0.";
+    std::string               qmom    = "-1. 0. 0. 0.";
+    std::string               mqmom   = "1. 0. 0. 0.";
+    std::vector<unsigned int> tKs     = {0};
+    unsigned int              dt_pi   = 16;
+    std::vector<unsigned int> tJs     = {8};
+    unsigned int              n_noise = 1;
+    unsigned int              nt      = 32;
+    bool                      do_disconnected(false);
+
+    // Global parameters.
+    Application::GlobalPar globalPar;
+    globalPar.trajCounter.start    = 1500;
+    globalPar.trajCounter.end      = 1520;
+    globalPar.trajCounter.step     = 20;
+    globalPar.seed                 = "1 2 3 4";
+    globalPar.genetic.maxGen       = 1000;
+    globalPar.genetic.maxCstGen    = 200;
+    globalPar.genetic.popSize      = 20;
+    globalPar.genetic.mutationRate = .1;
+    application.setPar(globalPar);
+
+    // gauge field
+    if (configStem == "None")
+    {
+        application.createModule<MGauge::Unit>("gauge");
+    }
+    else
+    {
+        MGauge::Load::Par gaugePar;
+        gaugePar.file = configStem;
+        application.createModule<MGauge::Load>("gauge", gaugePar);
+    }
+    for (unsigned int i = 0; i < flavour.size(); ++i)
+    {
+        // actions
+        MAction::DWF::Par actionPar;
+        actionPar.gauge = "gauge";
+        actionPar.Ls    = 16;
+        actionPar.M5    = 1.8;
+        actionPar.mass  = mass[i];
+        application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
+
+        // solvers
+        // RBPrecCG -> CG
+        MSolver::RBPrecCG::Par solverPar;
+        solverPar.action   = "DWF_" + flavour[i];
+        solverPar.residual = 1.0e-8;
+        application.createModule<MSolver::RBPrecCG>(solvers[i],
+                                                    solverPar);
+    }
+
+    // Create noise propagators for loops.
+    std::vector<std::string> noiseSrcs;
+    std::vector<std::vector<std::string>> noiseRes;
+    std::vector<std::vector<std::string>> noiseProps;
+    if (n_noise > 0)
+    {
+        MSource::Z2::Par noisePar;
+        noisePar.tA = 0;
+        noisePar.tB = nt - 1;
+        std::string loop_stem = "loop_";
+
+        noiseRes.resize(flavour.size());
+        noiseProps.resize(flavour.size());
+        for (unsigned int nn = 0; nn < n_noise; ++nn)
+        {
+            std::string eta = INIT_INDEX("noise", nn);
+            application.createModule<MSource::Z2>(eta, noisePar);
+            noiseSrcs.push_back(eta);
+
+            for (unsigned int f = 0; f < flavour.size(); ++f)
+            {
+                std::string loop_prop = INIT_INDEX(loop_stem + flavour[f], nn);
+                std::string loop_res  = loop_prop + "_res";
+                makePropagator(application, loop_res, eta, solvers[f]);
+                makeLoop(application, loop_prop, eta, loop_res);
+                noiseRes[f].push_back(loop_res);
+                noiseProps[f].push_back(loop_prop);
+            }
+        }
+    }
+
+    // Translate rare kaon decay across specified timeslices.
+    for (unsigned int i = 0; i < tKs.size(); ++i)
+    {
+        // Zero-momentum wall source propagators for kaon and pion.
+        unsigned int tK     = tKs[i];
+        unsigned int tpi    = (tK + dt_pi) % nt;
+        std::string q_Kl_0  = INIT_INDEX("Q_l_0", tK);
+        std::string q_pil_0 = INIT_INDEX("Q_l_0", tpi);
+        MAKE_WALL_PROP(tK, q_Kl_0, solvers[light]);
+        MAKE_WALL_PROP(tpi, q_pil_0, solvers[light]);
+
+        // Wall sources for kaon and pion with momentum insertion. If either
+        // p or k are zero, or p = k, re-use the existing name to avoid 
+        // duplicating a propagator.
+        std::string q_Ks_k  = INIT_INDEX("Q_Ks_k", tK);
+        std::string q_Ks_p  = INIT_INDEX((kmom == pmom) ? "Q_Ks_k" : "Q_Ks_p", tK);
+        std::string q_pil_k = INIT_INDEX((kmom == ZERO_MOM) ? "Q_l_0" : "Q_l_k", tpi);
+        std::string q_pil_p = INIT_INDEX((pmom == kmom) ? q_pil_k : ((pmom == ZERO_MOM) ? "Q_l_0" : "Q_l_p"), tpi);
+        MAKE_3MOM_WALL_PROP(tK, kmom, q_Ks_k, solvers[strange]);
+        MAKE_3MOM_WALL_PROP(tK, pmom, q_Ks_p, solvers[strange]);
+        MAKE_3MOM_WALL_PROP(tpi, kmom, q_pil_k, solvers[light]);
+        MAKE_3MOM_WALL_PROP(tpi, pmom, q_pil_p, solvers[light]);
+
+        /***********************************************************************
+         * CONTRACTIONS: pi and K 2pt contractions with mom = p, k.
+         **********************************************************************/
+        // Wall-Point
+        std::string PW_K_k = INIT_INDEX("PW_K_k", tK);
+        std::string PW_K_p = INIT_INDEX("PW_K_p", tK);
+        std::string PW_pi_k = INIT_INDEX("PW_pi_k", tpi);
+        std::string PW_pi_p = INIT_INDEX("PW_pi_p", tpi);
+        mesonContraction(application, 2, q_Kl_0, q_Ks_k, PW_K_k, kmom);
+        mesonContraction(application, 2, q_Kl_0, q_Ks_p, PW_K_p, pmom);
+        mesonContraction(application, 2, q_pil_k, q_pil_0, PW_pi_k, kmom);
+        mesonContraction(application, 2, q_pil_p, q_pil_0, PW_pi_p, pmom);
+        // Wall-Wall, to be done - requires modification of meson module.
+
+        /***********************************************************************
+         * CONTRACTIONS: 3pt Weak Hamiltonian, C & W (non-Eye type) classes.
+         **********************************************************************/
+        std::string HW_CW_k = LABEL_3PT("HW_CW_k", tK, tpi);
+        std::string HW_CW_p = LABEL_3PT("HW_CW_p", tK, tpi);
+        weakContractionNonEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, q_pil_0, HW_CW_k);
+        weakContractionNonEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, q_pil_0, HW_CW_p);
+
+        /***********************************************************************
+         * CONTRACTIONS: 3pt sd insertion.
+         **********************************************************************/
+        // Note: eventually will use wall sink smeared q_Kl_0 instead.
+        std::string sd_k = LABEL_3PT("sd_k", tK, tpi);
+        std::string sd_p = LABEL_3PT("sd_p", tK, tpi);
+        gamma3ptContraction(application, 3, q_Kl_0, q_Ks_k, q_pil_k, sd_k);
+        gamma3ptContraction(application, 3, q_Kl_0, q_Ks_p, q_pil_p, sd_p);
+
+        for (unsigned int nn = 0; nn < n_noise; ++nn)
+        {
+            /*******************************************************************
+             * CONTRACTIONS: 3pt Weak Hamiltonian, S and E (Eye type) classes.
+             ******************************************************************/
+            // Note: eventually will use wall sink smeared q_Kl_0 instead.
+            for (unsigned int f = 0; f < flavour.size(); ++f)
+            {
+                if ((f != strange) || do_disconnected)
+                {
+                    std::string HW_SE_k = LABEL_3PT("HW_SE_k_" + flavour[f], tK, tpi);
+                    std::string HW_SE_p = LABEL_3PT("HW_SE_p_" + flavour[f], tK, tpi);
+                    std::string loop_q  = noiseProps[f][nn];
+                    weakContractionEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, loop_q, HW_CW_k);
+                    weakContractionEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, loop_q, HW_CW_p);
+                }
+            }
+        }
+
+        // Perform separate contractions for each t_J position.
+        for (unsigned int j = 0; j < tJs.size(); ++j)
+        {
+            // Sequential sources for current insertions. Local for now,
+            // gamma_0 only.
+            unsigned int tJ = (tJs[j] + tK) % nt;
+            MSource::SeqGamma::Par seqPar;
+            std::string q_KlCl_q   = LABEL_3PT("Q_KlCl_q", tK, tJ);
+            std::string q_KsCs_mq  = LABEL_3PT("Q_KsCs_mq", tK, tJ);
+            std::string q_pilCl_q  = LABEL_3PT("Q_pilCl_q", tpi, tJ);
+            std::string q_pilCl_mq = LABEL_3PT("Q_pilCl_mq", tpi, tJ);
+            MAKE_SEQUENTIAL_PROP(tJ, q_Kl_0, qmom, q_KlCl_q, solvers[light]);
+            MAKE_SEQUENTIAL_PROP(tJ, q_Ks_k, mqmom, q_KsCs_mq, solvers[strange]);
+            MAKE_SEQUENTIAL_PROP(tJ, q_pil_p, qmom, q_pilCl_q, solvers[light]);
+            MAKE_SEQUENTIAL_PROP(tJ, q_pil_0, mqmom, q_pilCl_mq, solvers[light]);
+
+            /*******************************************************************
+             * CONTRACTIONS: pi and K 3pt contractions with current insertion.
+             ******************************************************************/
+            // Wall-Point
+            std::string C_PW_Kl   = LABEL_3PT("C_PW_Kl", tK, tJ);
+            std::string C_PW_Ksb  = LABEL_3PT("C_PW_Ksb", tK, tJ);
+            std::string C_PW_pilb = LABEL_3PT("C_PW_pilb", tK, tJ);
+            std::string C_PW_pil  = LABEL_3PT("C_PW_pil", tK, tJ);
+            mesonContraction(application, 3, q_KlCl_q, q_Ks_k, C_PW_Kl, pmom);
+            mesonContraction(application, 3, q_Kl_0, q_KsCs_mq, C_PW_Ksb, pmom);
+            mesonContraction(application, 3, q_pil_0, q_pilCl_q, C_PW_pilb, kmom);
+            mesonContraction(application, 3, q_pilCl_mq, q_pil_p, C_PW_pil, kmom);
+            // Wall-Wall, to be done.
+
+            /*******************************************************************
+             * CONTRACTIONS: 4pt contractions, C & W classes.
+             ******************************************************************/
+            std::string CW_Kl   = LABEL_4PT("CW_Kl", tK, tJ, tpi);
+            std::string CW_Ksb  = LABEL_4PT("CW_Ksb", tK, tJ, tpi);
+            std::string CW_pilb = LABEL_4PT("CW_pilb", tK, tJ, tpi);
+            std::string CW_pil  = LABEL_4PT("CW_pil", tK, tJ, tpi);
+            weakContractionNonEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, q_pil_0, CW_Kl);
+            weakContractionNonEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, q_pil_0, CW_Ksb);
+            weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, q_pil_0, CW_pilb);
+            weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, q_pilCl_mq, CW_pil);
+
+            /*******************************************************************
+             * CONTRACTIONS: 4pt contractions, sd insertions.
+             ******************************************************************/
+            // Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead.
+            std::string sd_Kl   = LABEL_4PT("sd_Kl", tK, tJ, tpi);
+            std::string sd_Ksb  = LABEL_4PT("sd_Ksb", tK, tJ, tpi);
+            std::string sd_pilb = LABEL_4PT("sd_pilb", tK, tJ, tpi);
+            gamma3ptContraction(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, sd_Kl);
+            gamma3ptContraction(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, sd_Ksb);
+            gamma3ptContraction(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, sd_pilb);
+
+            // Sequential sources for each noise propagator.
+            for (unsigned int nn = 0; nn < n_noise; ++nn)
+            {
+                std::string loop_stem = "loop_";
+
+                // Contraction required for each quark flavour - alternatively
+                // drop the strange loop if not performing disconnected
+                // contractions or neglecting H_W operators Q_3 -> Q_10.
+                for (unsigned int f = 0; f < flavour.size(); ++f)
+                {
+                    if ((f != strange) || do_disconnected)
+                    {
+                        std::string eta      = noiseSrcs[nn];
+                        std::string loop_q   = noiseProps[f][nn];
+                        std::string loop_qCq = LABEL_3PT(loop_stem + flavour[f], tJ, nn);
+                        std::string loop_qCq_res = loop_qCq + "_res";
+                        MAKE_SEQUENTIAL_PROP(tJ, noiseRes[f][nn], qmom, 
+                                             loop_qCq_res, solvers[f]);
+                        makeLoop(application, loop_qCq, eta, loop_qCq_res);
+
+                        /*******************************************************
+                         * CONTRACTIONS: 4pt contractions, S & E classes.
+                         ******************************************************/
+                        // Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead.
+                        std::string SE_Kl   = LABEL_4PT_NOISE("SE_Kl", tK, tJ, tpi, nn);
+                        std::string SE_Ksb  = LABEL_4PT_NOISE("SE_Ksb", tK, tJ, tpi, nn);
+                        std::string SE_pilb = LABEL_4PT_NOISE("SE_pilb", tK, tJ, tpi, nn);
+                        std::string SE_loop = LABEL_4PT_NOISE("SE_loop", tK, tJ, tpi, nn);
+                        weakContractionEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, loop_q, SE_Kl);
+                        weakContractionEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, loop_q, SE_Ksb);
+                        weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, SE_pilb);
+                        weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, loop_qCq, SE_loop);
+
+                        /*******************************************************
+                         * CONTRACTIONS: 4pt contractions, pi0 disconnected 
+                         * loop.
+                         ******************************************************/
+                        std::string disc0 = LABEL_4PT_NOISE("disc0", tK, tJ, tpi, nn);
+                        disc0Contraction(application, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, disc0);
+
+                        /*******************************************************
+                         * CONTRACTIONS: Disconnected loop.
+                         ******************************************************/
+                        std::string discLoop = "disc_" + loop_qCq;
+                        discLoopContraction(application, loop_qCq, discLoop);
+                    }
+                }
+            }
+        }
+    }
+    // execution
+    std::string par_file_name = "rarekaon_000_100_tK0_tpi16_tJ8_noloop_mc0.2.xml";
+    application.saveParameterFile(par_file_name);
+    application.run();
+
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons_spectrum.cc

@@ -96,12 +96,16 @@ int main(int argc, char *argv[])
         mesPar.output = "mesons/pt_" + flavour[i] + flavour[j];
         mesPar.q1     = "Qpt_" + flavour[i];
         mesPar.q2     = "Qpt_" + flavour[j];
+        mesPar.gammas = "all";
+        mesPar.mom    = "0. 0. 0. 0.";
         application.createModule<MContraction::Meson>("meson_pt_"
                                                       + flavour[i] + flavour[j],
                                                       mesPar);
         mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
         mesPar.q1     = "QZ2_" + flavour[i];
         mesPar.q2     = "QZ2_" + flavour[j];
+        mesPar.gammas = "all";
+        mesPar.mom    = "0. 0. 0. 0.";
         application.createModule<MContraction::Meson>("meson_Z2_"
                                                       + flavour[i] + flavour[j],
                                                       mesPar);

tests/qdpxx/Test_qdpxx_loops_staples.cc

@@ -317,6 +317,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
   Grid::QCD::LatticeComplex rect(UGrid);
   Grid::QCD::TComplex trect;
   Grid::QCD::Complex  crect;
+  Grid::RealD rrect;
   Grid::RealD vol = UGrid->gSites();
   for(int mu=0;mu<Grid::QCD::Nd;mu++){
   for(int nu=0;nu<Grid::QCD::Nd;nu++){
@@ -325,7 +326,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
       Grid::QCD::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu);
       trect = Grid::sum(rect);
       crect = Grid::TensorRemove(trect);
-      std::cout<< "mu/nu = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
+      rrect = real(crect);
+      std::cout<< "mu/nu = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/2.0/3.0<<std::endl;
 
       Grid::GridStopWatch Peter;
       Grid::GridStopWatch Azusa;
@@ -355,7 +357,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 
 	trect = Grid::sum(TrStap);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect=real(crect);
+	std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 
 	//              __ 
@@ -370,7 +373,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 
 	trect = Grid::sum(TrStap);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect=real(crect);
+	std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 	//           __ 
 	//          |__ __ |
@@ -384,7 +388,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 
 	trect = Grid::sum(TrStap);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 
 	//           __ ___ 
@@ -399,7 +404,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	TrStap = Grid::trace (U[mu]*Stap);
 	trect = Grid::sum(TrStap);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
   
 	//       --    
@@ -423,7 +429,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	SumTrStap += TrStap;
 	trect = Grid::sum(TrStap);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
   
 
 
@@ -441,11 +448,13 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	TrStap = Grid::trace (U[mu]*Stap);
 	trect = Grid::sum(TrStap);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 	trect = Grid::sum(SumTrStap);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline trace 2x1+1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline trace 2x1+1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/2.0/3.0<<std::endl;
       }
       Peter.Stop();
       Azusa.Start();
@@ -489,7 +498,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	RectPlaq_d = Grid::trace(U[mu]*ds_U);
 	trect = Grid::sum(RectPlaq_d);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 	//  __  __
 	// |__    |
@@ -501,7 +511,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	RectPlaq_d = Grid::trace(U[mu]*ds_U);
 	trect = Grid::sum(RectPlaq_d);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 	//      __
 	//  |__ __ |
@@ -513,7 +524,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	RectPlaq_d = Grid::trace(U[mu]*ds_U);
 	trect = Grid::sum(RectPlaq_d);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 
 	//   __
@@ -526,7 +538,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	RectPlaq_d = Grid::trace(U[mu]*ds_U);
 	trect = Grid::sum(RectPlaq_d);
 	crect = Grid::TensorRemove(trect);
-	std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 	   
 	// 1(mu) x 2 (nu)  ** this was ok
@@ -542,7 +555,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	trect = Grid::sum(RectPlaq_d);
 	crect = Grid::TensorRemove(trect);
 
-	std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	rrect = real(crect);
+	std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 
 	// 1(mu) x 2 (nu) ** this was ok
 	//   
@@ -570,8 +584,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
 	RectPlaq_d = Grid::trace(U[mu]*ds_U);
 	trect = Grid::sum(RectPlaq_d);
 	crect = Grid::TensorRemove(trect);
-
+	rrect = real(crect);
-	std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
+	std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
 	
       }
       Azusa.Stop();

tests/solver/Test_staggered_cg_unprec.cc

@@ -76,6 +76,7 @@ int main (int argc, char ** argv)
   ImprovedStaggeredFermionR Ds(Umu,Umu,Grid,RBGrid,mass);
 
   MdagMLinearOperator<ImprovedStaggeredFermionR,FermionField> HermOp(Ds);
+  ConjugateGradient<FermionField> CG(1.0e-6,10000);
   CG(HermOp,src,result);
 
   Grid_finalize();