From 07f2ebea1b8e3644d84af15e2be5649d0c34d58a Mon Sep 17 00:00:00 2001 From: Lanny91 Date: Thu, 19 Jan 2017 22:18:42 +0000 Subject: [PATCH] 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);