Merge branch 'feature/rare_kaon' of https://github.com/Lanny91/Grid into feature/hadrons

extras/Hadrons/Modules.hpp

@@ -19,6 +19,7 @@
 #include <Grid/Hadrons/Modules/MScalar/FreeProp.hpp>
 #include <Grid/Hadrons/Modules/MScalar/Scalar.hpp>
 #include <Grid/Hadrons/Modules/MSink/Point.hpp>
+#include <Grid/Hadrons/Modules/MSink/Smear.hpp>
 #include <Grid/Hadrons/Modules/MSolver/RBPrecCG.hpp>
 #include <Grid/Hadrons/Modules/MSource/Point.hpp>
 #include <Grid/Hadrons/Modules/MSource/SeqConserved.hpp>

extras/Hadrons/Modules/MContraction/Gamma3pt.hpp

@@ -51,6 +51,14 @@ BEGIN_HADRONS_NAMESPACE
  *                   q1
  *
  *      trace(g5*q1*adj(q2)*g5*gamma*q3)
+ * 
+ *  options:
+ *   - q1: sink smeared propagator, source at i
+ *   - q2: propagator, source at i
+ *   - q3: propagator, source at f
+ *   - gamma: gamma matrix to insert
+ *   - tSnk: sink position for propagator q1.
+ *
  */
 
 /******************************************************************************
@@ -66,6 +74,7 @@ public:
                                     std::string,    q2,
                                     std::string,    q3,
                                     Gamma::Algebra, gamma,
+                                    unsigned int,   tSnk,
                                     std::string,    output);
 };
 
@@ -140,17 +149,22 @@ void TGamma3pt<FImpl1, FImpl2, FImpl3>::execute(void)
                  << par().q3 << "', with " << par().gamma << " insertion." 
                  << std::endl;
 
+    // Initialise variables. q2 and q3 are normal propagators, q1 may be 
+    // sink smeared.
     CorrWriter            writer(par().output);
-    PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
+    SlicedPropagator1     &q1 = *env().template getObject<SlicedPropagator1>(par().q1);
     PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
-    PropagatorField3      &q3 = *env().template getObject<PropagatorField3>(par().q3);
+    PropagatorField3      &q3 = *env().template getObject<PropagatorField2>(par().q3);
     LatticeComplex        c(env().getGrid());
     Gamma                 g5(Gamma::Algebra::Gamma5);
     Gamma                 gamma(par().gamma);
     std::vector<TComplex> buf;
     Result                result;
-
+    
-    c = trace(g5*q1*adj(q2)*(g5*gamma)*q3);
+    // Extract relevant timeslice of sinked propagator q1, then contract &
+    // sum over all spacial positions of gamma insertion.
+    SitePropagator1 q1Snk = q1[par().tSnk];
+    c = trace(g5*q1Snk*adj(q2)*(g5*gamma)*q3);
     sliceSum(c, buf, Tp);
 
     result.gamma = par().gamma;

extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp

@@ -76,6 +76,7 @@ public:
                                     std::string, q2,
                                     std::string, q3,
                                     std::string, q4,
+                                    unsigned int, tSnk,
                                     std::string, output);
 };
 

extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc

@@ -54,6 +54,8 @@ using namespace MContraction;
  * 
  * 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])
+ * 
+ * Note q1 must be sink smeared.
  */
 
 /******************************************************************************
@@ -94,15 +96,15 @@ void TWeakHamiltonianEye::execute(void)
                  << "'." << std::endl;
 
     CorrWriter             writer(par().output);
-    PropagatorField &q1 = *env().template getObject<PropagatorField>(par().q1);
+    SlicedPropagator &q1 = *env().template getObject<SlicedPropagator>(par().q1);
-    PropagatorField &q2 = *env().template getObject<PropagatorField>(par().q2);
+    PropagatorField  &q2 = *env().template getObject<PropagatorField>(par().q2);
-    PropagatorField &q3 = *env().template getObject<PropagatorField>(par().q3);
+    PropagatorField  &q3 = *env().template getObject<PropagatorField>(par().q3);
-    PropagatorField &q4 = *env().template getObject<PropagatorField>(par().q4);
+    PropagatorField  &q4 = *env().template getObject<PropagatorField>(par().q4);
-    Gamma g5            = Gamma(Gamma::Algebra::Gamma5);
+    Gamma g5             = Gamma(Gamma::Algebra::Gamma5);
-    LatticeComplex        expbuf(env().getGrid());
+    LatticeComplex         expbuf(env().getGrid());
-    std::vector<TComplex> corrbuf;
+    std::vector<TComplex>  corrbuf;
-    std::vector<Result>   result(n_eye_diag);
+    std::vector<Result>    result(n_eye_diag);
-    unsigned int ndim   = env().getNd();
+    unsigned int ndim    = env().getNd();
 
     PropagatorField              tmp1(env().getGrid());
     LatticeComplex               tmp2(env().getGrid());
@@ -111,10 +113,13 @@ void TWeakHamiltonianEye::execute(void)
     std::vector<LatticeComplex>  E_body(ndim, tmp2);
     std::vector<LatticeComplex>  E_loop(ndim, tmp2);
 
+    // Get sink timeslice of q1.
+    SitePropagator q1Snk = q1[par().tSnk];
+
     // Setup for S-type contractions.
     for (int mu = 0; mu < ndim; ++mu)
     {
-        S_body[mu] = MAKE_SE_BODY(q1, q2, q3, GammaL(Gamma::gmu[mu]));
+        S_body[mu] = MAKE_SE_BODY(q1Snk, q2, q3, GammaL(Gamma::gmu[mu]));
         S_loop[mu] = MAKE_SE_LOOP(q4, GammaL(Gamma::gmu[mu]));
     }
 

extras/Hadrons/Modules/MSink/Smear.hpp

@@ -0,0 +1,99 @@
+#ifndef Hadrons_MSink_Smear_hpp_
+#define Hadrons_MSink_Smear_hpp_
+
+#include <Grid/Hadrons/Global.hpp>
+#include <Grid/Hadrons/Module.hpp>
+#include <Grid/Hadrons/ModuleFactory.hpp>
+
+BEGIN_HADRONS_NAMESPACE
+
+/******************************************************************************
+ *                                 Smear                                      *
+ ******************************************************************************/
+BEGIN_MODULE_NAMESPACE(MSink)
+
+class SmearPar: Serializable
+{
+public:
+    GRID_SERIALIZABLE_CLASS_MEMBERS(SmearPar,
+                                    std::string, q,
+                                    std::string, sink);
+};
+
+template <typename FImpl>
+class TSmear: public Module<SmearPar>
+{
+public:
+    FERM_TYPE_ALIASES(FImpl,);
+    SINK_TYPE_ALIASES();
+public:
+    // constructor
+    TSmear(const std::string name);
+    // destructor
+    virtual ~TSmear(void) = default;
+    // dependency relation
+    virtual std::vector<std::string> getInput(void);
+    virtual std::vector<std::string> getOutput(void);
+    // setup
+    virtual void setup(void);
+    // execution
+    virtual void execute(void);
+};
+
+MODULE_REGISTER_NS(Smear, TSmear<FIMPL>, MSink);
+
+/******************************************************************************
+ *                          TSmear implementation                             *
+ ******************************************************************************/
+// constructor /////////////////////////////////////////////////////////////////
+template <typename FImpl>
+TSmear<FImpl>::TSmear(const std::string name)
+: Module<SmearPar>(name)
+{}
+
+// dependencies/products ///////////////////////////////////////////////////////
+template <typename FImpl>
+std::vector<std::string> TSmear<FImpl>::getInput(void)
+{
+    std::vector<std::string> in = {par().q, par().sink};
+    
+    return in;
+}
+
+template <typename FImpl>
+std::vector<std::string> TSmear<FImpl>::getOutput(void)
+{
+    std::vector<std::string> out = {getName()};
+    
+    return out;
+}
+
+// setup ///////////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSmear<FImpl>::setup(void)
+{
+    unsigned int nt = env().getDim(Tp);
+    unsigned int size = nt * sizeof(SitePropagator);
+    env().registerObject(getName(), size);
+}
+
+// execution ///////////////////////////////////////////////////////////////////
+template <typename FImpl>
+void TSmear<FImpl>::execute(void)
+{
+    LOG(Message) << "Sink smearing propagator '" << par().q
+                 << "' using sink function '" << par().sink << "'."
+                 << std::endl;
+
+    SinkFn          &sink = *env().template getObject<SinkFn>(par().sink);
+    PropagatorField &q    = *env().template getObject<PropagatorField>(par().q);
+    SlicedPropagator *out = new SlicedPropagator(env().getDim(Tp));
+    *out  = sink(q);
+    env().setObject(getName(), out);
+}
+
+END_MODULE_NAMESPACE
+
+END_HADRONS_NAMESPACE
+
+#endif // Hadrons_MSink_Smear_hpp_

extras/Hadrons/modules.inc

@@ -31,6 +31,7 @@ modules_hpp =\
   Modules/MScalar/FreeProp.hpp \
   Modules/MScalar/Scalar.hpp \
   Modules/MSink/Point.hpp \
+  Modules/MSink/Smear.hpp \
   Modules/MSolver/RBPrecCG.hpp \
   Modules/MSource/Point.hpp \
   Modules/MSource/SeqConserved.hpp \

tests/hadrons/Test_hadrons.hpp

@@ -269,6 +269,26 @@ inline void makeConservedSequentialSource(Application &application,
     }
 }
 
+/*******************************************************************************
+ * Name: makeNoiseSource
+ * Parameters: application - main application that stores modules.
+ *             srcName     - name of source module to create.
+ *             tA          - lower source timeslice limit.
+ *             tB          - upper source timeslice limit.
+ * Returns: None.
+ ******************************************************************************/
+inline void makeNoiseSource(Application &application, std::string &srcName,
+                            unsigned int tA, unsigned int tB)
+{
+    if (!(Environment::getInstance().hasModule(srcName)))
+    {
+        MSource::Z2::Par noisePar;
+        noisePar.tA = tA;
+        noisePar.tB = tB;
+        application.createModule<MSource::Z2>(srcName, noisePar);
+    }
+ }
+
 /*******************************************************************************
  * Name: makeWallSource
  * Purpose: Construct wall source and add to application module.
@@ -292,26 +312,46 @@ inline void makeWallSource(Application &application, std::string &srcName,
 }
 
 /*******************************************************************************
- * Name: makeWallSink
+ * Name: makePointSink
- * Purpose: Wall sink smearing of a propagator.
+ * Purpose: Create function for point sink smearing of a propagator.
  * Parameters: application - main application that stores modules.
  *             propName    - name of input propagator.
- *             wallName    - name of smeared propagator.
+ *             sinkFnct    - name of output sink smearing module.
  *             mom         - momentum insertion (default is zero).
  * Returns: None.
  ******************************************************************************/
-inline void makeWallSink(Application &application, std::string &propName,
+inline void makePointSink(Application &application, std::string &sinkFnct,
-                         std::string &wallName, std::string mom = ZERO_MOM)
+                          std::string mom = ZERO_MOM)
+{
+    // If the sink function already exists, don't make it again.
+    if (!(Environment::getInstance().hasModule(sinkFnct)))
+    {
+        MSink::Point::Par pointPar;
+        pointPar.mom = mom;
+        application.createModule<MSink::Point>(sinkFnct, pointPar);
+    }
+}
+
+/*******************************************************************************
+ * Name: sinkSmear
+ * Purpose: Perform sink smearing of a propagator.
+ * Parameters: application - main application that stores modules.
+ *             sinkFnct    - sink smearing module.
+ *             propName    - propagator to smear.
+ *             smearedProp - name of output smeared propagator.
+ * Returns: None.
+ ******************************************************************************/
+inline void sinkSmear(Application &application, std::string &sinkFnct,
+                      std::string &propName, std::string &smearedProp)
 {
     // If the propagator has already been smeared, don't smear it again.
-    // Temporarily removed, strategy for sink smearing likely to change.
+    if (!(Environment::getInstance().hasModule(smearedProp)))
-    /*if (!(Environment::getInstance().hasModule(wallName)))
     {
-        MSink::Wall::Par wallPar;
+        MSink::Smear::Par smearPar;
-        wallPar.q   = propName;
+        smearPar.q    = propName;
-        wallPar.mom = mom;
+        smearPar.sink = sinkFnct;
-        application.createModule<MSink::Wall>(wallName, wallPar);
+        application.createModule<MSink::Smear>(smearedProp, smearPar);
-    }*/
+    }
 }
 
 /*******************************************************************************
@@ -398,16 +438,18 @@ inline void mesonContraction(Application &application,
  * Purpose: Create gamma3pt 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.
+ *             q1          - quark propagator 1, sink smeared.
  *             q2          - quark propagator 2.
  *             q3          - quark propagator 3.
  *             label       - unique label to construct module name.
+ *             tSnk        - sink position of sink for q1.
  *             gamma       - gamma insertions between q2 and q3.
  * Returns: None.
  ******************************************************************************/
 inline void gamma3ptContraction(Application &application, unsigned int npt, 
                                 std::string &q1, std::string &q2,
-                                std::string &q3, std::string &label, 
+                                std::string &q3, std::string &label,
+                                unsigned int tSnk = 0,
                                 Gamma::Algebra gamma = Gamma::Algebra::Identity)
 {
     std::string modName = std::to_string(npt) + "pt_" + label;
@@ -418,6 +460,7 @@ inline void gamma3ptContraction(Application &application, unsigned int npt,
         gamma3ptPar.q1 = q1;
         gamma3ptPar.q2 = q2;
         gamma3ptPar.q3 = q3;
+        gamma3ptPar.tSnk = tSnk;
         gamma3ptPar.gamma = gamma;
         application.createModule<MContraction::Gamma3pt>(modName, gamma3ptPar);
     }
@@ -434,13 +477,14 @@ inline void gamma3ptContraction(Application &application, unsigned int npt,
  *             q3          - quark propagator 3.
  *             q4          - quark propagator 4.
  *             label       - unique label to construct module name.
+ *             tSnk        - time position of sink (for sink smearing).
  * 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 &label, unsigned int tSnk = 0)\
 {\
     std::string modName = std::to_string(npt) + "pt_" + label;\
     if (!(Environment::getInstance().hasModule(modName)))\
@@ -451,6 +495,7 @@ inline void weakContraction##top(Application &application, unsigned int npt,\
         weakPar.q2 = q2;\
         weakPar.q3 = q3;\
         weakPar.q4 = q4;\
+        weakPar.tSnk = tSnk;\
         application.createModule<MContraction::WeakHamiltonian##top>(modName, weakPar);\
     }\
 }

tests/hadrons/Test_hadrons_3pt_contractions.cc

@@ -0,0 +1,122 @@
+/*******************************************************************************
+ Grid physics library, www.github.com/paboyle/Grid
+ 
+ Source file: tests/hadrons/Test_hadrons_3pt_contractions.cc
+ 
+ Copyright (C) 2017
+ 
+ Author: 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
+ 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 "Test_hadrons.hpp"
+
+using namespace Grid;
+using namespace Hadrons;
+
+int main(int argc, char *argv[])
+{
+    // initialization //////////////////////////////////////////////////////////
+    HADRONS_DEFAULT_INIT;
+
+    // run setup ///////////////////////////////////////////////////////////////
+    Application  application;
+    double       mass = 0.04;
+    double       M5   = 1.8;
+    unsigned int Ls   = 12;
+    unsigned int nt   = GridDefaultLatt()[Tp];
+    unsigned int t_i  = 0;
+    unsigned int t_f  = nt / 2;
+    std::string  mom  = "1. 0. 0. 0.";
+
+    // global parameters
+    HADRONS_DEFAULT_GLOBALS(application);
+
+    // gauge field
+    std::string gaugeField = "gauge";
+    application.createModule<MGauge::Unit>(gaugeField);
+
+    // Action & solver setup.
+    std::string action = "DWF";
+    std::string solver = "CG";
+    makeDWFAction(application, action, gaugeField, mass, M5, Ls);
+    makeRBPrecCGSolver(application, solver, action);
+
+    /***************************************************************************
+     * Weak Contraction test: Non-Eye class.
+     **************************************************************************/
+    // Make wall source propagators for each leg of 4-quark vertex.
+    std::string q_i_0 = "q_i_0";
+    std::string q_i_p = "q_i_p";
+    std::string q_f_0 = "q_f_0";
+    std::string q_f_p = "q_f_p";
+    MAKE_WALL_PROP(t_i, q_i_0, solver);
+    MAKE_WALL_PROP(t_f, q_f_0, solver);
+    MAKE_3MOM_WALL_PROP(t_i, mom, q_i_p, solver);
+    MAKE_3MOM_WALL_PROP(t_f, mom, q_f_p, solver);
+
+    // Perform contractions, zero and non-zero momentum.
+    std::string HW_CW_0 = LABEL_3PT("HW_CW_0", t_i, t_f);
+    std::string HW_CW_p = LABEL_3PT("HW_CW_p", t_i, t_f);
+    weakContractionNonEye(application, 3, q_i_0, q_i_0, q_f_0, q_f_0, HW_CW_0);
+    weakContractionNonEye(application, 3, q_i_0, q_i_p, q_f_p, q_f_0, HW_CW_p);
+
+    /***************************************************************************
+     * Weak Contraction test: Eye-class.
+     **************************************************************************/
+    // Create random propagator for loop.
+    std::string eta = "noise_source";
+    makeNoiseSource(application, eta, 0, nt - 1);
+    std::string loopProp = "loop";
+    std::string loopRes  = loopProp + "_res";
+    makePropagator(application, loopRes, eta, solver);
+    makeLoop(application, loopProp, eta, loopRes);
+
+    // Wall sink smear the propagator directly connecting the source & sink.
+    // (i.e. make point sink but smear before the contraction)
+    std::string wallSink = "wall_sink";
+    std::string qWall    = "q_wall";
+    makePointSink(application, wallSink);
+    sinkSmear(application, wallSink, q_i_0, qWall);
+
+    // Perform contractions, zero and non-zero momentum.
+    std::string HW_SE_0 = LABEL_3PT("HW_SE_0", t_i, t_f);
+    std::string HW_SE_p = LABEL_3PT("HW_SE_p", t_i, t_f);
+    weakContractionEye(application, 3, qWall, q_i_0, q_f_p, loopProp, HW_SE_0, t_f);
+    weakContractionEye(application, 3, qWall, q_i_p, q_f_p, loopProp, HW_SE_p, t_f);
+
+    /***************************************************************************
+     * Gamma insertion test.
+     **************************************************************************/
+    Gamma::Algebra gamma = Gamma::Algebra::GammaT;
+    std::string sd_0 = LABEL_3PT("sd_0", t_i, t_f);
+    std::string sd_p = LABEL_3PT("sd_p", t_i, t_f);
+    gamma3ptContraction(application, 3, qWall, q_i_0, q_f_0, sd_0, t_f, gamma);
+    gamma3ptContraction(application, 3, qWall, q_i_p, q_f_p, sd_p, t_f, gamma);
+
+    // execution
+    application.saveParameterFile("ContractionTest3pt.xml");
+    application.run();
+
+    // epilogue
+    LOG(Message) << "Grid is finalizing now" << std::endl;
+    Grid_finalize();
+
+    return EXIT_SUCCESS;
+}

tests/hadrons/Test_hadrons_quark.cc

@@ -26,7 +26,7 @@
  *******************************************************************************/
 
 #include "Test_hadrons.hpp"
-#include <Grid/Hadrons/Modules/Quark.hpp>
+#include <Grid/Hadrons/Modules/MFermion/GaugeProp.hpp>
 
 using namespace Grid;
 using namespace QCD;