Meson module now takes list of gamma matrices to insert at source and sink.

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -8,6 +8,7 @@ Copyright (C) 2015
 Copyright (C) 2016
 
 Author: Antonin Portelli <antonin.portelli@me.com>
+        Andrew Lawson    <andrew.lawson1991@gmail.com>
 
 This program is free software; you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -36,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<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
+    virtual void parseGammaString(SpinMatrix*,
+                                  unsigned int &, 
+                                  unsigned int &,
+                                  std::vector<std::vector<bool>> &);
     // execution
     virtual void execute(void);
 };
@@ -102,6 +124,79 @@ std::vector<std::string> TMeson<FImpl1, FImpl2>::getOutput(void)
     return output;
 }
 
+template <typename T>
+std::vector<std::pair<T, T>> strToVecPair(const std::string s)
+{
+    std::vector<std::pair<T, T>> v;
+    return v;
+}
+
+template <typename FImpl1, typename FImpl2>
+void TMeson<FImpl1, FImpl2>::parseGammaString(SpinMatrix *g,
+                                              unsigned int &n_snk, 
+                                              unsigned int &n_src,
+                                              std::vector<std::vector<bool>> &toDo)
+{
+    // Initialise counters for parsing gamma insertions at source & sink.
+    int empty = -1;
+    std::vector<int> 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<std::pair<int, int>> gamma_pairs;
+        gamma_pairs = strToVecPair<int>(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 <typename FImpl1, typename FImpl2>
 void TMeson<FImpl1, FImpl2>::execute(void)
@@ -115,26 +210,29 @@ void TMeson<FImpl1, FImpl2>::execute(void)
     PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
     LatticeComplex        c(env().getGrid());
     SpinMatrix            g[Ns*Ns], g5;
+    std::vector<std::vector<bool>> toDo;
     std::vector<TComplex> buf;
     Result                result;
-    
+    unsigned int          n_snk, n_src;
+
     g5 = makeGammaProd(Ns*Ns - 1);
-    result.corr.resize(Ns*Ns);
+    parseGammaString(g, n_snk, n_src, toDo);
-    for (unsigned int i = 0; i < Ns*Ns; ++i)
+
+    result.corr.resize(n_snk);
+    for (unsigned int iSink = 0; iSink < toDo.size(); ++iSink)
     {
-        g[i] = makeGammaProd(i);
+        result.corr[iSink].resize(n_src);
-    }
+        for (unsigned int iSrc = 0; iSrc < toDo.size(); ++iSrc)
-    for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink)
-    {
-        result.corr[iSink].resize(Ns*Ns);
-        for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc)
         {
-            c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5);
+            if (toDo[iSink][iSrc])
-            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]);
+                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]);
+                }
             }
         }
     }

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<MContraction::Meson>("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<MContraction::Meson>("3pt_Z2_"
                                                           + std::to_string(t)
                                                           + "_"

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<MContraction::Meson>("meson_pt_"
                                                       + flavour[i] + flavour[j],
                                                       mesPar);
         mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
         mesPar.q1     = "QZ2_" + flavour[i];
         mesPar.q2     = "QZ2_" + flavour[j];
+        mesPar.gammas = "all";
         application.createModule<MContraction::Meson>("meson_Z2_"
                                                       + flavour[i] + flavour[j],
                                                       mesPar);