HADRONS: 4pt Weak + current disconnected topology (e.g. for rare neutral kaon decays)

extras/Hadrons/Modules.hpp

@@ -35,6 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #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>

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,84 @@
+/*************************************************************************************
+
+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))
+
+class TWeakNeutral4ptDisc: 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
+    TWeakNeutral4ptDisc(const std::string name);
+    // destructor
+    virtual ~TWeakNeutral4ptDisc(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(WeakNeutral4ptDisc, TWeakNeutral4ptDisc, MContraction);
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_WeakNeutral4ptDisc_hpp_

extras/Hadrons/modules.inc

@@ -1,6 +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
@@ -15,6 +16,7 @@ modules_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 \