Commiting reation of meson field code before a merge with the upstream branch feature/hadrons

extras/Hadrons/Modules.hpp

@@ -29,13 +29,14 @@
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
 #include <Grid/Hadrons/Modules/MContraction/DiscLoop.hpp>
 #include <Grid/Hadrons/Modules/MContraction/Meson.hpp>
+#include <Grid/Hadrons/Modules/MContraction/MesonFieldGmu.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WardIdentity.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp>
 #include <Grid/Hadrons/Modules/MContraction/Gamma3pt.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp>
+#include <Grid/Hadrons/Modules/MContraction/MesonFieldGmu.hpp>
 #include <Grid/Hadrons/Modules/MContraction/Baryon.hpp>
 #include <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>
-#include <Grid/Hadrons/Modules/MContraction/A2AMeson.hpp>
 #include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
 #include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
 #include <Grid/Hadrons/Modules/MScalar/ChargedProp.hpp>

extras/Hadrons/Modules/MContraction/A2AMeson.cc

@@ -1,7 +0,0 @@
-#include <Grid/Hadrons/Modules/MContraction/A2AMeson.hpp>
-
-using namespace Grid;
-using namespace Hadrons;
-using namespace MContraction;
-
-template class Grid::Hadrons::MContraction::TA2AMeson<FIMPL>;

extras/Hadrons/Modules/MContraction/MesonFieldGmu.cc

@@ -0,0 +1,7 @@
+#include <Grid/Hadrons/Modules/MContraction/MesonFieldGmu.hpp>
+
+using namespace Grid;
+using namespace Hadrons;
+using namespace MContraction;
+
+template class Grid::Hadrons::MContraction::TMesonFieldGmu<FIMPL>;

extras/Hadrons/Modules/MContraction/MesonFieldGmu.hpp

@@ -13,12 +13,10 @@ BEGIN_HADRONS_NAMESPACE
  ******************************************************************************/
 BEGIN_MODULE_NAMESPACE(MContraction)
 
-typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaPair;
+class MesonFieldPar : Serializable
-
-class A2AMesonPar : Serializable
 {
   public:
-    GRID_SERIALIZABLE_CLASS_MEMBERS(A2AMesonPar,
+    GRID_SERIALIZABLE_CLASS_MEMBERS(MesonFieldPar,
                                     int, Nl,
                                     int, N,
                                     std::string, A2A1,
@@ -26,12 +24,12 @@ class A2AMesonPar : Serializable
                                     std::string, action,
                                     std::string, epack1,
                                     std::string, epack2,
-                                    std::string, gammas,
+                                    std::string, gamma,
                                     std::string, output);
 };
 
 template <typename FImpl>
-class TA2AMeson : public Module<A2AMesonPar>
+class TMesonFieldGmu : public Module<MesonFieldPar>
 {
   public:
     FERM_TYPE_ALIASES(FImpl, );
@@ -42,41 +40,39 @@ class TA2AMeson : public Module<A2AMesonPar>
     {
       public:
         GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
-                                        Gamma::Algebra, gamma_snk,
+                                        Gamma::Algebra, gamma_mu,
-                                        Gamma::Algebra, gamma_src,
+                                        std::vector<std::vector<std::vector<ComplexD>>>, MesonField);
-                                        std::vector<Complex>, corr);
     };
 
   public:
     // constructor
-    TA2AMeson(const std::string name);
+    TMesonFieldGmu(const std::string name);
     // destructor
-    virtual ~TA2AMeson(void){};
+    virtual ~TMesonFieldGmu(void){};
     // dependency relation
     virtual std::vector<std::string> getInput(void);
     virtual std::vector<std::string> getOutput(void);
-    virtual void parseGammaString(std::vector<GammaPair> &gammaList);
     // setup
     virtual void setup(void);
     // execution
     virtual void execute(void);
 };
 
-MODULE_REGISTER(A2AMeson, ARG(TA2AMeson<FIMPL>), MContraction);
+MODULE_REGISTER(A2AMeson, ARG(TMesonFieldGmu<FIMPL>), MContraction);
 
 /******************************************************************************
-*                  TA2AMeson implementation                             *
+*                  TMesonFieldGmu implementation                             *
 ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
 template <typename FImpl>
-TA2AMeson<FImpl>::TA2AMeson(const std::string name)
+TMesonFieldGmu<FImpl>::TMesonFieldGmu(const std::string name)
-    : Module<A2AMesonPar>(name)
+    : Module<MesonFieldPar>(name)
 {
 }
 
 // dependencies/products ///////////////////////////////////////////////////////
 template <typename FImpl>
-std::vector<std::string> TA2AMeson<FImpl>::getInput(void)
+std::vector<std::string> TMesonFieldGmu<FImpl>::getInput(void)
 {
     std::vector<std::string> in = {par().A2A1, par().A2A2, par().action};
     in.push_back(par().A2A1 + "_ret");
@@ -92,117 +88,72 @@ std::vector<std::string> TA2AMeson<FImpl>::getInput(void)
 }
 
 template <typename FImpl>
-std::vector<std::string> TA2AMeson<FImpl>::getOutput(void)
+std::vector<std::string> TMesonFieldGmu<FImpl>::getOutput(void)
 {
     std::vector<std::string> out = {};
 
     return out;
 }
 
-template <typename FImpl>
-void TA2AMeson<FImpl>::parseGammaString(std::vector<GammaPair> &gammaList)
-{
-    gammaList.clear();
-    // Parse individual contractions from input string.
-    gammaList = strToVec<GammaPair>(par().gammas);
-}
-
 // setup ///////////////////////////////////////////////////////////////////////
 template <typename FImpl>
-void TA2AMeson<FImpl>::setup(void)
+void TMesonFieldGmu<FImpl>::setup(void)
 {
     int nt = env().getDim(Tp);
     int N = par().N;
 
     int Ls_ = env().getObjectLs(par().A2A1 + "_ret");
 
-    envTmp(std::vector<FermionField>, "w1", 1, N, FermionField(env().getGrid(1)));
+    envTmpLat(FermionField, "w", Ls_);
-    envTmp(std::vector<FermionField>, "v1", 1, N, FermionField(env().getGrid(1)));
+    envTmpLat(FermionField, "v", Ls_);
     envTmpLat(FermionField, "tmpv_5d", Ls_);
     envTmpLat(FermionField, "tmpw_5d", Ls_);
-
-    envTmp(std::vector<ComplexD>, "MF_x", 1, nt);
-    envTmp(std::vector<ComplexD>, "MF_y", 1, nt);
-    envTmp(std::vector<ComplexD>, "tmp", 1, nt);
 }
 
 // execution ///////////////////////////////////////////////////////////////////
 template <typename FImpl>
-void TA2AMeson<FImpl>::execute(void)
+void TMesonFieldGmu<FImpl>::execute(void)
 {
-    LOG(Message) << "Computing A2A meson contractions" << std::endl;
+    LOG(Message) << "Computing A2A meson field for gamma_mu  = " << par().gamma << ", taking w from " << par().A2A1 << " and v from " << par().A2A2 << std::endl;
 
     Result result;
-    Gamma g5(Gamma::Algebra::Gamma5);
+    std::istringstream sstr(par().gamma);
-    std::vector<GammaPair> gammaList;
+    Gamma::Algebra g_mu;
-    int nt = env().getDim(Tp);
+    sstr >> g_mu;
+    Gamma gamma_mu(g_mu);
 
-    parseGammaString(gammaList);
-
-    result.gamma_snk = gammaList[0].first;
-    result.gamma_src = gammaList[0].second;
-    result.corr.resize(nt);
-
-    int Nl = par().Nl;
     int N = par().N;
     LOG(Message) << "N for A2A cont: " << N << std::endl;
+    result.gamma_mu = g_mu;
+    result.MesonField.resize(N);
 
-    envGetTmp(std::vector<ComplexD>, MF_x);
+    int nt = env().getDim(Tp);
-    envGetTmp(std::vector<ComplexD>, MF_y);
+    std::vector<ComplexD> MesonField_ij;
-    envGetTmp(std::vector<ComplexD>, tmp);
+    MesonField_ij.resize(nt);
 
-    for (unsigned int t = 0; t < nt; ++t)
+    auto &a2a1 = envGet(A2ABase, par().A2A1 + "_ret");
-    {
+    auto &a2a2 = envGet(A2ABase, par().A2A2 + "_ret");
-        tmp[t] = TensorRemove(MF_x[t] * MF_y[t] * 0.0);
-    }
 
-    Gamma gSnk(gammaList[0].first);
+    envGetTmp(FermionField, w);
-    Gamma gSrc(gammaList[0].second);
+    envGetTmp(FermionField, v);
-
-    auto &a2a1_fn = envGet(A2ABase, par().A2A1 + "_ret");
-
-    envGetTmp(std::vector<FermionField>, w1);
-    envGetTmp(std::vector<FermionField>, v1);
     envGetTmp(FermionField, tmpv_5d);
     envGetTmp(FermionField, tmpw_5d);
 
-    LOG(Message) << "Finding v and w vectors for N =  " << N << std::endl;
-    for (int i = 0; i < N; i++)
-    {
-        a2a1_fn.return_v(i, tmpv_5d, v1[i]);
-        a2a1_fn.return_w(i, tmpw_5d, w1[i]);
-    }
-    LOG(Message) << "Found v and w vectors for N =  " << N << std::endl;
-    for (unsigned int i = 0; i < N; i++)
-    {
-        v1[i] = gSnk * v1[i];
-    }
-    int ty;
     for (unsigned int i = 0; i < N; i++)
     {
+        a2a1.return_w(i, tmpw_5d, w);
         for (unsigned int j = 0; j < N; j++)
         {
-            sliceInnerProductVector(MF_x, w1[i], v1[j], Tp);
+            a2a2.return_v(j, tmpv_5d, v);
-            sliceInnerProductVector(MF_y, w1[j], v1[i], Tp);
+            v = gamma_mu*v;
-            for (unsigned int t = 0; t < nt; ++t)
+            sliceInnerProductVector(MesonField_ij, w, v, Tp);
-            {
+            result.MesonField[j][i] = MesonField_ij;
-                for (unsigned int tx = 0; tx < nt; tx++)
-                {
-                    ty = (tx + t) % nt;
-                    tmp[t] += TensorRemove((MF_x[tx]) * (MF_y[ty]));
-                }
-            }
         }
         if (i % 10 == 0)
         {
             LOG(Message) << "MF for i = " << i << " of " << N << std::endl;
         }
     }
-    double NTinv = 1.0 / static_cast<double>(nt);
-    for (unsigned int t = 0; t < nt; ++t)
-    {
-        result.corr[t] = NTinv*tmp[t];
-    }
 
     saveResult(par().output, "meson", result);
 }

extras/Hadrons/modules.inc

@@ -30,8 +30,8 @@ modules_cc =\
   Modules/MContraction/WardIdentity.cc \
   Modules/MContraction/WeakHamiltonianEye.cc \
   Modules/MContraction/DiscLoop.cc \
-  Modules/MContraction/A2AMeson.cc \
   Modules/MContraction/Baryon.cc \
+  Modules/MContraction/MesonFieldGmu.cc \
   Modules/MContraction/Gamma3pt.cc \
   Modules/MContraction/WeakNeutral4ptDisc.cc \
   Modules/MContraction/Meson.cc \
@@ -86,9 +86,9 @@ modules_hpp =\
   Modules/MContraction/WeakHamiltonianEye.hpp \
   Modules/MContraction/Gamma3pt.hpp \
   Modules/MContraction/WeakHamiltonianNonEye.hpp \
+  Modules/MContraction/MesonFieldGmu.hpp \
   Modules/MContraction/Baryon.hpp \
   Modules/MContraction/WeakNeutral4ptDisc.hpp \
-  Modules/MContraction/A2AMeson.hpp \
   Modules/MScalar/Scalar.hpp \
   Modules/MScalar/FreeProp.hpp \
   Modules/MScalar/ChargedProp.hpp \