diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 1ab53a56..05ad1697 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -35,6 +35,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc new file mode 100644 index 00000000..6685f292 --- /dev/null +++ b/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 + +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 + +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(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TWeakNeutral4ptDisc::getInput(void) +{ + std::vector in = {par().q1, par().q2, par().q3, par().q4}; + + return in; +} + +std::vector TWeakNeutral4ptDisc::getOutput(void) +{ + std::vector 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(par().q1); + PropagatorField &q2 = *env().template getObject(par().q2); + PropagatorField &q3 = *env().template getObject(par().q3); + PropagatorField &q4 = *env().template getObject(par().q4); + Gamma g5 = Gamma(Gamma::Algebra::Gamma5); + LatticeComplex expbuf(env().getGrid()); + std::vector corrbuf; + std::vector result(n_neut_disc_diag); + unsigned int ndim = env().getNd(); + + PropagatorField tmp(env().getGrid()); + std::vector meson(ndim, tmp); + std::vector 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); +} diff --git a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp new file mode 100644 index 00000000..1c014410 --- /dev/null +++ b/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 + +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 + +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 +{ +public: + TYPE_ALIASES(FIMPL,) + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, name, + std::vector, corr); + }; +public: + // constructor + TWeakNeutral4ptDisc(const std::string name); + // destructor + virtual ~TWeakNeutral4ptDisc(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector 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_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 911e20cd..af291631 100644 --- a/extras/Hadrons/modules.inc +++ b/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 \