Merge pull request #84 from Lanny91/feature/rare_kaon

Rare Kaon decay contraction code

extras/Hadrons/Modules.hpp

@@ -30,11 +30,16 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Grid/Hadrons/Modules/MAction/Wilson.hpp>
 #include <Grid/Hadrons/Modules/MContraction/Baryon.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/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/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,86 @@
+/*************************************************************************************
+
+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);
+};
+
+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,83 @@
+/*************************************************************************************
+
+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)
+
+class TWeakHamiltonianEye: 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
+    TWeakHamiltonianEye(const std::string name);
+    // destructor
+    virtual ~TWeakHamiltonianEye(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(WeakHamiltonianEye, TWeakHamiltonianEye, MContraction);
+
+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,82 @@
+/*************************************************************************************
+
+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)
+
+class TWeakHamiltonianNonEye: 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
+    TWeakHamiltonianNonEye(const std::string name);
+    // destructor
+    virtual ~TWeakHamiltonianNonEye(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(WeakHamiltonianNonEye, TWeakHamiltonianNonEye, MContraction);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WeakHamiltonianNonEye_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.inc

@@ -1,4 +1,6 @@
 modules_cc =\
+  Modules/MContraction/WeakHamiltonianEye.cc \
+  Modules/MContraction/WeakHamiltonianNonEye.cc \
   Modules/MGauge/Load.cc \
   Modules/MGauge/Random.cc \
   Modules/MGauge/Unit.cc
@@ -8,12 +10,17 @@ modules_hpp =\
   Modules/MAction/Wilson.hpp \
   Modules/MContraction/Baryon.hpp \
   Modules/MContraction/Meson.hpp \
+  Modules/MContraction/WeakHamiltonian.hpp \
+  Modules/MContraction/WeakHamiltonianEye.hpp \
+  Modules/MContraction/WeakHamiltonianNonEye.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
 

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;

tests/IO/Test_serialisation.cc

@@ -113,6 +113,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);
@@ -120,6 +121,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;
@@ -127,21 +130,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/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_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);