From 4db82da0db456b85db49bdc2e499de248b7a37f2 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 12 Jan 2017 11:41:10 +0000 Subject: [PATCH 01/36] Wall sources --- extras/Hadrons/Modules/MSource/Wall.hpp | 91 +++++++++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 2 files changed, 92 insertions(+) create mode 100644 extras/Hadrons/Modules/MSource/Wall.hpp diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp new file mode 100644 index 00000000..36a58ada --- /dev/null +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -0,0 +1,91 @@ +#ifndef Hadrons_Wall_hpp_ +#define Hadrons_Wall_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Wall * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class WallPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WallPar, + unsigned int, tW); +}; + +template +class TWall: public Module +{ +public: + // constructor + TWall(const std::string name); + // destructor + virtual ~TWall(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(Wall, TWall, MSource); + +/****************************************************************************** + * TWall implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TWall::TWall(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TWall::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TWall::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TWall::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWall::execute(void) +{ + Lattice> t(env().getGrid()); + + LOG(Message) << "Generating wall source at t = " << par().tW << std::endl; + PropagatorField &src = *env().template createLattice(getName()); + LatticeCoordinate(t, Tp); + src = 1.; + src = where((t == par().tW), src, 0.*src); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Wall_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 4251ffa3..6270a88e 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -14,6 +14,7 @@ modules_hpp =\ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqGamma.hpp \ + Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ Modules/Quark.hpp From 3855673ebf4c5588da9efda615153b88294f7d0f Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 12 Jan 2017 11:42:37 +0000 Subject: [PATCH 02/36] Added header for wall source --- extras/Hadrons/Modules.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 77ae08b7..c017766e 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -36,5 +36,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include #include #include From f22b79da8fb9b7bcb35b25612e4630f5a53ec6ba Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 12 Jan 2017 12:52:12 +0000 Subject: [PATCH 03/36] Added missing type aliases --- extras/Hadrons/Modules/MSource/Wall.hpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp index 36a58ada..7a491fd2 100644 --- a/extras/Hadrons/Modules/MSource/Wall.hpp +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -22,6 +22,8 @@ public: template class TWall: public Module { +public: + TYPE_ALIASES(FImpl,); public: // constructor TWall(const std::string name); From b7f90aa011d17357df3bf3a82e2433f0f51c4f54 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 13 Jan 2017 15:54:19 +0000 Subject: [PATCH 04/36] Added momentum choice for wall source --- extras/Hadrons/Modules/MSource/Wall.hpp | 37 +++++++++++++++++++++---- 1 file changed, 32 insertions(+), 5 deletions(-) diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp index 7a491fd2..c299a7bf 100644 --- a/extras/Hadrons/Modules/MSource/Wall.hpp +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -7,6 +7,19 @@ BEGIN_HADRONS_NAMESPACE +/* + + Wall source + ----------------------------- + * src_x = theta(x_3 - tA) * theta(tB - x_3) * exp(i x.mom) + + * options: + - tA: begin timeslice (integer) + - tB: end timeslice (integer) + - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + + */ + /****************************************************************************** * Wall * ******************************************************************************/ @@ -16,7 +29,8 @@ class WallPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(WallPar, - unsigned int, tW); + unsigned int, tW, + std::string, mom); }; template @@ -76,14 +90,27 @@ void TWall::setup(void) // execution /////////////////////////////////////////////////////////////////// template void TWall::execute(void) -{ - Lattice> t(env().getGrid()); +{ + LOG(Message) << "Generating wall source at t = " << par().tW + << " with momentum " << par().mom << std::endl; - LOG(Message) << "Generating wall source at t = " << par().tW << std::endl; PropagatorField &src = *env().template createLattice(getName()); + Lattice> t(env().getGrid()); + LatticeComplex ph(env().getGrid()), coor(env().getGrid()); + std::vector p; + Complex i(0.0,1.0); + + p = strToVec(par().mom); + ph = zero; + for(unsigned int mu = 0; mu < Nd; mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + p[mu]*coor; + } + ph = exp(i*ph); LatticeCoordinate(t, Tp); src = 1.; - src = where((t == par().tW), src, 0.*src); + src = where((t == par().tW), src*ph, 0.*src); } END_MODULE_NAMESPACE From 78774fbdc08c93ceec42aebcc1beb3e3ba1ef222 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 17 Jan 2017 15:29:45 +0000 Subject: [PATCH 05/36] Construct loop propagator --- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MLoop/NoiseLoop.hpp | 91 ++++++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 93 insertions(+) create mode 100644 extras/Hadrons/Modules/MLoop/NoiseLoop.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index c017766e..bb5fba39 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -33,6 +33,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/MLoop/NoiseLoop.hpp b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp new file mode 100644 index 00000000..2354ff1d --- /dev/null +++ b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp @@ -0,0 +1,91 @@ +#ifndef Hadrons_NoiseLoop_hpp_ +#define Hadrons_NoiseLoop_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * NoiseLoop * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MLoop) + +class NoiseLoopPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(NoiseLoopPar, + std::string, q, + std::string, eta); +}; + +template +class TNoiseLoop: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TNoiseLoop(const std::string name); + // destructor + virtual ~TNoiseLoop(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(NoiseLoop, TNoiseLoop, MLoop); + +/****************************************************************************** + * TNoiseLoop implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TNoiseLoop::TNoiseLoop(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TNoiseLoop::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TNoiseLoop::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TNoiseLoop::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TNoiseLoop::execute(void) +{ + PropagatorField &loop = *env().template createLattice(getName()); + PropagatorField &q = *env().template getObject(par().q); + PropagatorField &eta = *env().template getObject(par().eta); + loop = q*adj(eta); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_NoiseLoop_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 6270a88e..ab46cbef 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -11,6 +11,7 @@ modules_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 \ From 82b7d4eaf03e977b4dc4261db026985c8a522c11 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 17 Jan 2017 15:58:32 +0000 Subject: [PATCH 06/36] Added noise loop dependencies --- extras/Hadrons/Modules/MLoop/NoiseLoop.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp index 2354ff1d..d1ede0c5 100644 --- a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp +++ b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp @@ -54,7 +54,7 @@ TNoiseLoop::TNoiseLoop(const std::string name) template std::vector TNoiseLoop::getInput(void) { - std::vector in; + std::vector in = {par().q, par().eta}; return in; } From 8ae1a95ec6545e0c29ac97908f2ce6f0bbbf43fd Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 17 Jan 2017 18:14:20 +0000 Subject: [PATCH 07/36] Legal banners and module descriptions --- extras/Hadrons/Modules/MLoop/NoiseLoop.hpp | 43 +++++++++++++++++++++- extras/Hadrons/Modules/MSource/Wall.hpp | 33 +++++++++++++++-- 2 files changed, 72 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp index d1ede0c5..3d2850d1 100644 --- a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp +++ b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MLoop/NoiseLoop.hpp + +Copyright (C) 2016 + +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_NoiseLoop_hpp_ #define Hadrons_NoiseLoop_hpp_ @@ -7,6 +35,19 @@ 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 * ******************************************************************************/ @@ -16,7 +57,7 @@ class NoiseLoopPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(NoiseLoopPar, - std::string, q, + std::string, q, std::string, eta); }; diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp index c299a7bf..9d22e23b 100644 --- a/extras/Hadrons/Modules/MSource/Wall.hpp +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -1,3 +1,31 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSource/Wall.hpp + +Copyright (C) 2016 + +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_Wall_hpp_ #define Hadrons_Wall_hpp_ @@ -11,11 +39,10 @@ BEGIN_HADRONS_NAMESPACE Wall source ----------------------------- - * src_x = theta(x_3 - tA) * theta(tB - x_3) * exp(i x.mom) + * src_x = delta(x_3 - tW) * exp(i x.mom) * options: - - tA: begin timeslice (integer) - - tB: end timeslice (integer) + - tW: source timeslice (integer) - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") */ From 07f2ebea1b8e3644d84af15e2be5649d0c34d58a Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 19 Jan 2017 22:18:42 +0000 Subject: [PATCH 08/36] Meson module now takes list of gamma matrices to insert at source and sink. --- extras/Hadrons/Modules/MContraction/Meson.hpp | 126 ++++++++++++++++-- tests/hadrons/Test_hadrons_meson_3pt.cc | 2 + tests/hadrons/Test_hadrons_spectrum.cc | 2 + 3 files changed, 116 insertions(+), 14 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index dc4f60ee..46c1f3a3 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -8,6 +8,7 @@ Copyright (C) 2015 Copyright (C) 2016 Author: Antonin Portelli + 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 @@ -36,6 +37,22 @@ 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 source & sink, pairs of gamma matrices + (space-separated integers) in square brackets, in a sequence + (e.g. "[15 7][7 15][7 7]"). + + Special values: "all" - perform all possible contractions. + + */ + /****************************************************************************** * TMeson * ******************************************************************************/ @@ -47,6 +64,7 @@ public: GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar, std::string, q1, std::string, q2, + std::string, gammas, std::string, output); }; @@ -70,6 +88,10 @@ public: // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); + virtual void parseGammaString(SpinMatrix*, + unsigned int &, + unsigned int &, + std::vector> &); // execution virtual void execute(void); }; @@ -102,6 +124,79 @@ std::vector TMeson::getOutput(void) return output; } +template +std::vector> strToVecPair(const std::string s) +{ + std::vector> v; + return v; +} + +template +void TMeson::parseGammaString(SpinMatrix *g, + unsigned int &n_snk, + unsigned int &n_src, + std::vector> &toDo) +{ + // Initialise counters for parsing gamma insertions at source & sink. + int empty = -1; + std::vector gamma_inds(Ns*Ns, empty); + unsigned int n_gam = 0; + + // Determine gamma matrices to insert at source/sink. + if (par().gammas.compare("all") == 0) + { + // Do all contractions. + toDo.resize(Ns*Ns); + for (int i = 0; i < Ns*Ns; ++i) + { + g[i] = makeGammaProd(i); + toDo[i].assign(Ns*Ns, true); + } + } + else + { + // Parse individual contractions from input string. + std::vector> gamma_pairs; + gamma_pairs = strToVecPair(par().gammas); + + // Function for gamma matrix counting & indexing at source/sink. + auto index_gamma = [&empty, &g, &gamma_inds, &n_gam](int i, + unsigned int &n) + { + if (i >= gamma_inds.size()) + { + HADRON_ERROR("Invalid gamma matrix index " << i); + } + if (gamma_inds[i] == empty) + { + g[n_gam] = makeGammaProd(i); + gamma_inds[i] = n_gam; + ++n_gam; + ++n; + } + }; + + // Count no. of unique gamma matrices, then construct matrix of + // contractions to do. + for (unsigned int i = 0; i < gamma_inds.size(); ++i) + { + index_gamma(gamma_pairs[i].first, n_snk); + index_gamma(gamma_pairs[i].second, n_src); + } + toDo.resize(n_gam); + for (int i = 0; i < n_gam; ++i) + { + toDo[i].assign(n_gam, false); + } + for (int i = 0; i < gamma_inds.size(); ++i) + { + toDo[gamma_inds[gamma_pairs[i].first]] + [gamma_inds[gamma_pairs[i].second]] = true; + } + } +} + + // execution /////////////////////////////////////////////////////////////////// template void TMeson::execute(void) @@ -115,26 +210,29 @@ void TMeson::execute(void) PropagatorField2 &q2 = *env().template getObject(par().q2); LatticeComplex c(env().getGrid()); SpinMatrix g[Ns*Ns], g5; + std::vector> toDo; std::vector buf; Result result; - + unsigned int n_snk, n_src; + g5 = makeGammaProd(Ns*Ns - 1); - result.corr.resize(Ns*Ns); - for (unsigned int i = 0; i < Ns*Ns; ++i) + parseGammaString(g, n_snk, n_src, toDo); + + result.corr.resize(n_snk); + for (unsigned int iSink = 0; iSink < toDo.size(); ++iSink) { - g[i] = makeGammaProd(i); - } - for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink) - { - result.corr[iSink].resize(Ns*Ns); - for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc) + result.corr[iSink].resize(n_src); + for (unsigned int iSrc = 0; iSrc < toDo.size(); ++iSrc) { - c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5); - sliceSum(c, buf, Tp); - result.corr[iSink][iSrc].resize(buf.size()); - for (unsigned int t = 0; t < buf.size(); ++t) + if (toDo[iSink][iSrc]) { - result.corr[iSink][iSrc][t] = TensorRemove(buf[t]); + c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5); + sliceSum(c, buf, Tp); + result.corr[iSink][iSrc].resize(buf.size()); + for (unsigned int t = 0; t < buf.size(); ++t) + { + result.corr[iSink][iSrc][t] = TensorRemove(buf[t]); + } } } } diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc index 56d0efa7..9166e6e2 100644 --- a/tests/hadrons/Test_hadrons_meson_3pt.cc +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -130,6 +130,7 @@ int main(int argc, char *argv[]) mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; mesPar.q1 = qName[i]; mesPar.q2 = qName[j]; + mesPar.gammas = "all"; application.createModule("meson_Z2_" + std::to_string(t) + "_" @@ -147,6 +148,7 @@ int main(int argc, char *argv[]) + std::to_string(mu); mesPar.q1 = qName[i]; mesPar.q2 = seqName[j][mu]; + mesPar.gammas = "all"; application.createModule("3pt_Z2_" + std::to_string(t) + "_" diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index ea6e4479..2c2d54c0 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -96,12 +96,14 @@ 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"; application.createModule("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"; application.createModule("meson_Z2_" + flavour[i] + flavour[j], mesPar); From af29be2c90c1b6ed08d8e1ba51bdc0afd33d3382 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 20 Jan 2017 16:38:50 +0000 Subject: [PATCH 09/36] Simplified operation of meson module. Result has been modified to output one contraction at a time for each pair of gamma insertions at source and sink. --- extras/Hadrons/Modules/MContraction/Meson.hpp | 154 +++++++----------- 1 file changed, 61 insertions(+), 93 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 46c1f3a3..a3e682ce 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -35,6 +35,22 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include +namespace Grid { + // Overload >> to extract gamma pair from "[g1 g2]" string. + template + inline std::istringstream &operator>>(std::istringstream &sstr, + std::pair &buf) + { + T1 buf1; + T2 buf2; + char c; + sstr >> c >> buf1 >> buf2 >> c; + sstr.peek(); + buf = std::make_pair(buf1, buf2); + return sstr; + } +} + BEGIN_HADRONS_NAMESPACE /* @@ -45,9 +61,9 @@ BEGIN_HADRONS_NAMESPACE * options: - q1: input propagator 1 (string) - q2: input propagator 2 (string) - - gammas: gamma products to insert at source & sink, pairs of gamma matrices - (space-separated integers) in square brackets, in a sequence - (e.g. "[15 7][7 15][7 7]"). + - gammas: gamma products to insert at sink & source, pairs of gamma matrices + (space-separated integers) in square brackets (i.e. [g_sink g_src]), + in a sequence (e.g. "[15 7][7 15][7 7]"). Special values: "all" - perform all possible contractions. @@ -58,6 +74,8 @@ BEGIN_HADRONS_NAMESPACE ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MContraction) +typedef std::pair GammaPair; + class MesonPar: Serializable { public: @@ -78,7 +96,9 @@ public: { public: GRID_SERIALIZABLE_CLASS_MEMBERS(Result, - std::vector>>, corr); + unsigned int, gamma_snk, + unsigned int, gamma_src, + std::vector, corr); }; public: // constructor @@ -88,10 +108,7 @@ public: // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); - virtual void parseGammaString(SpinMatrix*, - unsigned int &, - unsigned int &, - std::vector> &); + virtual void parseGammaString(std::vector &gammaList); // execution virtual void execute(void); }; @@ -124,75 +141,27 @@ std::vector TMeson::getOutput(void) return output; } -template -std::vector> strToVecPair(const std::string s) -{ - std::vector> v; - return v; -} - template -void TMeson::parseGammaString(SpinMatrix *g, - unsigned int &n_snk, - unsigned int &n_src, - std::vector> &toDo) +void TMeson::parseGammaString(std::vector &gammaList) { - // Initialise counters for parsing gamma insertions at source & sink. - int empty = -1; - std::vector gamma_inds(Ns*Ns, empty); - unsigned int n_gam = 0; - // Determine gamma matrices to insert at source/sink. if (par().gammas.compare("all") == 0) { // Do all contractions. - toDo.resize(Ns*Ns); - for (int i = 0; i < Ns*Ns; ++i) + unsigned int n_gam = Ns*Ns; + gammaList.resize(n_gam*n_gam); + for (unsigned int i = 0; i < n_gam; ++i) { - g[i] = makeGammaProd(i); - toDo[i].assign(Ns*Ns, true); + for (unsigned int j = 0; j < n_gam; ++j) + { + gammaList.push_back(std::make_pair(i, j)); + } } } else { // Parse individual contractions from input string. - std::vector> gamma_pairs; - gamma_pairs = strToVecPair(par().gammas); - - // Function for gamma matrix counting & indexing at source/sink. - auto index_gamma = [&empty, &g, &gamma_inds, &n_gam](int i, - unsigned int &n) - { - if (i >= gamma_inds.size()) - { - HADRON_ERROR("Invalid gamma matrix index " << i); - } - if (gamma_inds[i] == empty) - { - g[n_gam] = makeGammaProd(i); - gamma_inds[i] = n_gam; - ++n_gam; - ++n; - } - }; - - // Count no. of unique gamma matrices, then construct matrix of - // contractions to do. - for (unsigned int i = 0; i < gamma_inds.size(); ++i) - { - index_gamma(gamma_pairs[i].first, n_snk); - index_gamma(gamma_pairs[i].second, n_src); - } - toDo.resize(n_gam); - for (int i = 0; i < n_gam; ++i) - { - toDo[i].assign(n_gam, false); - } - for (int i = 0; i < gamma_inds.size(); ++i) - { - toDo[gamma_inds[gamma_pairs[i].first]] - [gamma_inds[gamma_pairs[i].second]] = true; - } + gammaList = strToVec(par().gammas); } } @@ -205,38 +174,37 @@ void TMeson::execute(void) << " quarks '" << par().q1 << "' and '" << par().q2 << "'" << std::endl; - XmlWriter writer(par().output); - PropagatorField1 &q1 = *env().template getObject(par().q1); - PropagatorField2 &q2 = *env().template getObject(par().q2); - LatticeComplex c(env().getGrid()); - SpinMatrix g[Ns*Ns], g5; - std::vector> toDo; - std::vector buf; - Result result; - unsigned int n_snk, n_src; + XmlWriter writer(par().output); + PropagatorField1 &q1 = *env().template getObject(par().q1); + PropagatorField2 &q2 = *env().template getObject(par().q2); + LatticeComplex c(env().getGrid()); + SpinMatrix g[Ns*Ns], g5; + std::vector gammaList; + std::vector buf; + std::vector result; g5 = makeGammaProd(Ns*Ns - 1); - parseGammaString(g, n_snk, n_src, toDo); - - result.corr.resize(n_snk); - for (unsigned int iSink = 0; iSink < toDo.size(); ++iSink) + for (int i = 0; i < Ns*Ns; ++i) { - result.corr[iSink].resize(n_src); - for (unsigned int iSrc = 0; iSrc < toDo.size(); ++iSrc) - { - if (toDo[iSink][iSrc]) - { - c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5); - sliceSum(c, buf, Tp); - result.corr[iSink][iSrc].resize(buf.size()); - for (unsigned int t = 0; t < buf.size(); ++t) - { - result.corr[iSink][iSrc][t] = TensorRemove(buf[t]); - } - } - } + g[i] = makeGammaProd(i); + } + parseGammaString(gammaList); + + result.resize(gammaList.size()); + for (unsigned int i = 0; i < result.size(); ++i) + { + c = trace(g[gammaList[i].first]*q1*g[gammaList[i].second]*g5*adj(q2)*g5); + 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[i]); } - write(writer, "meson", result); } END_MODULE_NAMESPACE From 7dd2764bb26da9bf1c14929c058d71f5df0fe467 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 23 Jan 2017 15:17:54 +0000 Subject: [PATCH 10/36] Wall sink smearing --- extras/Hadrons/Modules.hpp | 1 + extras/Hadrons/Modules/MSink/Wall.hpp | 144 ++++++++++++++++++++++++++ extras/Hadrons/modules.inc | 1 + 3 files changed, 146 insertions(+) create mode 100644 extras/Hadrons/Modules/MSink/Wall.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index bb5fba39..e379bd77 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -34,6 +34,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/MSink/Wall.hpp b/extras/Hadrons/Modules/MSink/Wall.hpp new file mode 100644 index 00000000..324afa3a --- /dev/null +++ b/extras/Hadrons/Modules/MSink/Wall.hpp @@ -0,0 +1,144 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSink/Wall.hpp + +Copyright (C) 2016 + +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_Wall_hpp_ +#define Hadrons_Wall_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Wall sink smearing + ----------------------------- + * prop_x_3 = sum_x_i (prop_x * exp(-i x.mom)) + + * options: + - q: input propagator + - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + + */ + +/****************************************************************************** + * Wall * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSink) + +class WallPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WallPar, + std::string, q, + std::string, mom); +}; + +template +class TWall: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TWall(const std::string name); + // destructor + virtual ~TWall(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(Wall, TWall, MSink); + +/****************************************************************************** + * TWall implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TWall::TWall(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TWall::getInput(void) +{ + std::vector in = {par().q}; + + return in; +} + +template +std::vector TWall::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TWall::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWall::execute(void) +{ + LOG(Message) << "Wall smearing " << par().q << std::endl; + + PropagatorField &q = *env().template getObject(par().q); + std::vector prop; + LatticeComplex ph(env().getGrid()), coor(env().getGrid()); + std::vector p; + Complex i(0.0,1.0); + + p = strToVec(par().mom); + ph = zero; + for(unsigned int mu = 0; mu < Nd; mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + p[mu]*coor; + } + ph = exp(-i*ph); + sliceSum(ph*q, prop, Tp); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Wall_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index ab46cbef..b2494236 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -12,6 +12,7 @@ modules_hpp =\ Modules/MGauge/Random.hpp \ Modules/MGauge/Unit.hpp \ Modules/MLoop/NoiseLoop.hpp \ + Modules/MSink/Wall.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqGamma.hpp \ From 977f34dca64fafd6119d8475a3d3ab8f8e23aebe Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 26 Jan 2017 13:18:33 +0000 Subject: [PATCH 11/36] Added missing typename --- extras/Hadrons/Modules/MSink/Wall.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/MSink/Wall.hpp b/extras/Hadrons/Modules/MSink/Wall.hpp index 324afa3a..64966218 100644 --- a/extras/Hadrons/Modules/MSink/Wall.hpp +++ b/extras/Hadrons/Modules/MSink/Wall.hpp @@ -121,7 +121,7 @@ void TWall::execute(void) LOG(Message) << "Wall smearing " << par().q << std::endl; PropagatorField &q = *env().template getObject(par().q); - std::vector prop; + std::vector prop; LatticeComplex ph(env().getGrid()), coor(env().getGrid()); std::vector p; Complex i(0.0,1.0); From 9bf4108d1f9c54c4b644d9da5a5aa7214db7dd0a Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Fri, 27 Jan 2017 16:58:11 +0000 Subject: [PATCH 12/36] Weak Hamiltonian contraction modules, for Eye and NonEye contraction topologies. Execution for NonEye type diagrams has been implemented, but not yet for Eye type. --- extras/Hadrons/Modules.hpp | 3 + .../Modules/MContraction/WeakHamiltonian.hpp | 86 ++++++++++ .../MContraction/WeakHamiltonianEye.cc | 92 +++++++++++ .../MContraction/WeakHamiltonianEye.hpp | 75 +++++++++ .../MContraction/WeakHamiltonianNonEye.cc | 148 ++++++++++++++++++ .../MContraction/WeakHamiltonianNonEye.hpp | 81 ++++++++++ extras/Hadrons/modules.inc | 5 + 7 files changed, 490 insertions(+) create mode 100644 extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp create mode 100644 extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc create mode 100644 extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp create mode 100644 extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc create mode 100644 extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index e379bd77..46c65cb2 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -30,6 +30,9 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include +#include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp new file mode 100644 index 00000000..cb3775ab --- /dev/null +++ b/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 + +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 +#include +#include + +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 < Nd; ++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_ diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc new file mode 100644 index 00000000..86920a7c --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -0,0 +1,92 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.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-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 | q4 + * /-<-¬ | /-<-¬ + * / \ | / \ + * \ / | \ / + * q2 \ / q3 | \_ _/ + * |-----<-----* *-----<----¬ | q2 * q3 + * i * H_W * f | |-----<------*------<----¬ + * \ / | i * H_W * 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(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TWeakHamiltonianEye::getInput(void) +{ + std::vector in; + + return in; +} + +std::vector TWeakHamiltonianEye::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TWeakHamiltonianEye::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +void TWeakHamiltonianEye::execute(void) +{ + +} diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp new file mode 100644 index 00000000..8f88a78b --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -0,0 +1,75 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.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_WeakHamiltonianEye_hpp_ +#define Hadrons_WeakHamiltonianEye_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * WeakHamiltonianEye * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +enum +{ + S_diag = 0, + E_diag = 1 +}; + +// 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 +{ +public: + TYPE_ALIASES(FIMPL,) +public: + // constructor + TWeakHamiltonianEye(const std::string name); + // destructor + virtual ~TWeakHamiltonianEye(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(WeakHamiltonianEye, TWeakHamiltonianEye, MContraction); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WeakHamiltonianEye_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc new file mode 100644 index 00000000..81d18f15 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -0,0 +1,148 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.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-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(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TWeakHamiltonianNonEye::getInput(void) +{ + std::vector in; + + return in; +} + +std::vector TWeakHamiltonianNonEye::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TWeakHamiltonianNonEye::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +void TWeakHamiltonianNonEye::execute(void) +{ + XmlWriter 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); + SpinMatrix g5 = makeGammaProd(Ns*Ns - 1); // Soon to be deprecated. + LatticeComplex expbuf(env().getGrid()); + std::vector corrbuf; + std::vector result; + + PropagatorField tmp1(env().getGrid()); + LatticeComplex tmp2(env().getGrid()); + std::vector C_i_side_loop(Nd, tmp1); + std::vector C_f_side_loop(Nd, tmp1); + std::vector W_i_side_loop(Nd, tmp2); + std::vector W_f_side_loop(Nd, tmp2); + + // Soon to be deprecated. Keeping V and A parts distinct to take advantage + // of zero-flop gamma products, when implemented. + std::vector> gL; + gL.resize(n_i); + gL[i_V].resize(Nd); + gL[i_A].resize(Nd); + for (unsigned int mu = 0; mu < Nd; ++mu) + { + gL[i_V][mu] = makeGammaProd(mu); + gL[i_A][mu] = g5*makeGammaProd(mu); + } + + // Setup for C-type contractions. + for (int mu = 0; mu < Nd; ++mu) + { + C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, gL[i_V][mu]) - + MAKE_CW_SUBDIAG(q1, q2, gL[i_A][mu]); + C_f_side_loop[mu] = MAKE_CW_SUBDIAG(q3, q4, gL[i_V][mu]) - + MAKE_CW_SUBDIAG(q3, q4, gL[i_A][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], "3pt_HW_C") + + // Recycle sub-expressions for W-type contractions. + for (unsigned int mu = 0; mu < Nd; ++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], "3pt_HW_W") + + write(writer, "HW_3pt_NonEye", result[C_diag]); + write(writer, "HW_3pt_NonEye", result[W_diag]); +} diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp new file mode 100644 index 00000000..8541df9e --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp @@ -0,0 +1,81 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.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_WeakHamiltonianNonEye_hpp_ +#define Hadrons_WeakHamiltonianNonEye_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * WeakHamiltonianNonEye * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +enum +{ + W_diag = 0, + C_diag = 1 +}; + +// 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 +{ +public: + TYPE_ALIASES(FIMPL,) + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, name, + std::vector, corr); + }; +public: + // constructor + TWeakHamiltonianNonEye(const std::string name); + // destructor + virtual ~TWeakHamiltonianNonEye(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(WeakHamiltonianNonEye, TWeakHamiltonianNonEye, MContraction); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WeakHamiltonianNonEye_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index b2494236..69794808 100644 --- a/extras/Hadrons/modules.inc +++ b/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,6 +10,9 @@ 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 \ From c4d367272001476f4c905a5c3dffd5efb880f273 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 30 Jan 2017 17:09:04 +0000 Subject: [PATCH 13/36] Hadrons: Momentum projection in meson module. --- extras/Hadrons/Modules/MContraction/Meson.hpp | 21 ++++++++++++++++--- tests/hadrons/Test_hadrons_meson_3pt.cc | 2 ++ tests/hadrons/Test_hadrons_spectrum.cc | 2 ++ 3 files changed, 22 insertions(+), 3 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index a3e682ce..3a8c2533 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MContraction/Meson.hpp Copyright (C) 2015 Copyright (C) 2016 +Copyright (C) 2017 Author: Antonin Portelli Andrew Lawson @@ -66,8 +67,9 @@ BEGIN_HADRONS_NAMESPACE in a sequence (e.g. "[15 7][7 15][7 7]"). 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 * @@ -83,6 +85,7 @@ public: std::string, q1, std::string, q2, std::string, gammas, + std::string, mom, std::string, output); }; @@ -182,7 +185,19 @@ void TMeson::execute(void) std::vector gammaList; std::vector buf; std::vector result; + std::vector p; + p = strToVec(par().mom); + LatticeComplex ph(env().getGrid()), coor(env().getGrid()); + Complex i(0.0,1.0); + ph = zero; + for(unsigned int mu = 0; mu < env().getNd(); mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); + } + ph = exp(-2*M_PI*i*ph); + g5 = makeGammaProd(Ns*Ns - 1); for (int i = 0; i < Ns*Ns; ++i) { @@ -193,7 +208,7 @@ void TMeson::execute(void) result.resize(gammaList.size()); for (unsigned int i = 0; i < result.size(); ++i) { - c = trace(g[gammaList[i].first]*q1*g[gammaList[i].second]*g5*adj(q2)*g5); + c = trace(g[gammaList[i].first]*q1*g[gammaList[i].second]*g5*adj(q2)*g5*ph); sliceSum(c, buf, Tp); result[i].gamma_snk = gammaList[i].first; diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc index 9166e6e2..efef6931 100644 --- a/tests/hadrons/Test_hadrons_meson_3pt.cc +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -131,6 +131,7 @@ int main(int argc, char *argv[]) mesPar.q1 = qName[i]; mesPar.q2 = qName[j]; mesPar.gammas = "all"; + mesPar.mom = "0. 0. 0. 0."; application.createModule("meson_Z2_" + std::to_string(t) + "_" @@ -149,6 +150,7 @@ int main(int argc, char *argv[]) mesPar.q1 = qName[i]; mesPar.q2 = seqName[j][mu]; mesPar.gammas = "all"; + mesPar.mom = "0. 0. 0. 0."; application.createModule("3pt_Z2_" + std::to_string(t) + "_" diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index 2c2d54c0..2d731ff4 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -97,6 +97,7 @@ int main(int argc, char *argv[]) mesPar.q1 = "Qpt_" + flavour[i]; mesPar.q2 = "Qpt_" + flavour[j]; mesPar.gammas = "all"; + mesPar.mom = "0. 0. 0. 0."; application.createModule("meson_pt_" + flavour[i] + flavour[j], mesPar); @@ -104,6 +105,7 @@ int main(int argc, char *argv[]) mesPar.q1 = "QZ2_" + flavour[i]; mesPar.q2 = "QZ2_" + flavour[j]; mesPar.gammas = "all"; + mesPar.mom = "0. 0. 0. 0."; application.createModule("meson_Z2_" + flavour[i] + flavour[j], mesPar); From 651e1a7cbc2d6dd865d37b13d768f0532f9d01f5 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 30 Jan 2017 17:14:33 +0000 Subject: [PATCH 14/36] Hadrons: Momentum inserted as multiples of 2*pi/L --- extras/Hadrons/Modules/MSink/Wall.hpp | 6 +++--- extras/Hadrons/Modules/MSource/SeqGamma.hpp | 5 +++-- extras/Hadrons/Modules/MSource/Wall.hpp | 6 +++--- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/extras/Hadrons/Modules/MSink/Wall.hpp b/extras/Hadrons/Modules/MSink/Wall.hpp index 64966218..761357a9 100644 --- a/extras/Hadrons/Modules/MSink/Wall.hpp +++ b/extras/Hadrons/Modules/MSink/Wall.hpp @@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSink/Wall.hpp -Copyright (C) 2016 +Copyright (C) 2017 Author: Andrew Lawson @@ -131,9 +131,9 @@ void TWall::execute(void) for(unsigned int mu = 0; mu < Nd; 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(-2*M_PI*i*ph); sliceSum(ph*q, prop, Tp); } diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp index 611b0108..814da273 100644 --- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp +++ b/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 @@ -150,9 +151,9 @@ void TSeqGamma::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(2*M_PI*i*ph); LatticeCoordinate(t, Tp); src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q); } diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp index 9d22e23b..e94829bc 100644 --- a/extras/Hadrons/Modules/MSource/Wall.hpp +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -4,7 +4,7 @@ Grid physics library, www.github.com/paboyle/Grid Source file: extras/Hadrons/Modules/MSource/Wall.hpp -Copyright (C) 2016 +Copyright (C) 2017 Author: Andrew Lawson @@ -132,9 +132,9 @@ void TWall::execute(void) for(unsigned int mu = 0; mu < Nd; 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(2*M_PI*i*ph); LatticeCoordinate(t, Tp); src = 1.; src = where((t == par().tW), src*ph, 0.*src); From 9e9f621d5dcc3e0f784b07548576cb3ed19f12e5 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 30 Jan 2017 17:54:21 +0000 Subject: [PATCH 15/36] Hadrons: added Weak Hamiltonian module dependencies, some reformatting. --- .../MContraction/WeakHamiltonianEye.cc | 20 +++++++++--------- .../MContraction/WeakHamiltonianEye.hpp | 7 +++++++ .../MContraction/WeakHamiltonianNonEye.cc | 21 ++++++++++--------- 3 files changed, 28 insertions(+), 20 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index 86920a7c..3da2c5a1 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -38,16 +38,16 @@ using namespace MContraction; * 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 | q4 - * /-<-¬ | /-<-¬ - * / \ | / \ - * \ / | \ / - * q2 \ / q3 | \_ _/ - * |-----<-----* *-----<----¬ | q2 * q3 - * i * H_W * f | |-----<------*------<----¬ - * \ / | i * H_W * f + * Schematics: q4 | + * /-<-¬ | + * / \ | q2 q3 + * \ / | /----<------*------<----¬ + * q2 \ / q3 | / /-*-¬ \ + * /-----<-----* *-----<----¬ | / / \ \ + * i * H_W * f | i * \ / q4 * f + * \ / | \ \->-/ / * \ / | \ / - * \---------->---------/ | \------------>---------/ + * \---------->---------/ | \----------->----------/ * q1 | q1 * | * Saucer (S) | Eye (E) @@ -67,7 +67,7 @@ TWeakHamiltonianEye::TWeakHamiltonianEye(const std::string name) // dependencies/products /////////////////////////////////////////////////////// std::vector TWeakHamiltonianEye::getInput(void) { - std::vector in; + std::vector in = {par().q1, par().q2, par().q3, par().q4}; return in; } diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp index 8f88a78b..2b4c8c31 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -52,6 +52,13 @@ class TWeakHamiltonianEye: public Module { public: TYPE_ALIASES(FIMPL,) + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, name, + std::vector, corr); + }; public: // constructor TWeakHamiltonianEye(const std::string name); diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index 81d18f15..31763dac 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -69,7 +69,7 @@ TWeakHamiltonianNonEye::TWeakHamiltonianNonEye(const std::string name) // dependencies/products /////////////////////////////////////////////////////// std::vector TWeakHamiltonianNonEye::getInput(void) { - std::vector in; + std::vector in = {par().q1, par().q2, par().q3, par().q4}; return in; } @@ -99,28 +99,29 @@ void TWeakHamiltonianNonEye::execute(void) LatticeComplex expbuf(env().getGrid()); std::vector corrbuf; std::vector result; + unsigned int ndim = env().getNd(); PropagatorField tmp1(env().getGrid()); LatticeComplex tmp2(env().getGrid()); - std::vector C_i_side_loop(Nd, tmp1); - std::vector C_f_side_loop(Nd, tmp1); - std::vector W_i_side_loop(Nd, tmp2); - std::vector W_f_side_loop(Nd, tmp2); + std::vector C_i_side_loop(ndim, tmp1); + std::vector C_f_side_loop(ndim, tmp1); + std::vector W_i_side_loop(ndim, tmp2); + std::vector W_f_side_loop(ndim, tmp2); // Soon to be deprecated. Keeping V and A parts distinct to take advantage // of zero-flop gamma products, when implemented. std::vector> gL; gL.resize(n_i); - gL[i_V].resize(Nd); - gL[i_A].resize(Nd); - for (unsigned int mu = 0; mu < Nd; ++mu) + gL[i_V].resize(ndim); + gL[i_A].resize(ndim); + for (unsigned int mu = 0; mu < ndim; ++mu) { gL[i_V][mu] = makeGammaProd(mu); gL[i_A][mu] = g5*makeGammaProd(mu); } // Setup for C-type contractions. - for (int mu = 0; mu < Nd; ++mu) + for (int mu = 0; mu < ndim; ++mu) { C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, gL[i_V][mu]) - MAKE_CW_SUBDIAG(q1, q2, gL[i_A][mu]); @@ -133,7 +134,7 @@ void TWeakHamiltonianNonEye::execute(void) MAKE_DIAG(expbuf, corrbuf, result[C_diag], "3pt_HW_C") // Recycle sub-expressions for W-type contractions. - for (unsigned int mu = 0; mu < Nd; ++mu) + 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]); From fa237401ffb96bfe44fc0a72751142399f7e4171 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 1 Feb 2017 12:56:55 +0000 Subject: [PATCH 16/36] Consistent variable name in macro --- extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp index cb3775ab..f9b0481f 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp @@ -56,7 +56,7 @@ for (unsigned int t = 0; t < buf.size(); ++t)\ //// Contraction of mu index: use 'mu' variable in exp. #define SUM_MU(buf,exp)\ buf = zero;\ -for (unsigned int mu = 0; mu < Nd; ++mu)\ +for (unsigned int mu = 0; mu < ndim; ++mu)\ {\ buf += exp;\ } From 35ac85aea89f567360c7df61335a1dd4207fed83 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 1 Feb 2017 12:57:34 +0000 Subject: [PATCH 17/36] Updated Weak Hamiltonian contractions to use zero-flop gamma matrices --- .../MContraction/WeakHamiltonianNonEye.cc | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index 31763dac..106c5964 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -95,7 +95,7 @@ void TWeakHamiltonianNonEye::execute(void) PropagatorField &q2 = *env().template getObject(par().q2); PropagatorField &q3 = *env().template getObject(par().q3); PropagatorField &q4 = *env().template getObject(par().q4); - SpinMatrix g5 = makeGammaProd(Ns*Ns - 1); // Soon to be deprecated. + Gamma g5 = Gamma(Gamma::Algebra::Gamma5); LatticeComplex expbuf(env().getGrid()); std::vector corrbuf; std::vector result; @@ -108,17 +108,14 @@ void TWeakHamiltonianNonEye::execute(void) std::vector W_i_side_loop(ndim, tmp2); std::vector W_f_side_loop(ndim, tmp2); - // Soon to be deprecated. Keeping V and A parts distinct to take advantage - // of zero-flop gamma products, when implemented. - std::vector> gL; - gL.resize(n_i); - gL[i_V].resize(ndim); - gL[i_A].resize(ndim); - for (unsigned int mu = 0; mu < ndim; ++mu) - { - gL[i_V][mu] = makeGammaProd(mu); - gL[i_A][mu] = g5*makeGammaProd(mu); - } + std::vector> gL = {{Gamma(Gamma::Algebra::GammaX), + Gamma(Gamma::Algebra::GammaY), + Gamma(Gamma::Algebra::GammaZ), + Gamma(Gamma::Algebra::GammaT)}, + {Gamma(Gamma::Algebra::GammaXGamma5), + Gamma(Gamma::Algebra::GammaYGamma5), + Gamma(Gamma::Algebra::GammaZGamma5), + Gamma(Gamma::Algebra::GammaTGamma5)}}; // Setup for C-type contractions. for (int mu = 0; mu < ndim; ++mu) From 5be05d85b8bd591ad9d2385ae9ce950b139add3a Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 1 Feb 2017 13:56:22 +0000 Subject: [PATCH 18/36] Fixed collision of Wall source and sink header ifdefs --- extras/Hadrons/Modules/MSink/Wall.hpp | 6 +++--- extras/Hadrons/Modules/MSource/Wall.hpp | 6 +++--- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/extras/Hadrons/Modules/MSink/Wall.hpp b/extras/Hadrons/Modules/MSink/Wall.hpp index 761357a9..81780f62 100644 --- a/extras/Hadrons/Modules/MSink/Wall.hpp +++ b/extras/Hadrons/Modules/MSink/Wall.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Wall_hpp_ -#define Hadrons_Wall_hpp_ +#ifndef Hadrons_WallSink_hpp_ +#define Hadrons_WallSink_hpp_ #include #include @@ -141,4 +141,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Wall_hpp_ +#endif // Hadrons_WallSink_hpp_ diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp index e94829bc..a8e88ee9 100644 --- a/extras/Hadrons/Modules/MSource/Wall.hpp +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo *************************************************************************************/ /* END LEGAL */ -#ifndef Hadrons_Wall_hpp_ -#define Hadrons_Wall_hpp_ +#ifndef Hadrons_WallSource_hpp_ +#define Hadrons_WallSource_hpp_ #include #include @@ -144,4 +144,4 @@ END_MODULE_NAMESPACE END_HADRONS_NAMESPACE -#endif // Hadrons_Wall_hpp_ +#endif // Hadrons_WallSource_hpp_ From 74a5cda84b2ad0b54c1441b4305a46f7ab53989e Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 1 Feb 2017 15:03:49 +0000 Subject: [PATCH 19/36] Removed unnecessary "3pt" labels --- .../Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index 106c5964..fed2e701 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -128,7 +128,7 @@ void TWeakHamiltonianNonEye::execute(void) // 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], "3pt_HW_C") + MAKE_DIAG(expbuf, corrbuf, result[C_diag], "HW_C") // Recycle sub-expressions for W-type contractions. for (unsigned int mu = 0; mu < ndim; ++mu) @@ -139,8 +139,8 @@ void TWeakHamiltonianNonEye::execute(void) // Perform W-type contractions. SUM_MU(expbuf, W_i_side_loop[mu]*W_f_side_loop[mu]) - MAKE_DIAG(expbuf, corrbuf, result[W_diag], "3pt_HW_W") + MAKE_DIAG(expbuf, corrbuf, result[W_diag], "HW_W") - write(writer, "HW_3pt_NonEye", result[C_diag]); - write(writer, "HW_3pt_NonEye", result[W_diag]); + write(writer, "HW_NonEye", result[C_diag]); + write(writer, "HW_NonEye", result[W_diag]); } From d35d87d2c208a7d05463923771a9247e0e0e6ebd Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 1 Feb 2017 16:33:24 +0000 Subject: [PATCH 20/36] Weak Hamiltonian Eye-type contraction execution --- .../MContraction/WeakHamiltonianEye.cc | 52 +++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index 3da2c5a1..758f2164 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -88,5 +88,57 @@ void TWeakHamiltonianEye::setup(void) // execution /////////////////////////////////////////////////////////////////// void TWeakHamiltonianEye::execute(void) { + XmlWriter 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; + unsigned int ndim = env().getNd(); + PropagatorField tmp1(env().getGrid()); + LatticeComplex tmp2(env().getGrid()); + std::vector S_body(ndim, tmp1); + std::vector S_loop(ndim, tmp1); + std::vector E_body(ndim, tmp2); + std::vector E_loop(ndim, tmp2); + + std::vector> gL = {{Gamma(Gamma::Algebra::GammaX), + Gamma(Gamma::Algebra::GammaY), + Gamma(Gamma::Algebra::GammaZ), + Gamma(Gamma::Algebra::GammaT)}, + {Gamma(Gamma::Algebra::GammaXGamma5), + Gamma(Gamma::Algebra::GammaYGamma5), + Gamma(Gamma::Algebra::GammaZGamma5), + Gamma(Gamma::Algebra::GammaTGamma5)}}; + + // Setup for S-type contractions. + for (int mu = 0; mu < ndim; ++mu) + { + S_body[mu] = MAKE_SE_BODY(q1, q2, q3, gL[i_V][mu]) - + MAKE_SE_BODY(q1, q2, q3, gL[i_A][mu]); + S_loop[mu] = MAKE_SE_LOOP(q4, gL[i_V][mu]) - + MAKE_SE_LOOP(q4, gL[i_A][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[S_diag]); + write(writer, "HW_Eye", result[E_diag]); } From 522f6bf91ae85ef33aea2dc51ec77b545944242a Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 1 Feb 2017 16:36:08 +0000 Subject: [PATCH 21/36] Add function to say whether or not a module exists in application class --- extras/Hadrons/Application.cc | 5 +++++ extras/Hadrons/Application.hpp | 1 + 2 files changed, 6 insertions(+) diff --git a/extras/Hadrons/Application.cc b/extras/Hadrons/Application.cc index 62674f30..c1d9d722 100644 --- a/extras/Hadrons/Application.cc +++ b/extras/Hadrons/Application.cc @@ -91,6 +91,11 @@ const Application::GlobalPar & Application::getPar(void) return par_; } +bool Application::hasModule(const std::string name) +{ + return env().hasModule(name); +} + // execute ///////////////////////////////////////////////////////////////////// void Application::run(void) { diff --git a/extras/Hadrons/Application.hpp b/extras/Hadrons/Application.hpp index fce9b6eb..2e3d7997 100644 --- a/extras/Hadrons/Application.hpp +++ b/extras/Hadrons/Application.hpp @@ -86,6 +86,7 @@ public: void createModule(const std::string name); template void createModule(const std::string name, const typename M::Par &par); + bool hasModule(const std::string name); // execute void run(void); // XML parameter file I/O From 1e257a12512d2eda08411337ddf341f9796ccc8a Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 1 Feb 2017 16:36:40 +0000 Subject: [PATCH 22/36] Hadrons: test for rare kaon contraction code. --- tests/hadrons/Test_hadrons_rarekaon.cc | 461 +++++++++++++++++++++++++ 1 file changed, 461 insertions(+) create mode 100644 tests/hadrons/Test_hadrons_rarekaon.cc diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc new file mode 100644 index 00000000..d48e2718 --- /dev/null +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -0,0 +1,461 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_rarekaon.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. + *******************************************************************************/ + +#include + +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) + +// 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 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);\ +} + +/******************************************************************************* + * Functions for propagator construction. + ******************************************************************************/ + +/******************************************************************************* + * 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 (!(application.hasModule(srcName))) + { + MSource::SeqGamma::Par seqPar; + seqPar.q = qSrc; + seqPar.tA = tS; + seqPar.tB = tS; + seqPar.mom = mom; + application.createModule(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 (!(application.hasModule(srcName))) + { + MSource::Wall::Par wallPar; + wallPar.tW = tW; + wallPar.mom = mom; + application.createModule(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. + if (!(application.hasModule(wallName))) + { + MSink::Wall::Par wallPar; + wallPar.q = propName; + wallPar.mom = mom; + application.createModule(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 (!(application.hasModule(propName))) + { + Quark::Par quarkPar; + quarkPar.source = srcName; + quarkPar.solver = solver; + application.createModule(propName, quarkPar); + } +} + +/******************************************************************************* + * 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) + * 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 modName = std::to_string(npt) + "pt_" + label; + if (!(application.hasModule(modName))) + { + MContraction::Meson::Par mesPar; + mesPar.output = std::to_string(npt) + "pt/" + label; + mesPar.q1 = q1; + mesPar.q2 = q2; + mesPar.gammas = "[1 1]"; + application.createModule(modName, mesPar); + } + } + +/******************************************************************************* + * 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 (!(application.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(modName, weakPar);\ + }\ +} +HW_CONTRACTION(Eye) // weakContractionEye +HW_CONTRACTION(NonEye) // weakContractionNonEye + +int main(int argc, char *argv[]) +{ + // parse command line ////////////////////////////////////////////////////// + std::string configStem; + + if (argc < 2) + { + std::cerr << "usage: " << argv[0] << " [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; + unsigned int light = 0; + unsigned int strange = 1; + unsigned int charm = 2; + std::vector mass = {.01, .04, .2}; + std::vector flavour = {"l", "s", "c"}; + std::vector 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 tKs = {0, 16}; + unsigned int dt_pi = 16; + std::vector tJs = {8}; + unsigned int n_noise = 0; + unsigned int nt = 32; + + // 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 + MGauge::Load::Par gaugePar; + gaugePar.file = configStem; + application.createModule("gauge", gaugePar); + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 12; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + application.createModule("DWF_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule(solvers[i], + solverPar); + } + + // Create noise propagators for loops. + if (n_noise > 0) + { + std::vector noiseProps_l(n_noise); + std::vector noiseProps_c(n_noise); + MSource::Z2::Par noisePar; + noisePar.tA = 0; + noisePar.tB = nt - 1; + for (unsigned int nn = 0; nn < n_noise; ++nn) + { + std::string eta = INIT_INDEX("noise", nn); + std::string loop_l = INIT_INDEX("loop_l", nn); + std::string loop_c = INIT_INDEX("loop_c", nn); + + application.createModule(eta, noisePar); + makePropagator(application, loop_l, eta, solvers[light]); + makePropagator(application, loop_c, eta, solvers[charm]); + noiseProps_l.push_back(loop_l); + noiseProps_c.push_back(loop_c); + } + } + + // 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]); + + // Wall sink-smeared propagators. + std::string q_Kl_0_W = q_Kl_0 + "_W"; + std::string q_Ks_k_W = q_Ks_k + "_W"; + std::string q_Ks_p_W = q_Ks_p + "_W"; + std::string q_pil_0_W = q_pil_0 + "_W"; + std::string q_pil_k_W = q_pil_k + "_W"; + std::string q_pil_p_W = q_pil_p + "_W"; + makeWallSink(application, q_Kl_0, q_Kl_0_W); + makeWallSink(application, q_Ks_k, q_Ks_k_W, kmom); + makeWallSink(application, q_Ks_p, q_Ks_p_W, pmom); + makeWallSink(application, q_pil_0, q_pil_0_W); + makeWallSink(application, q_pil_k, q_pil_k_W, kmom); + makeWallSink(application, q_pil_p, q_pil_p_W, pmom); + + /*********************************************************************** + * 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. + **********************************************************************/ + + for (unsigned int nn = 0; nn < n_noise; ++nn) + { + /******************************************************************* + * CONTRACTIONS: 3pt Weak Hamiltonian, S and E (Eye type) classes. + ******************************************************************/ + + } + + // Perform separate contractions for each t_J position. + for (unsigned int i = 0; i < tJs.size(); ++i) + { + // Sequential sources for current insertions. Local for now, + // gamma_0 only. + unsigned int tJ = (tJs[i] + 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. + ******************************************************************/ + + // Sequential sources for each noise propagator. + for (unsigned int nn = 0; nn < n_noise; ++nn) + { + /*************************************************************** + * CONTRACTIONS: 4pt contractions, S & E classes. + **************************************************************/ + + /*************************************************************** + * CONTRACTIONS: 4pt contractions, pi0 disconnected loop. + **************************************************************/ + } + } + } + // execution + application.saveParameterFile("rarekaon_000_100_tK0_tpi28_tJ14_noloop_mc0.2.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} From b7cd1a19e3ada7f9e3006402c0965a4690397e17 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 6 Feb 2017 14:08:59 +0000 Subject: [PATCH 23/36] Utilities for reading and writing "pair" objects. --- extras/Hadrons/Modules/MContraction/Meson.hpp | 20 +------- lib/serialisation/BaseIO.h | 48 +++++++++++++++++++ tests/IO/Test_serialisation.cc | 12 +++++ 3 files changed, 62 insertions(+), 18 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 270875a6..a46dd445 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -36,22 +36,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include -namespace Grid { - // Overload >> to extract gamma pair from "" string. - template - inline std::istringstream &operator>>(std::istringstream &sstr, - std::pair &buf) - { - unsigned int buf1; - unsigned int buf2; - char c; - sstr >> c >> buf1 >> buf2 >> c; - sstr.peek(); - buf = std::make_pair((T1)buf1, (T2)buf2); - return sstr; - } -} - BEGIN_HADRONS_NAMESPACE /* @@ -63,8 +47,8 @@ BEGIN_HADRONS_NAMESPACE - q1: input propagator 1 (string) - q2: input propagator 2 (string) - gammas: gamma products to insert at sink & source, pairs of gamma matrices - (space-separated integers) in square brackets (i.e. [g_sink g_src]), - in a sequence (e.g. "[15 7][7 15][7 7]"). + (space-separated strings) in angled brackets (i.e. ), + in a sequence (e.g. ""). Special values: "all" - perform all possible contractions. - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0."), diff --git a/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 0357915d..36e6fd77 100644 --- a/lib/serialisation/BaseIO.h +++ b/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 "" + template + inline std::istream & operator>>(std::istream &is, std::pair &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 + inline std::ostream & operator<<(std::ostream &os, const std::pair &p) + { + os << "<" << p.first << " " << p.second << ">"; + return os; + } + // Abstract writer/reader classes //////////////////////////////////////////// // static polymorphism implemented using CRTP idiom class Serializable; diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 8204b05b..acafc900 100644 --- a/tests/IO/Test_serialisation.cc +++ b/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 vec; + std::pair 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("iotest.xml", obj, "XML (object) "); ioTest("iotest.xml", vec, "XML (vector of objects)"); + ioTest("iotest.xml", pair, "XML (pair of objects)"); //// binary ioTest("iotest.bin", obj, "binary (object) "); ioTest("iotest.bin", vec, "binary (vector of objects)"); + ioTest("iotest.bin", pair, "binary (pair of objects)"); //// text ioTest("iotest.dat", obj, "text (object) "); ioTest("iotest.dat", vec, "text (vector of objects)"); + ioTest("iotest.dat", pair, "text (pair of objects)"); //// HDF5 +#undef HAVE_HDF5 #ifdef HAVE_HDF5 ioTest("iotest.h5", obj, "HDF5 (object) "); ioTest("iotest.h5", vec, "HDF5 (vector of objects)"); + ioTest("iotest.h5", pair, "HDF5 (pair of objects)"); #endif std::cout << "\n==== vector flattening/reconstruction" << std::endl; From 1a122a0dd819d2131c1e2ccda82608730ad29df1 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 6 Feb 2017 17:35:41 +0000 Subject: [PATCH 24/36] Hadrons: corrected gamma matrix inputs in rare kaon test --- tests/hadrons/Test_hadrons_rarekaon.cc | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc index d48e2718..3976d034 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -88,6 +88,7 @@ inline void makeSequentialSource(Application &application, std::string srcName, seqPar.tA = tS; seqPar.tB = tS; seqPar.mom = mom; + seqPar.gamma = Gamma::Algebra::GammaT; application.createModule(srcName, seqPar); } } @@ -184,7 +185,7 @@ inline void mesonContraction(Application &application, unsigned int npt, mesPar.output = std::to_string(npt) + "pt/" + label; mesPar.q1 = q1; mesPar.q2 = q2; - mesPar.gammas = "[1 1]"; + mesPar.gammas = ""; application.createModule(modName, mesPar); } } From d6a7d7d1e0005b6f6eae311f095e39b612464bc6 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 6 Feb 2017 18:15:49 +0000 Subject: [PATCH 25/36] Hadrons: added missing momentum parameter in rare kaon contraction test --- tests/hadrons/Test_hadrons_rarekaon.cc | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc index 3976d034..0601eeb6 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -185,6 +185,7 @@ inline void mesonContraction(Application &application, unsigned int npt, mesPar.output = std::to_string(npt) + "pt/" + label; mesPar.q1 = q1; mesPar.q2 = q2; + mesPar.mom = mom; mesPar.gammas = ""; application.createModule(modName, mesPar); } From 4a45c06dd77eb8402cf1ee6b37702e966a6ac14c Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Mon, 6 Feb 2017 20:12:30 +0000 Subject: [PATCH 26/36] Code cleaning and addition of Weak Hamiltonian contraction log message --- .../MContraction/WeakHamiltonianEye.cc | 22 ++++++++----------- .../MContraction/WeakHamiltonianEye.hpp | 2 +- .../MContraction/WeakHamiltonianNonEye.cc | 22 ++++++++----------- .../MContraction/WeakHamiltonianNonEye.hpp | 2 +- 4 files changed, 20 insertions(+), 28 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index 758f2164..0df6ec2f 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -88,6 +88,11 @@ 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; + XmlWriter writer(par().output); PropagatorField &q1 = *env().template getObject(par().q1); PropagatorField &q2 = *env().template getObject(par().q2); @@ -106,22 +111,13 @@ void TWeakHamiltonianEye::execute(void) std::vector E_body(ndim, tmp2); std::vector E_loop(ndim, tmp2); - std::vector> gL = {{Gamma(Gamma::Algebra::GammaX), - Gamma(Gamma::Algebra::GammaY), - Gamma(Gamma::Algebra::GammaZ), - Gamma(Gamma::Algebra::GammaT)}, - {Gamma(Gamma::Algebra::GammaXGamma5), - Gamma(Gamma::Algebra::GammaYGamma5), - Gamma(Gamma::Algebra::GammaZGamma5), - Gamma(Gamma::Algebra::GammaTGamma5)}}; - // Setup for S-type contractions. for (int mu = 0; mu < ndim; ++mu) { - S_body[mu] = MAKE_SE_BODY(q1, q2, q3, gL[i_V][mu]) - - MAKE_SE_BODY(q1, q2, q3, gL[i_A][mu]); - S_loop[mu] = MAKE_SE_LOOP(q4, gL[i_V][mu]) - - MAKE_SE_LOOP(q4, gL[i_A][mu]); + S_body[mu] = MAKE_SE_BODY(q1, q2, q3, Gamma::gmu[mu]) - + MAKE_SE_BODY(q1, q2, q3, Gamma::gmu[mu]*g5); + S_loop[mu] = MAKE_SE_LOOP(q4, Gamma::gmu[mu]) - + MAKE_SE_LOOP(q4, Gamma::gmu[mu]*g5); } // Perform S-type contractions. diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp index 2b4c8c31..6f5a95c0 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -45,7 +45,7 @@ enum }; // 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_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 diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index fed2e701..148b9afa 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -90,6 +90,11 @@ 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; + XmlWriter writer(par().output); PropagatorField &q1 = *env().template getObject(par().q1); PropagatorField &q2 = *env().template getObject(par().q2); @@ -108,22 +113,13 @@ void TWeakHamiltonianNonEye::execute(void) std::vector W_i_side_loop(ndim, tmp2); std::vector W_f_side_loop(ndim, tmp2); - std::vector> gL = {{Gamma(Gamma::Algebra::GammaX), - Gamma(Gamma::Algebra::GammaY), - Gamma(Gamma::Algebra::GammaZ), - Gamma(Gamma::Algebra::GammaT)}, - {Gamma(Gamma::Algebra::GammaXGamma5), - Gamma(Gamma::Algebra::GammaYGamma5), - Gamma(Gamma::Algebra::GammaZGamma5), - Gamma(Gamma::Algebra::GammaTGamma5)}}; - // Setup for C-type contractions. for (int mu = 0; mu < ndim; ++mu) { - C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, gL[i_V][mu]) - - MAKE_CW_SUBDIAG(q1, q2, gL[i_A][mu]); - C_f_side_loop[mu] = MAKE_CW_SUBDIAG(q3, q4, gL[i_V][mu]) - - MAKE_CW_SUBDIAG(q3, q4, gL[i_A][mu]); + C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, Gamma::gmu[mu]) - + MAKE_CW_SUBDIAG(q1, q2, Gamma::gmu[mu]*g5); + C_f_side_loop[mu] = MAKE_CW_SUBDIAG(q3, q4, Gamma::gmu[mu]) - + MAKE_CW_SUBDIAG(q3, q4, Gamma::gmu[mu]*g5); } // Perform C-type contractions. diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp index 8541df9e..f0d26b54 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp @@ -45,7 +45,7 @@ enum }; // Wing and Connected subdiagram contractions -#define MAKE_CW_SUBDIAG(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma) +#define MAKE_CW_SUBDIAG(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*(g5*gamma)) class TWeakHamiltonianNonEye: public Module { From bdd2765461d839f6d167d368b1a2c4cc311bf45f Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 7 Feb 2017 13:06:42 +0000 Subject: [PATCH 27/36] Added missing allocation of Weak Hamiltonian result vector --- extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc | 2 +- extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp | 3 ++- extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc | 2 +- extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp | 3 ++- 4 files changed, 6 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index 0df6ec2f..6c13aee5 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -101,7 +101,7 @@ void TWeakHamiltonianEye::execute(void) Gamma g5 = Gamma(Gamma::Algebra::Gamma5); LatticeComplex expbuf(env().getGrid()); std::vector corrbuf; - std::vector result; + std::vector result(n_eye_diag); unsigned int ndim = env().getNd(); PropagatorField tmp1(env().getGrid()); diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp index 6f5a95c0..a07d814f 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -41,7 +41,8 @@ BEGIN_MODULE_NAMESPACE(MContraction) enum { S_diag = 0, - E_diag = 1 + E_diag = 1, + n_eye_diag = 2 }; // Saucer and Eye subdiagram contractions. diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index 148b9afa..98a5e104 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -103,7 +103,7 @@ void TWeakHamiltonianNonEye::execute(void) Gamma g5 = Gamma(Gamma::Algebra::Gamma5); LatticeComplex expbuf(env().getGrid()); std::vector corrbuf; - std::vector result; + std::vector result(n_noneye_diag); unsigned int ndim = env().getNd(); PropagatorField tmp1(env().getGrid()); diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp index f0d26b54..a594cfda 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp @@ -41,7 +41,8 @@ BEGIN_MODULE_NAMESPACE(MContraction) enum { W_diag = 0, - C_diag = 1 + C_diag = 1, + n_noneye_diag = 2 }; // Wing and Connected subdiagram contractions From 2b2fc6453f34210b99350eaec7de85e02fdd0784 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Tue, 7 Feb 2017 13:59:29 +0000 Subject: [PATCH 28/36] Fixed single precision compatibility issues --- extras/Hadrons/Modules/MContraction/Meson.hpp | 2 +- extras/Hadrons/Modules/MSink/Wall.hpp | 2 +- extras/Hadrons/Modules/MSource/SeqGamma.hpp | 2 +- extras/Hadrons/Modules/MSource/Wall.hpp | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index a46dd445..1ea46eb6 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -181,7 +181,7 @@ void TMeson::execute(void) LatticeCoordinate(coor, mu); ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); } - ph = exp(-2*M_PI*i*ph); + ph = exp(-(Real)(2*M_PI)*i*ph); parseGammaString(gammaList); diff --git a/extras/Hadrons/Modules/MSink/Wall.hpp b/extras/Hadrons/Modules/MSink/Wall.hpp index 81780f62..9c39528b 100644 --- a/extras/Hadrons/Modules/MSink/Wall.hpp +++ b/extras/Hadrons/Modules/MSink/Wall.hpp @@ -133,7 +133,7 @@ void TWall::execute(void) LatticeCoordinate(coor, mu); ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); } - ph = exp(-2*M_PI*i*ph); + ph = exp(-(Real)(2*M_PI)*i*ph); sliceSum(ph*q, prop, Tp); } diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp index fab3751f..366ebee7 100644 --- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp +++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp @@ -152,7 +152,7 @@ void TSeqGamma::execute(void) LatticeCoordinate(coor, mu); ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); } - ph = exp(2*M_PI*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); } diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp index a8e88ee9..8722876f 100644 --- a/extras/Hadrons/Modules/MSource/Wall.hpp +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -134,7 +134,7 @@ void TWall::execute(void) LatticeCoordinate(coor, mu); ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); } - ph = exp(2*M_PI*i*ph); + ph = exp((Real)(2*M_PI)*i*ph); LatticeCoordinate(t, Tp); src = 1.; src = where((t == par().tW), src*ph, 0.*src); From 06f7ee202e8413fca6bd1da97579700acbd05576 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 8 Feb 2017 10:08:18 +0000 Subject: [PATCH 29/36] Revert "Add function to say whether or not a module exists in application class" This reverts commit 522f6bf91ae85ef33aea2dc51ec77b545944242a. --- extras/Hadrons/Application.cc | 5 ----- extras/Hadrons/Application.hpp | 1 - 2 files changed, 6 deletions(-) diff --git a/extras/Hadrons/Application.cc b/extras/Hadrons/Application.cc index c1d9d722..62674f30 100644 --- a/extras/Hadrons/Application.cc +++ b/extras/Hadrons/Application.cc @@ -91,11 +91,6 @@ const Application::GlobalPar & Application::getPar(void) return par_; } -bool Application::hasModule(const std::string name) -{ - return env().hasModule(name); -} - // execute ///////////////////////////////////////////////////////////////////// void Application::run(void) { diff --git a/extras/Hadrons/Application.hpp b/extras/Hadrons/Application.hpp index 2e3d7997..fce9b6eb 100644 --- a/extras/Hadrons/Application.hpp +++ b/extras/Hadrons/Application.hpp @@ -86,7 +86,6 @@ public: void createModule(const std::string name); template void createModule(const std::string name, const typename M::Par &par); - bool hasModule(const std::string name); // execute void run(void); // XML parameter file I/O From b9f7ea47c3b2d4eace16677b052e5eb36d1ae693 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 8 Feb 2017 10:10:06 +0000 Subject: [PATCH 30/36] Access hasModule function directly from Environment instance. --- tests/hadrons/Test_hadrons_rarekaon.cc | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc index 0601eeb6..a6feacad 100644 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -81,7 +81,7 @@ inline void makeSequentialSource(Application &application, std::string srcName, std::string mom = ZERO_MOM) { // If the source already exists, don't make the module again. - if (!(application.hasModule(srcName))) + if (!(Environment::getInstance().hasModule(srcName))) { MSource::SeqGamma::Par seqPar; seqPar.q = qSrc; @@ -106,7 +106,7 @@ 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 (!(application.hasModule(srcName))) + if (!(Environment::getInstance().hasModule(srcName))) { MSource::Wall::Par wallPar; wallPar.tW = tW; @@ -128,7 +128,7 @@ 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. - if (!(application.hasModule(wallName))) + if (!(Environment::getInstance().hasModule(wallName))) { MSink::Wall::Par wallPar; wallPar.q = propName; @@ -150,7 +150,7 @@ 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 (!(application.hasModule(propName))) + if (!(Environment::getInstance().hasModule(propName))) { Quark::Par quarkPar; quarkPar.source = srcName; @@ -179,7 +179,7 @@ inline void mesonContraction(Application &application, unsigned int npt, std::string &label, std::string mom = ZERO_MOM) { std::string modName = std::to_string(npt) + "pt_" + label; - if (!(application.hasModule(modName))) + if (!(Environment::getInstance().hasModule(modName))) { MContraction::Meson::Par mesPar; mesPar.output = std::to_string(npt) + "pt/" + label; @@ -211,7 +211,7 @@ inline void weakContraction##top(Application &application, unsigned int npt,\ std::string &label)\ {\ std::string modName = std::to_string(npt) + "pt_" + label;\ - if (!(application.hasModule(modName)))\ + if (!(Environment::getInstance().hasModule(modName)))\ {\ MContraction::WeakHamiltonian##top::Par weakPar;\ weakPar.output = std::to_string(npt) + "pt/" + label;\ @@ -452,7 +452,7 @@ int main(int argc, char *argv[]) } } // execution - application.saveParameterFile("rarekaon_000_100_tK0_tpi28_tJ14_noloop_mc0.2.xml"); + application.saveParameterFile("rarekaon_000_100_tK0_tpi16_tJ8_noloop_mc0.2.xml"); application.run(); // epilogue From e5a7ed43628516e675b1a470702dc625a4edfad4 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 8 Feb 2017 12:29:33 +0000 Subject: [PATCH 31/36] Moved write outside of loop, some physics corrections --- extras/Hadrons/Modules/MContraction/Meson.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 1ea46eb6..9581b8d5 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -181,7 +181,7 @@ void TMeson::execute(void) LatticeCoordinate(coor, mu); ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); } - ph = exp(-(Real)(2*M_PI)*i*ph); + ph = exp((Real)(2*M_PI)*i*ph); parseGammaString(gammaList); @@ -190,7 +190,7 @@ void TMeson::execute(void) { Gamma gSnk(gammaList[i].first); Gamma gSrc(gammaList[i].second); - c = trace((g5*gSnk)*q1*(gSrc*g5)*adj(q2))*ph; + c = trace((g5*gSnk)*q1*(adj(gSrc)*g5)*adj(q2))*ph; sliceSum(c, buf, Tp); result[i].gamma_snk = gammaList[i].first; @@ -200,8 +200,8 @@ void TMeson::execute(void) { result[i].corr[t] = TensorRemove(buf[t]); } - write(writer, "meson", result[i]); } + write(writer, "meson", result); } END_MODULE_NAMESPACE From 6ebf8b12b6bc180902c00e4d5235a6f4c8f6aa67 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 8 Feb 2017 12:43:13 +0000 Subject: [PATCH 32/36] Removed unnecessary repeat of write in Weak Hamiltonian contractions --- extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc | 3 +-- extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index 6c13aee5..19808bfa 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -135,6 +135,5 @@ void TWeakHamiltonianEye::execute(void) SUM_MU(expbuf, E_body[mu]*E_loop[mu]) MAKE_DIAG(expbuf, corrbuf, result[E_diag], "HW_E") - write(writer, "HW_Eye", result[S_diag]); - write(writer, "HW_Eye", result[E_diag]); + write(writer, "HW_Eye", result); } diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index 98a5e104..b94311bc 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -137,6 +137,5 @@ void TWeakHamiltonianNonEye::execute(void) 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[C_diag]); - write(writer, "HW_NonEye", result[W_diag]); + write(writer, "HW_NonEye", result); } From beba824136e7ed753d8d10fde987fcd138895f8c Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 8 Feb 2017 12:45:39 +0000 Subject: [PATCH 33/36] Make use of GammaL class in Weak Hamiltonian contractions --- extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc | 6 ++---- extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp | 2 +- .../Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc | 6 ++---- .../Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp | 2 +- 4 files changed, 6 insertions(+), 10 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index 19808bfa..4a1c9e39 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -114,10 +114,8 @@ void TWeakHamiltonianEye::execute(void) // Setup for S-type contractions. for (int mu = 0; mu < ndim; ++mu) { - S_body[mu] = MAKE_SE_BODY(q1, q2, q3, Gamma::gmu[mu]) - - MAKE_SE_BODY(q1, q2, q3, Gamma::gmu[mu]*g5); - S_loop[mu] = MAKE_SE_LOOP(q4, Gamma::gmu[mu]) - - MAKE_SE_LOOP(q4, Gamma::gmu[mu]*g5); + 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. diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp index a07d814f..2b855ce0 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -46,7 +46,7 @@ enum }; // 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_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 diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index b94311bc..7945f9cb 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -116,10 +116,8 @@ void TWeakHamiltonianNonEye::execute(void) // Setup for C-type contractions. for (int mu = 0; mu < ndim; ++mu) { - C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, Gamma::gmu[mu]) - - MAKE_CW_SUBDIAG(q1, q2, Gamma::gmu[mu]*g5); - C_f_side_loop[mu] = MAKE_CW_SUBDIAG(q3, q4, Gamma::gmu[mu]) - - MAKE_CW_SUBDIAG(q3, q4, Gamma::gmu[mu]*g5); + 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. diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp index a594cfda..b7fbdf28 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp @@ -46,7 +46,7 @@ enum }; // Wing and Connected subdiagram contractions -#define MAKE_CW_SUBDIAG(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*(g5*gamma)) +#define MAKE_CW_SUBDIAG(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma) class TWeakHamiltonianNonEye: public Module { From fc19503673d580d5a41dce916838360349ca88d5 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 8 Feb 2017 13:17:39 +0000 Subject: [PATCH 34/36] Removed MSink namespace. --- extras/Hadrons/Modules.hpp | 1 - extras/Hadrons/Modules/MSink/Wall.hpp | 144 -------------------------- extras/Hadrons/modules.inc | 1 - 3 files changed, 146 deletions(-) delete mode 100644 extras/Hadrons/Modules/MSink/Wall.hpp diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 46c65cb2..e8cdf2af 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -37,7 +37,6 @@ 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/MSink/Wall.hpp b/extras/Hadrons/Modules/MSink/Wall.hpp deleted file mode 100644 index 9c39528b..00000000 --- a/extras/Hadrons/Modules/MSink/Wall.hpp +++ /dev/null @@ -1,144 +0,0 @@ -/************************************************************************************* - -Grid physics library, www.github.com/paboyle/Grid - -Source file: extras/Hadrons/Modules/MSink/Wall.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_WallSink_hpp_ -#define Hadrons_WallSink_hpp_ - -#include -#include -#include - -BEGIN_HADRONS_NAMESPACE - -/* - - Wall sink smearing - ----------------------------- - * prop_x_3 = sum_x_i (prop_x * exp(-i x.mom)) - - * options: - - q: input propagator - - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") - - */ - -/****************************************************************************** - * Wall * - ******************************************************************************/ -BEGIN_MODULE_NAMESPACE(MSink) - -class WallPar: Serializable -{ -public: - GRID_SERIALIZABLE_CLASS_MEMBERS(WallPar, - std::string, q, - std::string, mom); -}; - -template -class TWall: public Module -{ -public: - TYPE_ALIASES(FImpl,); -public: - // constructor - TWall(const std::string name); - // destructor - virtual ~TWall(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(Wall, TWall, MSink); - -/****************************************************************************** - * TWall implementation * - ******************************************************************************/ -// constructor ///////////////////////////////////////////////////////////////// -template -TWall::TWall(const std::string name) -: Module(name) -{} - -// dependencies/products /////////////////////////////////////////////////////// -template -std::vector TWall::getInput(void) -{ - std::vector in = {par().q}; - - return in; -} - -template -std::vector TWall::getOutput(void) -{ - std::vector out = {getName()}; - - return out; -} - -// setup /////////////////////////////////////////////////////////////////////// -template -void TWall::setup(void) -{ - -} - -// execution /////////////////////////////////////////////////////////////////// -template -void TWall::execute(void) -{ - LOG(Message) << "Wall smearing " << par().q << std::endl; - - PropagatorField &q = *env().template getObject(par().q); - std::vector prop; - LatticeComplex ph(env().getGrid()), coor(env().getGrid()); - std::vector p; - Complex i(0.0,1.0); - - p = strToVec(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); - sliceSum(ph*q, prop, Tp); -} - -END_MODULE_NAMESPACE - -END_HADRONS_NAMESPACE - -#endif // Hadrons_WallSink_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 69794808..3e2f1020 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -17,7 +17,6 @@ modules_hpp =\ Modules/MGauge/Random.hpp \ Modules/MGauge/Unit.hpp \ Modules/MLoop/NoiseLoop.hpp \ - Modules/MSink/Wall.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqGamma.hpp \ From d7a1dc85beadda24ac7475ee34424b3ba9e4b5ae Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 8 Feb 2017 13:23:05 +0000 Subject: [PATCH 35/36] Revert "Hadrons: test for rare kaon contraction code." This reverts commit 1e257a12512d2eda08411337ddf341f9796ccc8a. --- tests/hadrons/Test_hadrons_rarekaon.cc | 463 ------------------------- 1 file changed, 463 deletions(-) delete mode 100644 tests/hadrons/Test_hadrons_rarekaon.cc diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc deleted file mode 100644 index a6feacad..00000000 --- a/tests/hadrons/Test_hadrons_rarekaon.cc +++ /dev/null @@ -1,463 +0,0 @@ -/******************************************************************************* - Grid physics library, www.github.com/paboyle/Grid - - Source file: tests/hadrons/Test_hadrons_rarekaon.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. - *******************************************************************************/ - -#include - -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) - -// 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 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);\ -} - -/******************************************************************************* - * Functions for propagator construction. - ******************************************************************************/ - -/******************************************************************************* - * 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(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(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. - if (!(Environment::getInstance().hasModule(wallName))) - { - MSink::Wall::Par wallPar; - wallPar.q = propName; - wallPar.mom = mom; - application.createModule(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(propName, quarkPar); - } -} - -/******************************************************************************* - * 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) - * 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 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 = ""; - application.createModule(modName, mesPar); - } - } - -/******************************************************************************* - * 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(modName, weakPar);\ - }\ -} -HW_CONTRACTION(Eye) // weakContractionEye -HW_CONTRACTION(NonEye) // weakContractionNonEye - -int main(int argc, char *argv[]) -{ - // parse command line ////////////////////////////////////////////////////// - std::string configStem; - - if (argc < 2) - { - std::cerr << "usage: " << argv[0] << " [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; - unsigned int light = 0; - unsigned int strange = 1; - unsigned int charm = 2; - std::vector mass = {.01, .04, .2}; - std::vector flavour = {"l", "s", "c"}; - std::vector 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 tKs = {0, 16}; - unsigned int dt_pi = 16; - std::vector tJs = {8}; - unsigned int n_noise = 0; - unsigned int nt = 32; - - // 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 - MGauge::Load::Par gaugePar; - gaugePar.file = configStem; - application.createModule("gauge", gaugePar); - for (unsigned int i = 0; i < flavour.size(); ++i) - { - // actions - MAction::DWF::Par actionPar; - actionPar.gauge = "gauge"; - actionPar.Ls = 12; - actionPar.M5 = 1.8; - actionPar.mass = mass[i]; - application.createModule("DWF_" + flavour[i], actionPar); - - // solvers - MSolver::RBPrecCG::Par solverPar; - solverPar.action = "DWF_" + flavour[i]; - solverPar.residual = 1.0e-8; - application.createModule(solvers[i], - solverPar); - } - - // Create noise propagators for loops. - if (n_noise > 0) - { - std::vector noiseProps_l(n_noise); - std::vector noiseProps_c(n_noise); - MSource::Z2::Par noisePar; - noisePar.tA = 0; - noisePar.tB = nt - 1; - for (unsigned int nn = 0; nn < n_noise; ++nn) - { - std::string eta = INIT_INDEX("noise", nn); - std::string loop_l = INIT_INDEX("loop_l", nn); - std::string loop_c = INIT_INDEX("loop_c", nn); - - application.createModule(eta, noisePar); - makePropagator(application, loop_l, eta, solvers[light]); - makePropagator(application, loop_c, eta, solvers[charm]); - noiseProps_l.push_back(loop_l); - noiseProps_c.push_back(loop_c); - } - } - - // 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]); - - // Wall sink-smeared propagators. - std::string q_Kl_0_W = q_Kl_0 + "_W"; - std::string q_Ks_k_W = q_Ks_k + "_W"; - std::string q_Ks_p_W = q_Ks_p + "_W"; - std::string q_pil_0_W = q_pil_0 + "_W"; - std::string q_pil_k_W = q_pil_k + "_W"; - std::string q_pil_p_W = q_pil_p + "_W"; - makeWallSink(application, q_Kl_0, q_Kl_0_W); - makeWallSink(application, q_Ks_k, q_Ks_k_W, kmom); - makeWallSink(application, q_Ks_p, q_Ks_p_W, pmom); - makeWallSink(application, q_pil_0, q_pil_0_W); - makeWallSink(application, q_pil_k, q_pil_k_W, kmom); - makeWallSink(application, q_pil_p, q_pil_p_W, pmom); - - /*********************************************************************** - * 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. - **********************************************************************/ - - for (unsigned int nn = 0; nn < n_noise; ++nn) - { - /******************************************************************* - * CONTRACTIONS: 3pt Weak Hamiltonian, S and E (Eye type) classes. - ******************************************************************/ - - } - - // Perform separate contractions for each t_J position. - for (unsigned int i = 0; i < tJs.size(); ++i) - { - // Sequential sources for current insertions. Local for now, - // gamma_0 only. - unsigned int tJ = (tJs[i] + 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. - ******************************************************************/ - - // Sequential sources for each noise propagator. - for (unsigned int nn = 0; nn < n_noise; ++nn) - { - /*************************************************************** - * CONTRACTIONS: 4pt contractions, S & E classes. - **************************************************************/ - - /*************************************************************** - * CONTRACTIONS: 4pt contractions, pi0 disconnected loop. - **************************************************************/ - } - } - } - // execution - application.saveParameterFile("rarekaon_000_100_tK0_tpi16_tJ8_noloop_mc0.2.xml"); - application.run(); - - // epilogue - LOG(Message) << "Grid is finalizing now" << std::endl; - Grid_finalize(); - - return EXIT_SUCCESS; -} From d7464aa0fe95bf286c8b7e0e88c309139f8b10f1 Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Wed, 8 Feb 2017 16:13:44 +0000 Subject: [PATCH 36/36] Switched from XmlWriter to CorrWriter in contraction code --- extras/Hadrons/Modules/MContraction/Baryon.hpp | 2 +- extras/Hadrons/Modules/MContraction/Meson.hpp | 2 +- extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc | 2 +- extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/extras/Hadrons/Modules/MContraction/Baryon.hpp b/extras/Hadrons/Modules/MContraction/Baryon.hpp index 6bccda2e..be7d919c 100644 --- a/extras/Hadrons/Modules/MContraction/Baryon.hpp +++ b/extras/Hadrons/Modules/MContraction/Baryon.hpp @@ -112,7 +112,7 @@ void TBaryon::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(par().q1); PropagatorField2 &q2 = *env().template getObject(par().q2); PropagatorField3 &q3 = *env().template getObject(par().q2); diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index 9581b8d5..4cbe1ac4 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -162,7 +162,7 @@ void TMeson::execute(void) << " quarks '" << par().q1 << "' and '" << par().q2 << "'" << std::endl; - XmlWriter writer(par().output); + CorrWriter writer(par().output); PropagatorField1 &q1 = *env().template getObject(par().q1); PropagatorField2 &q2 = *env().template getObject(par().q2); LatticeComplex c(env().getGrid()); diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index 4a1c9e39..a44c2534 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -93,7 +93,7 @@ void TWeakHamiltonianEye::execute(void) << par().q2 << ", '" << par().q3 << "' and '" << par().q4 << "'." << std::endl; - XmlWriter writer(par().output); + 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); diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc index 7945f9cb..2c4df68a 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -95,7 +95,7 @@ void TWeakHamiltonianNonEye::execute(void) << par().q2 << ", '" << par().q3 << "' and '" << par().q4 << "'." << std::endl; - XmlWriter writer(par().output); + 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);