diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index b08a7468..72807d8b 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -7,5 +7,6 @@ #include #include #include +#include #include #include diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp new file mode 100644 index 00000000..e3dad26a --- /dev/null +++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp @@ -0,0 +1,137 @@ +#ifndef Hadrons_SeqGamma_hpp_ +#define Hadrons_SeqGamma_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Sequential source + ----------------------------- + * src_x = q_x * theta(x_3 - tA) * theta(tB - x_3) * gamma * exp(i x.mom) + + * options: + - q: input propagator (string) + - tA: begin timeslice (integer) + - tB: end timesilce (integer) + - gamma: gamma product to insert (integer) + - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + + */ + +/****************************************************************************** + * SeqGamma * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class SeqGammaPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(SeqGammaPar, + std::string, q, + unsigned int, tA, + unsigned int, tB, + unsigned int, gamma, + std::string, mom); +}; + +template +class TSeqGamma: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TSeqGamma(const std::string name); + // destructor + virtual ~TSeqGamma(void) = default; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +typedef TSeqGamma SeqGamma; + +/****************************************************************************** + * TSeqGamma implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TSeqGamma::TSeqGamma(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TSeqGamma::getInput(void) +{ + std::vector in = {par().q}; + + return in; +} + +template +std::vector TSeqGamma::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TSeqGamma::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TSeqGamma::execute(void) +{ + if (par().tA == par().tB) + { + LOG(Message) << "Generating gamma_" << par().gamma + << " sequential source at t= " << par().tA << std::endl; + } + else + { + LOG(Message) << "Generating gamma_" << par().gamma + << " sequential source for " + << par().tA << " <= t <= " << par().tB << std::endl; + } + PropagatorField &src = *env().template createLattice(getName()); + PropagatorField &q = *env().template getObject(par().q); + Lattice> t(env().getGrid()); + LatticeComplex ph(env().getGrid()), coor(env().getGrid()); + SpinMatrix g; + std::vector p; + Complex i(0.0,1.0); + + g = makeGammaProd(par().gamma); + p = strToVec(par().mom); + ph = zero; + for(unsigned int mu = 0; mu < Nd; mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + p[mu]*coor; + } + ph = exp(i*ph); + LatticeCoordinate(t, Tp); + src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q); +} + +END_MODULE_NAMESPACE + +MODULE_REGISTER_NS(SeqGamma, MSource); + +END_HADRONS_NAMESPACE + +#endif // Hadrons_SeqGamma_hpp_ diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index aad6c73b..4251ffa3 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -13,6 +13,7 @@ modules_hpp =\ Modules/MGauge/Unit.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ + Modules/MSource/SeqGamma.hpp \ Modules/MSource/Z2.hpp \ Modules/Quark.hpp diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc new file mode 100644 index 00000000..159981d0 --- /dev/null +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -0,0 +1,171 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_meson_3pt.cc + + Copyright (C) 2015 + + Author: Antonin Portelli + + 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 + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + + std::vector flavour = {"l", "s", "c1", "c2", "c3"}; + std::vector mass = {.01, .04, .2 , .25 , .3 }; + unsigned int nt = GridDefaultLatt()[Tp]; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.seed = "1 2 3 4"; + globalPar.genetic.maxGen = 5000; + globalPar.genetic.maxCstGen = 100; + globalPar.genetic.popSize = 50; + globalPar.genetic.mutationRate = .25; + application.setPar(globalPar); + + // gauge field + application.createModule("gauge"); + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 12; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + application.createModule("DWF_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule("CG_" + flavour[i], + solverPar); + } + for (unsigned int t = 0; t < nt; t += 8) + { + std::string srcName; + std::vector qName; + std::vector> seqName; + + // Z2 source + MSource::Z2::Par z2Par; + z2Par.tA = t; + z2Par.tB = t; + srcName = "z2_" + std::to_string(t); + application.createModule(srcName, z2Par); + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // sequential sources + MSource::SeqGamma::Par seqPar; + qName.push_back("QZ2_" + flavour[i] + srcName); + seqPar.q = qName[i]; + seqPar.tA = (t + nt/4) % nt; + seqPar.tB = (t + nt/4) % nt; + seqPar.mom = "1. 0. 0. 0."; + seqName.push_back(std::vector(Nd)); + for (unsigned int mu = 0; mu < Nd; ++mu) + { + seqPar.gamma = 0x1 << mu; + seqName[i][mu] = "G" + std::to_string(seqPar.gamma) + + "_" + std::to_string(seqPar.tA) + "-" + + qName[i]; + application.createModule(seqName[i][mu], seqPar); + } + + // propagators + Quark::Par quarkPar; + quarkPar.solver = "CG_" + flavour[i]; + quarkPar.source = srcName; + application.createModule(qName[i], quarkPar); + for (unsigned int mu = 0; mu < Nd; ++mu) + { + quarkPar.source = seqName[i][mu]; + seqName[i][mu] = qName[i] + "-" + seqName[i][mu]; + application.createModule(seqName[i][mu], quarkPar); + } + } + + // contractions + MContraction::Meson::Par mesPar; + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = i; j < flavour.size(); ++j) + { + mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; + mesPar.q1 = qName[i]; + mesPar.q2 = qName[j]; + application.createModule("meson_Z2_" + + std::to_string(t) + + "_" + + flavour[i] + + flavour[j], + mesPar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + for (unsigned int j = 0; j < flavour.size(); ++j) + for (unsigned int mu = 0; mu < Nd; ++mu) + { + MContraction::Meson::Par mesPar; + + mesPar.output = "3pt/Z2_" + flavour[i] + flavour[j] + "_" + + std::to_string(mu); + mesPar.q1 = qName[i]; + mesPar.q2 = seqName[j][mu]; + application.createModule("3pt_Z2_" + + std::to_string(t) + + "_" + + flavour[i] + + flavour[j] + + "_" + + std::to_string(mu), + mesPar); + } + } + + Environment::getInstance().printContent(); + // execution + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +}