diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 6e1b8823..d0d0d80d 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include diff --git a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp index 7f643d49..162ab786 100644 --- a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp +++ b/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::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(par().q1); + SlicedPropagator1 &q1 = *env().template getObject(par().q1); PropagatorField2 &q2 = *env().template getObject(par().q2); - PropagatorField3 &q3 = *env().template getObject(par().q3); + PropagatorField3 &q3 = *env().template getObject(par().q3); LatticeComplex c(env().getGrid()); Gamma g5(Gamma::Algebra::Gamma5); Gamma gamma(par().gamma); std::vector 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; diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp index 0a3c2e31..302b207e 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp +++ b/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); }; diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc index a44c2534..314b080a 100644 --- a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc +++ b/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(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(n_eye_diag); - unsigned int ndim = env().getNd(); + SlicedPropagator &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(n_eye_diag); + unsigned int ndim = env().getNd(); PropagatorField tmp1(env().getGrid()); LatticeComplex tmp2(env().getGrid()); @@ -111,10 +113,13 @@ void TWeakHamiltonianEye::execute(void) std::vector E_body(ndim, tmp2); std::vector 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])); } diff --git a/extras/Hadrons/Modules/MSink/Smear.hpp b/extras/Hadrons/Modules/MSink/Smear.hpp new file mode 100644 index 00000000..9327001f --- /dev/null +++ b/extras/Hadrons/Modules/MSink/Smear.hpp @@ -0,0 +1,99 @@ +#ifndef Hadrons_MSink_Smear_hpp_ +#define Hadrons_MSink_Smear_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Smear * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSink) + +class SmearPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(SmearPar, + std::string, q, + std::string, sink); +}; + +template +class TSmear: public Module +{ +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 getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_NS(Smear, TSmear, MSink); + +/****************************************************************************** + * TSmear implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TSmear::TSmear(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TSmear::getInput(void) +{ + std::vector in = {par().q, par().sink}; + + return in; +} + +template +std::vector TSmear::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TSmear::setup(void) +{ + unsigned int nt = env().getDim(Tp); + unsigned int size = nt * sizeof(SitePropagator); + env().registerObject(getName(), size); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TSmear::execute(void) +{ + LOG(Message) << "Sink smearing propagator '" << par().q + << "' using sink function '" << par().sink << "'." + << std::endl; + + SinkFn &sink = *env().template getObject(par().sink); + PropagatorField &q = *env().template getObject(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_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 91d0bbe1..fbbb2eb9 100644 --- a/extras/Hadrons/modules.inc +++ b/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 \ diff --git a/tests/hadrons/Test_hadrons.hpp b/tests/hadrons/Test_hadrons.hpp index 6ea51d72..9bd3ee0a 100644 --- a/tests/hadrons/Test_hadrons.hpp +++ b/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(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 - * Purpose: Wall sink smearing of a propagator. + * Name: makePointSink + * 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, - std::string &wallName, std::string mom = ZERO_MOM) +inline void makePointSink(Application &application, std::string &sinkFnct, + 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(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(wallName))) + if (!(Environment::getInstance().hasModule(smearedProp))) { - MSink::Wall::Par wallPar; - wallPar.q = propName; - wallPar.mom = mom; - application.createModule(wallName, wallPar); - }*/ + MSink::Smear::Par smearPar; + smearPar.q = propName; + smearPar.sink = sinkFnct; + application.createModule(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(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(modName, weakPar);\ }\ } diff --git a/tests/hadrons/Test_hadrons_3pt_contractions.cc b/tests/hadrons/Test_hadrons_3pt_contractions.cc new file mode 100644 index 00000000..452fc34d --- /dev/null +++ b/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 + + 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(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; +} \ No newline at end of file diff --git a/tests/hadrons/Test_hadrons_quark.cc b/tests/hadrons/Test_hadrons_quark.cc index 5b9d0ce1..eac065e9 100644 --- a/tests/hadrons/Test_hadrons_quark.cc +++ b/tests/hadrons/Test_hadrons_quark.cc @@ -26,7 +26,7 @@ *******************************************************************************/ #include "Test_hadrons.hpp" -#include +#include using namespace Grid; using namespace QCD;