diff --git a/Makefile.am b/Makefile.am index d3401c48..1a75f2eb 100644 --- a/Makefile.am +++ b/Makefile.am @@ -3,6 +3,8 @@ SUBDIRS = lib benchmarks tests extras include $(top_srcdir)/doxygen.inc +bin_SCRIPTS=grid-config + tests: all $(MAKE) -C tests tests diff --git a/benchmarks/Benchmark_dwf.cc b/benchmarks/Benchmark_dwf.cc index 61c06800..a071c050 100644 --- a/benchmarks/Benchmark_dwf.cc +++ b/benchmarks/Benchmark_dwf.cc @@ -145,6 +145,7 @@ int main (int argc, char ** argv) RealD M5 =1.8; RealD NP = UGrid->_Nprocessors; + RealD NN = UGrid->NodeCount(); std::cout << GridLogMessage<< "*****************************************************************" < grid.configure.summary + +GRID_SUMMARY="`cat grid.configure.summary`" +AM_SUBST_NOTMAKE([GRID_SUMMARY]) +AC_SUBST([GRID_SUMMARY]) + +AC_CONFIG_FILES([grid-config], [chmod +x grid-config]) +AC_CONFIG_FILES(Makefile) +AC_CONFIG_FILES(lib/Makefile) +AC_CONFIG_FILES(tests/Makefile) +AC_CONFIG_FILES(tests/IO/Makefile) +AC_CONFIG_FILES(tests/core/Makefile) +AC_CONFIG_FILES(tests/debug/Makefile) +AC_CONFIG_FILES(tests/forces/Makefile) +AC_CONFIG_FILES(tests/hadrons/Makefile) +AC_CONFIG_FILES(tests/hmc/Makefile) +AC_CONFIG_FILES(tests/solver/Makefile) +AC_CONFIG_FILES(tests/qdpxx/Makefile) +AC_CONFIG_FILES(tests/testu01/Makefile) +AC_CONFIG_FILES(benchmarks/Makefile) +AC_CONFIG_FILES(extras/Makefile) +AC_CONFIG_FILES(extras/Hadrons/Makefile) +AC_OUTPUT + echo "" cat grid.configure.summary echo "" diff --git a/extras/Hadrons/Application.cc b/extras/Hadrons/Application.cc index 62674f30..90ebcfd7 100644 --- a/extras/Hadrons/Application.cc +++ b/extras/Hadrons/Application.cc @@ -162,7 +162,8 @@ void Application::saveParameterFile(const std::string parameterFileName) sizeString((size)*locVol_) << " (" << sizeString(size) << "/site)" #define DEFINE_MEMPEAK \ -auto memPeak = [this](const std::vector &program)\ +GeneticScheduler::ObjFunc memPeak = \ +[this](const std::vector &program)\ {\ unsigned int memPeak;\ bool msg;\ diff --git a/extras/Hadrons/Global.hpp b/extras/Hadrons/Global.hpp index 81afab13..3e11ddf8 100644 --- a/extras/Hadrons/Global.hpp +++ b/extras/Hadrons/Global.hpp @@ -145,6 +145,15 @@ std::string typeName(void) return typeName(typeIdPt()); } +// default writers/readers +#ifdef HAVE_HDF5 +typedef Hdf5Reader CorrReader; +typedef Hdf5Writer CorrWriter; +#else +typedef XmlReader CorrReader; +typedef XmlWriter CorrWriter; +#endif + END_HADRONS_NAMESPACE #endif // Hadrons_Global_hpp_ diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 77ae08b7..05ad1697 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -29,12 +29,20 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include +#include #include +#include +#include +#include +#include #include #include #include +#include #include #include #include +#include #include #include diff --git a/extras/Hadrons/Modules/MContraction/Baryon.hpp b/extras/Hadrons/Modules/MContraction/Baryon.hpp index 6bccda2e..be7d919c 100644 --- a/extras/Hadrons/Modules/MContraction/Baryon.hpp +++ b/extras/Hadrons/Modules/MContraction/Baryon.hpp @@ -112,7 +112,7 @@ void TBaryon::execute(void) << " quarks '" << par().q1 << "', '" << par().q2 << "', and '" << par().q3 << "'" << std::endl; - XmlWriter writer(par().output); + CorrWriter writer(par().output); PropagatorField1 &q1 = *env().template getObject(par().q1); PropagatorField2 &q2 = *env().template getObject(par().q2); PropagatorField3 &q3 = *env().template getObject(par().q2); diff --git a/extras/Hadrons/Modules/MContraction/DiscLoop.hpp b/extras/Hadrons/Modules/MContraction/DiscLoop.hpp new file mode 100644 index 00000000..4ad12e90 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/DiscLoop.hpp @@ -0,0 +1,144 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/DiscLoop.hpp + +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_DiscLoop_hpp_ +#define Hadrons_DiscLoop_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * DiscLoop * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +class DiscLoopPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(DiscLoopPar, + std::string, q_loop, + Gamma::Algebra, gamma, + std::string, output); +}; + +template +class TDiscLoop: public Module +{ + TYPE_ALIASES(FImpl,); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + Gamma::Algebra, gamma, + std::vector, corr); + }; +public: + // constructor + TDiscLoop(const std::string name); + // destructor + virtual ~TDiscLoop(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(DiscLoop, TDiscLoop, MContraction); + +/****************************************************************************** + * TDiscLoop implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TDiscLoop::TDiscLoop(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TDiscLoop::getInput(void) +{ + std::vector in = {par().q_loop}; + + return in; +} + +template +std::vector TDiscLoop::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TDiscLoop::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TDiscLoop::execute(void) +{ + LOG(Message) << "Computing disconnected loop contraction '" << getName() + << "' using '" << par().q_loop << "' with " << par().gamma + << " insertion." << std::endl; + + CorrWriter writer(par().output); + PropagatorField &q_loop = *env().template getObject(par().q_loop); + LatticeComplex c(env().getGrid()); + Gamma gamma(par().gamma); + std::vector buf; + Result result; + + c = trace(gamma*q_loop); + sliceSum(c, buf, Tp); + + result.gamma = par().gamma; + result.corr.resize(buf.size()); + for (unsigned int t = 0; t < buf.size(); ++t) + { + result.corr[t] = TensorRemove(buf[t]); + } + + write(writer, "disc", result); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_DiscLoop_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp new file mode 100644 index 00000000..e5e73fa6 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/Gamma3pt.hpp @@ -0,0 +1,170 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/Gamma3pt.hpp + +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_Gamma3pt_hpp_ +#define Hadrons_Gamma3pt_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + * 3pt contraction with gamma matrix insertion. + * + * Schematic: + * + * q2 q3 + * /----<------*------<----¬ + * / gamma \ + * / \ + * i * * f + * \ / + * \ / + * \----------->----------/ + * q1 + * + * trace(g5*q1*adj(q2)*g5*gamma*q3) + */ + +/****************************************************************************** + * Gamma3pt * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +class Gamma3ptPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Gamma3ptPar, + std::string, q1, + std::string, q2, + std::string, q3, + Gamma::Algebra, gamma, + std::string, output); +}; + +template +class TGamma3pt: public Module +{ + TYPE_ALIASES(FImpl1, 1); + TYPE_ALIASES(FImpl2, 2); + TYPE_ALIASES(FImpl3, 3); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + Gamma::Algebra, gamma, + std::vector, corr); + }; +public: + // constructor + TGamma3pt(const std::string name); + // destructor + virtual ~TGamma3pt(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(Gamma3pt, ARG(TGamma3pt), MContraction); + +/****************************************************************************** + * TGamma3pt implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TGamma3pt::TGamma3pt(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TGamma3pt::getInput(void) +{ + std::vector in = {par().q1, par().q2, par().q3}; + + return in; +} + +template +std::vector TGamma3pt::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TGamma3pt::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TGamma3pt::execute(void) +{ + LOG(Message) << "Computing 3pt contractions '" << getName() << "' using" + << " quarks '" << par().q1 << "', '" << par().q2 << "' and '" + << par().q3 << "', with " << par().gamma << " insertion." + << std::endl; + + CorrWriter writer(par().output); + PropagatorField1 &q1 = *env().template getObject(par().q1); + PropagatorField2 &q2 = *env().template getObject(par().q2); + 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); + sliceSum(c, buf, Tp); + + result.gamma = par().gamma; + result.corr.resize(buf.size()); + for (unsigned int t = 0; t < buf.size(); ++t) + { + result.corr[t] = TensorRemove(buf[t]); + } + + write(writer, "gamma3pt", result); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Gamma3pt_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/Meson.hpp b/extras/Hadrons/Modules/MContraction/Meson.hpp index c8ea3371..4cbe1ac4 100644 --- a/extras/Hadrons/Modules/MContraction/Meson.hpp +++ b/extras/Hadrons/Modules/MContraction/Meson.hpp @@ -6,8 +6,10 @@ Source file: extras/Hadrons/Modules/MContraction/Meson.hpp Copyright (C) 2015 Copyright (C) 2016 +Copyright (C) 2017 Author: Antonin Portelli + 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 @@ -36,20 +38,39 @@ 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 sink & source, pairs of gamma matrices + (space-separated strings) in angled brackets (i.e. ), + in a sequence (e.g. ""). + + Special values: "all" - perform all possible contractions. + - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0."), + given as multiples of (2*pi) / L. +*/ + /****************************************************************************** * TMeson * ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MContraction) +typedef std::pair GammaPair; + class MesonPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar, - std::string, q1, - std::string, q2, - std::string, output, - Gamma::Algebra, gammaSource, - Gamma::Algebra, gammaSink); + std::string, q1, + std::string, q2, + std::string, gammas, + std::string, mom, + std::string, output); }; template @@ -61,7 +82,10 @@ public: class Result: Serializable { public: - GRID_SERIALIZABLE_CLASS_MEMBERS(Result, std::vector, corr); + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + Gamma::Algebra, gamma_snk, + Gamma::Algebra, gamma_src, + std::vector, corr); }; public: // constructor @@ -71,6 +95,7 @@ public: // dependencies/products virtual std::vector getInput(void); virtual std::vector getOutput(void); + virtual void parseGammaString(std::vector &gammaList); // execution virtual void execute(void); }; @@ -103,6 +128,32 @@ std::vector TMeson::getOutput(void) return output; } +template +void TMeson::parseGammaString(std::vector &gammaList) +{ + // Determine gamma matrices to insert at source/sink. + if (par().gammas.compare("all") == 0) + { + // Do all contractions. + unsigned int n_gam = Ns * Ns; + gammaList.resize(n_gam*n_gam); + for (unsigned int i = 1; i < Gamma::nGamma; i += 2) + { + for (unsigned int j = 1; j < Gamma::nGamma; j += 2) + { + gammaList.push_back(std::make_pair((Gamma::Algebra)i, + (Gamma::Algebra)j)); + } + } + } + else + { + // Parse individual contractions from input string. + gammaList = strToVec(par().gammas); + } +} + + // execution /////////////////////////////////////////////////////////////////// template void TMeson::execute(void) @@ -111,21 +162,44 @@ void TMeson::execute(void) << " quarks '" << par().q1 << "' and '" << par().q2 << "'" << std::endl; - XmlWriter writer(par().output); - PropagatorField1 &q1 = *env().template getObject(par().q1); - PropagatorField2 &q2 = *env().template getObject(par().q2); - LatticeComplex c(env().getGrid()); - Gamma gSrc(par().gammaSource), gSnk(par().gammaSink); - Gamma g5(Gamma::Algebra::Gamma5); - std::vector buf; - Result result; - - c = trace(gSnk*q1*adj(gSrc)*g5*adj(q2)*g5); - sliceSum(c, buf, Tp); - result.corr.resize(buf.size()); - for (unsigned int t = 0; t < buf.size(); ++t) + CorrWriter writer(par().output); + PropagatorField1 &q1 = *env().template getObject(par().q1); + PropagatorField2 &q2 = *env().template getObject(par().q2); + LatticeComplex c(env().getGrid()); + Gamma g5(Gamma::Algebra::Gamma5); + std::vector gammaList; + std::vector buf; + std::vector result; + std::vector p; + + p = strToVec(par().mom); + LatticeComplex ph(env().getGrid()), coor(env().getGrid()); + Complex i(0.0,1.0); + ph = zero; + for(unsigned int mu = 0; mu < env().getNd(); mu++) { - result.corr[t] = TensorRemove(buf[t]); + LatticeCoordinate(coor, mu); + ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); + } + ph = exp((Real)(2*M_PI)*i*ph); + + parseGammaString(gammaList); + + result.resize(gammaList.size()); + for (unsigned int i = 0; i < result.size(); ++i) + { + Gamma gSnk(gammaList[i].first); + Gamma gSrc(gammaList[i].second); + c = trace((g5*gSnk)*q1*(adj(gSrc)*g5)*adj(q2))*ph; + sliceSum(c, buf, Tp); + + result[i].gamma_snk = gammaList[i].first; + result[i].gamma_src = gammaList[i].second; + result[i].corr.resize(buf.size()); + for (unsigned int t = 0; t < buf.size(); ++t) + { + result[i].corr[t] = TensorRemove(buf[t]); + } } write(writer, "meson", result); } diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp new file mode 100644 index 00000000..23482feb --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp @@ -0,0 +1,114 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp + +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_WeakHamiltonian_hpp_ +#define Hadrons_WeakHamiltonian_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * WeakHamiltonian * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +/******************************************************************************* + * Utilities for contractions involving the Weak Hamiltonian. + ******************************************************************************/ +//// Sum and store correlator. +#define MAKE_DIAG(exp, buf, res, n)\ +sliceSum(exp, buf, Tp);\ +res.name = (n);\ +res.corr.resize(buf.size());\ +for (unsigned int t = 0; t < buf.size(); ++t)\ +{\ + res.corr[t] = TensorRemove(buf[t]);\ +} + +//// Contraction of mu index: use 'mu' variable in exp. +#define SUM_MU(buf,exp)\ +buf = zero;\ +for (unsigned int mu = 0; mu < ndim; ++mu)\ +{\ + buf += exp;\ +} + +enum +{ + i_V = 0, + i_A = 1, + n_i = 2 +}; + +class WeakHamiltonianPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WeakHamiltonianPar, + std::string, q1, + std::string, q2, + std::string, q3, + std::string, q4, + std::string, output); +}; + +#define MAKE_WEAK_MODULE(modname)\ +class T##modname: public Module\ +{\ +public:\ + TYPE_ALIASES(FIMPL,)\ + class Result: Serializable\ + {\ + public:\ + GRID_SERIALIZABLE_CLASS_MEMBERS(Result,\ + std::string, name,\ + std::vector, corr);\ + };\ +public:\ + /* constructor */ \ + T##modname(const std::string name);\ + /* destructor */ \ + virtual ~T##modname(void) = default;\ + /* dependency relation */ \ + virtual std::vector getInput(void);\ + virtual std::vector getOutput(void);\ + /* setup */ \ + virtual void setup(void);\ + /* execution */ \ + virtual void execute(void);\ + std::vector VA_label = {"V", "A"};\ +};\ +MODULE_REGISTER_NS(modname, T##modname, MContraction); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WeakHamiltonian_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc new file mode 100644 index 00000000..a44c2534 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.cc @@ -0,0 +1,137 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MContraction; + +/* + * Weak Hamiltonian current-current contractions, Eye-type. + * + * These contractions are generated by the Q1 and Q2 operators in the physical + * basis (see e.g. Fig 3 of arXiv:1507.03094). + * + * Schematics: q4 | + * /-<-¬ | + * / \ | q2 q3 + * \ / | /----<------*------<----¬ + * q2 \ / q3 | / /-*-¬ \ + * /-----<-----* *-----<----¬ | / / \ \ + * i * H_W * f | i * \ / q4 * f + * \ / | \ \->-/ / + * \ / | \ / + * \---------->---------/ | \----------->----------/ + * q1 | q1 + * | + * Saucer (S) | Eye (E) + * + * 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]) + */ + +/****************************************************************************** + * TWeakHamiltonianEye implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TWeakHamiltonianEye::TWeakHamiltonianEye(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TWeakHamiltonianEye::getInput(void) +{ + std::vector in = {par().q1, par().q2, par().q3, par().q4}; + + return in; +} + +std::vector TWeakHamiltonianEye::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TWeakHamiltonianEye::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +void TWeakHamiltonianEye::execute(void) +{ + LOG(Message) << "Computing Weak Hamiltonian (Eye type) contractions '" + << getName() << "' using quarks '" << par().q1 << "', '" + << par().q2 << ", '" << par().q3 << "' and '" << par().q4 + << "'." << 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(); + + PropagatorField tmp1(env().getGrid()); + LatticeComplex tmp2(env().getGrid()); + std::vector S_body(ndim, tmp1); + std::vector S_loop(ndim, tmp1); + std::vector E_body(ndim, tmp2); + std::vector E_loop(ndim, tmp2); + + // 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_loop[mu] = MAKE_SE_LOOP(q4, GammaL(Gamma::gmu[mu])); + } + + // Perform S-type contractions. + SUM_MU(expbuf, trace(S_body[mu]*S_loop[mu])) + MAKE_DIAG(expbuf, corrbuf, result[S_diag], "HW_S") + + // Recycle sub-expressions for E-type contractions. + for (unsigned int mu = 0; mu < ndim; ++mu) + { + E_body[mu] = trace(S_body[mu]); + E_loop[mu] = trace(S_loop[mu]); + } + + // Perform E-type contractions. + SUM_MU(expbuf, E_body[mu]*E_loop[mu]) + MAKE_DIAG(expbuf, corrbuf, result[E_diag], "HW_E") + + write(writer, "HW_Eye", result); +} diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp new file mode 100644 index 00000000..2ee87895 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -0,0 +1,58 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp + +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_WeakHamiltonianEye_hpp_ +#define Hadrons_WeakHamiltonianEye_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * WeakHamiltonianEye * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +enum +{ + S_diag = 0, + E_diag = 1, + n_eye_diag = 2 +}; + +// Saucer and Eye subdiagram contractions. +#define MAKE_SE_BODY(Q_1, Q_2, Q_3, gamma) (Q_3*g5*Q_1*adj(Q_2)*g5*gamma) +#define MAKE_SE_LOOP(Q_loop, gamma) (Q_loop*gamma) + +MAKE_WEAK_MODULE(WeakHamiltonianEye) + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WeakHamiltonianEye_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc new file mode 100644 index 00000000..2c4df68a --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.cc @@ -0,0 +1,139 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MContraction; + +/* + * Weak Hamiltonian current-current contractions, Non-Eye-type. + * + * These contractions are generated by the Q1 and Q2 operators in the physical + * basis (see e.g. Fig 3 of arXiv:1507.03094). + * + * Schematic: + * q2 q3 | q2 q3 + * /--<--¬ /--<--¬ | /--<--¬ /--<--¬ + * / \ / \ | / \ / \ + * / \ / \ | / \ / \ + * / \ / \ | / \ / \ + * i * * H_W * f | i * * * H_W * f + * \ * | | \ / \ / + * \ / \ / | \ / \ / + * \ / \ / | \ / \ / + * \ / \ / | \-->--/ \-->--/ + * \-->--/ \-->--/ | q1 q4 + * q1 q4 | + * Connected (C) | Wing (W) + * + * C: trace(q1*adj(q2)*g5*gL[mu]*q3*adj(q4)*g5*gL[mu]) + * W: trace(q1*adj(q2)*g5*gL[mu])*trace(q3*adj(q4)*g5*gL[mu]) + * + */ + +/****************************************************************************** + * TWeakHamiltonianNonEye implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TWeakHamiltonianNonEye::TWeakHamiltonianNonEye(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TWeakHamiltonianNonEye::getInput(void) +{ + std::vector in = {par().q1, par().q2, par().q3, par().q4}; + + return in; +} + +std::vector TWeakHamiltonianNonEye::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TWeakHamiltonianNonEye::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +void TWeakHamiltonianNonEye::execute(void) +{ + LOG(Message) << "Computing Weak Hamiltonian (Non-Eye type) contractions '" + << getName() << "' using quarks '" << par().q1 << "', '" + << par().q2 << ", '" << par().q3 << "' and '" << par().q4 + << "'." << 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_noneye_diag); + unsigned int ndim = env().getNd(); + + PropagatorField tmp1(env().getGrid()); + LatticeComplex tmp2(env().getGrid()); + std::vector C_i_side_loop(ndim, tmp1); + std::vector C_f_side_loop(ndim, tmp1); + std::vector W_i_side_loop(ndim, tmp2); + std::vector W_f_side_loop(ndim, tmp2); + + // Setup for C-type contractions. + for (int mu = 0; mu < ndim; ++mu) + { + C_i_side_loop[mu] = MAKE_CW_SUBDIAG(q1, q2, GammaL(Gamma::gmu[mu])); + C_f_side_loop[mu] = MAKE_CW_SUBDIAG(q3, q4, GammaL(Gamma::gmu[mu])); + } + + // Perform C-type contractions. + SUM_MU(expbuf, trace(C_i_side_loop[mu]*C_f_side_loop[mu])) + MAKE_DIAG(expbuf, corrbuf, result[C_diag], "HW_C") + + // Recycle sub-expressions for W-type contractions. + for (unsigned int mu = 0; mu < ndim; ++mu) + { + W_i_side_loop[mu] = trace(C_i_side_loop[mu]); + W_f_side_loop[mu] = trace(C_f_side_loop[mu]); + } + + // Perform W-type contractions. + SUM_MU(expbuf, W_i_side_loop[mu]*W_f_side_loop[mu]) + MAKE_DIAG(expbuf, corrbuf, result[W_diag], "HW_W") + + write(writer, "HW_NonEye", result); +} diff --git a/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp new file mode 100644 index 00000000..69bb8005 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp @@ -0,0 +1,57 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp + +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_WeakHamiltonianNonEye_hpp_ +#define Hadrons_WeakHamiltonianNonEye_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * WeakHamiltonianNonEye * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +enum +{ + W_diag = 0, + C_diag = 1, + n_noneye_diag = 2 +}; + +// Wing and Connected subdiagram contractions +#define MAKE_CW_SUBDIAG(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma) + +MAKE_WEAK_MODULE(WeakHamiltonianNonEye) + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WeakHamiltonianNonEye_hpp_ diff --git a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc new file mode 100644 index 00000000..6685f292 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc @@ -0,0 +1,135 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.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 +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MContraction; + +/* + * Weak Hamiltonian + current contractions, disconnected topology for neutral + * mesons. + * + * These contractions are generated by operators Q_1,...,10 of the dS=1 Weak + * Hamiltonian in the physical basis and an additional current J (see e.g. + * Fig 11 of arXiv:1507.03094). + * + * Schematic: + * + * q2 q4 q3 + * /--<--¬ /---<--¬ /---<--¬ + * / \ / \ / \ + * i * * H_W | J * * f + * \ / \ / \ / + * \--->---/ \-------/ \------/ + * q1 + * + * options + * - q1: input propagator 1 (string) + * - q2: input propagator 2 (string) + * - q3: input propagator 3 (string), assumed to be sequential propagator + * - q4: input propagator 4 (string), assumed to be a loop + * + * type 1: trace(q1*adj(q2)*g5*gL[mu])*trace(loop*gL[mu])*trace(q3*g5) + * type 2: trace(q1*adj(q2)*g5*gL[mu]*loop*gL[mu])*trace(q3*g5) + */ + +/******************************************************************************* + * TWeakNeutral4ptDisc implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +TWeakNeutral4ptDisc::TWeakNeutral4ptDisc(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +std::vector TWeakNeutral4ptDisc::getInput(void) +{ + std::vector in = {par().q1, par().q2, par().q3, par().q4}; + + return in; +} + +std::vector TWeakNeutral4ptDisc::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +void TWeakNeutral4ptDisc::setup(void) +{ + +} + +// execution /////////////////////////////////////////////////////////////////// +void TWeakNeutral4ptDisc::execute(void) +{ + LOG(Message) << "Computing Weak Hamiltonian neutral disconnected contractions '" + << getName() << "' using quarks '" << par().q1 << "', '" + << par().q2 << ", '" << par().q3 << "' and '" << par().q4 + << "'." << 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_neut_disc_diag); + unsigned int ndim = env().getNd(); + + PropagatorField tmp(env().getGrid()); + std::vector meson(ndim, tmp); + std::vector loop(ndim, tmp); + LatticeComplex curr(env().getGrid()); + + // Setup for type 1 contractions. + for (int mu = 0; mu < ndim; ++mu) + { + meson[mu] = MAKE_DISC_MESON(q1, q2, GammaL(Gamma::gmu[mu])); + loop[mu] = MAKE_DISC_LOOP(q4, GammaL(Gamma::gmu[mu])); + } + curr = MAKE_DISC_CURR(q3, GammaL(Gamma::Algebra::Gamma5)); + + // Perform type 1 contractions. + SUM_MU(expbuf, trace(meson[mu]*loop[mu])) + expbuf *= curr; + MAKE_DIAG(expbuf, corrbuf, result[neut_disc_1_diag], "HW_disc0_1") + + // Perform type 2 contractions. + SUM_MU(expbuf, trace(meson[mu])*trace(loop[mu])) + expbuf *= curr; + MAKE_DIAG(expbuf, corrbuf, result[neut_disc_2_diag], "HW_disc0_2") + + write(writer, "HW_disc0", result); +} diff --git a/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp new file mode 100644 index 00000000..c0d8f829 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp @@ -0,0 +1,59 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp + +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_WeakNeutral4ptDisc_hpp_ +#define Hadrons_WeakNeutral4ptDisc_hpp_ + +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * WeakNeutral4ptDisc * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +enum +{ + neut_disc_1_diag = 0, + neut_disc_2_diag = 1, + n_neut_disc_diag = 2 +}; + +// Neutral 4pt disconnected subdiagram contractions. +#define MAKE_DISC_MESON(Q_1, Q_2, gamma) (Q_1*adj(Q_2)*g5*gamma) +#define MAKE_DISC_LOOP(Q_LOOP, gamma) (Q_LOOP*gamma) +#define MAKE_DISC_CURR(Q_c, gamma) (trace(Q_c*gamma)) + +MAKE_WEAK_MODULE(WeakNeutral4ptDisc) + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WeakNeutral4ptDisc_hpp_ diff --git a/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp new file mode 100644 index 00000000..3d2850d1 --- /dev/null +++ b/extras/Hadrons/Modules/MLoop/NoiseLoop.hpp @@ -0,0 +1,132 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MLoop/NoiseLoop.hpp + +Copyright (C) 2016 + +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_NoiseLoop_hpp_ +#define Hadrons_NoiseLoop_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Noise loop propagator + ----------------------------- + * loop_x = q_x * adj(eta_x) + + * options: + - q = Result of inversion on noise source. + - eta = noise source. + + */ + + +/****************************************************************************** + * NoiseLoop * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MLoop) + +class NoiseLoopPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(NoiseLoopPar, + std::string, q, + std::string, eta); +}; + +template +class TNoiseLoop: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TNoiseLoop(const std::string name); + // destructor + virtual ~TNoiseLoop(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(NoiseLoop, TNoiseLoop, MLoop); + +/****************************************************************************** + * TNoiseLoop implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TNoiseLoop::TNoiseLoop(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TNoiseLoop::getInput(void) +{ + std::vector in = {par().q, par().eta}; + + return in; +} + +template +std::vector TNoiseLoop::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TNoiseLoop::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TNoiseLoop::execute(void) +{ + PropagatorField &loop = *env().template createLattice(getName()); + PropagatorField &q = *env().template getObject(par().q); + PropagatorField &eta = *env().template getObject(par().eta); + loop = q*adj(eta); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_NoiseLoop_hpp_ diff --git a/extras/Hadrons/Modules/MSource/SeqGamma.hpp b/extras/Hadrons/Modules/MSource/SeqGamma.hpp index 5cb8483f..366ebee7 100644 --- a/extras/Hadrons/Modules/MSource/SeqGamma.hpp +++ b/extras/Hadrons/Modules/MSource/SeqGamma.hpp @@ -6,6 +6,7 @@ Source file: extras/Hadrons/Modules/MSource/SeqGamma.hpp Copyright (C) 2015 Copyright (C) 2016 +Copyright (C) 2017 Author: Antonin Portelli @@ -149,9 +150,9 @@ void TSeqGamma::execute(void) for(unsigned int mu = 0; mu < env().getNd(); mu++) { LatticeCoordinate(coor, mu); - ph = ph + p[mu]*coor; + ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); } - ph = exp(i*ph); + ph = exp((Real)(2*M_PI)*i*ph); LatticeCoordinate(t, Tp); src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q); } diff --git a/extras/Hadrons/Modules/MSource/Wall.hpp b/extras/Hadrons/Modules/MSource/Wall.hpp new file mode 100644 index 00000000..8722876f --- /dev/null +++ b/extras/Hadrons/Modules/MSource/Wall.hpp @@ -0,0 +1,147 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MSource/Wall.hpp + +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 +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_WallSource_hpp_ +#define Hadrons_WallSource_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* + + Wall source + ----------------------------- + * src_x = delta(x_3 - tW) * exp(i x.mom) + + * options: + - tW: source timeslice (integer) + - mom: momentum insertion, space-separated float sequence (e.g ".1 .2 1. 0.") + + */ + +/****************************************************************************** + * Wall * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MSource) + +class WallPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(WallPar, + unsigned int, tW, + std::string, mom); +}; + +template +class TWall: public Module +{ +public: + TYPE_ALIASES(FImpl,); +public: + // constructor + TWall(const std::string name); + // destructor + virtual ~TWall(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(Wall, TWall, MSource); + +/****************************************************************************** + * TWall implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TWall::TWall(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TWall::getInput(void) +{ + std::vector in; + + return in; +} + +template +std::vector TWall::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TWall::setup(void) +{ + env().template registerLattice(getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TWall::execute(void) +{ + LOG(Message) << "Generating wall source at t = " << par().tW + << " with momentum " << par().mom << std::endl; + + PropagatorField &src = *env().template createLattice(getName()); + Lattice> t(env().getGrid()); + LatticeComplex ph(env().getGrid()), coor(env().getGrid()); + std::vector p; + Complex i(0.0,1.0); + + p = strToVec(par().mom); + ph = zero; + for(unsigned int mu = 0; mu < Nd; mu++) + { + LatticeCoordinate(coor, mu); + ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu]))); + } + ph = exp((Real)(2*M_PI)*i*ph); + LatticeCoordinate(t, Tp); + src = 1.; + src = where((t == par().tW), src*ph, 0.*src); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WallSource_hpp_ diff --git a/extras/Hadrons/Modules/Quark.hpp b/extras/Hadrons/Modules/Quark.hpp index e441a096..be7426ab 100644 --- a/extras/Hadrons/Modules/Quark.hpp +++ b/extras/Hadrons/Modules/Quark.hpp @@ -173,7 +173,7 @@ void TQuark::execute(void) *env().template getObject(getName()); axpby_ssp_pminus(sol, 0., sol, 1., sol, 0, 0); - axpby_ssp_pplus(sol, 0., sol, 1., sol, 0, Ls_-1); + axpby_ssp_pplus(sol, 1., sol, 1., sol, 0, Ls_-1); ExtractSlice(tmp, sol, 0, 0); FermToProp(p4d, tmp, s, c); } diff --git a/extras/Hadrons/modules.inc b/extras/Hadrons/modules.inc index 4251ffa3..af291631 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,4 +1,7 @@ modules_cc =\ + Modules/MContraction/WeakHamiltonianEye.cc \ + Modules/MContraction/WeakHamiltonianNonEye.cc \ + Modules/MContraction/WeakNeutral4ptDisc.cc \ Modules/MGauge/Load.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/Unit.cc @@ -7,13 +10,21 @@ modules_hpp =\ Modules/MAction/DWF.hpp \ Modules/MAction/Wilson.hpp \ Modules/MContraction/Baryon.hpp \ + Modules/MContraction/DiscLoop.hpp \ + Modules/MContraction/Gamma3pt.hpp \ Modules/MContraction/Meson.hpp \ + Modules/MContraction/WeakHamiltonian.hpp \ + Modules/MContraction/WeakHamiltonianEye.hpp \ + Modules/MContraction/WeakHamiltonianNonEye.hpp \ + Modules/MContraction/WeakNeutral4ptDisc.hpp \ Modules/MGauge/Load.hpp \ Modules/MGauge/Random.hpp \ Modules/MGauge/Unit.hpp \ + Modules/MLoop/NoiseLoop.hpp \ Modules/MSolver/RBPrecCG.hpp \ Modules/MSource/Point.hpp \ Modules/MSource/SeqGamma.hpp \ + Modules/MSource/Wall.hpp \ Modules/MSource/Z2.hpp \ Modules/Quark.hpp diff --git a/grid-config.in b/grid-config.in new file mode 100755 index 00000000..bd340846 --- /dev/null +++ b/grid-config.in @@ -0,0 +1,86 @@ +#! /bin/sh + +prefix=@prefix@ +exec_prefix=@exec_prefix@ +includedir=@includedir@ + +usage() +{ + cat < #ifndef GRID_BASE_H #define GRID_BASE_H -/////////////////// -// Std C++ dependencies -/////////////////// -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -/////////////////// -// Grid headers -/////////////////// -#include "Config.h" +#include #include #include diff --git a/lib/GridStd.h b/lib/GridStd.h new file mode 100644 index 00000000..fb5e5b21 --- /dev/null +++ b/lib/GridStd.h @@ -0,0 +1,27 @@ +#ifndef GRID_STD_H +#define GRID_STD_H + +/////////////////// +// Std C++ dependencies +/////////////////// +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/////////////////// +// Grid config +/////////////////// +#include "Config.h" + +#endif /* GRID_STD_H */ diff --git a/lib/algorithms/approx/Remez.h b/lib/algorithms/approx/Remez.h index 31938779..71b1093b 100644 --- a/lib/algorithms/approx/Remez.h +++ b/lib/algorithms/approx/Remez.h @@ -16,7 +16,7 @@ #define INCLUDED_ALG_REMEZ_H #include -#include +#include #ifdef HAVE_LIBGMP #include "bigfloat.h" diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index 51677522..ff814eb6 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -320,7 +320,7 @@ void CayleyFermion5D::MDeriv (GaugeField &mat,const FermionField &U,const this->DhopDeriv(mat,U,Din,dag); } else { // U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call - MeooeDag5D(U,Din); + Meooe5D(U,Din); this->DhopDeriv(mat,Din,V,dag); } }; @@ -335,8 +335,8 @@ void CayleyFermion5D::MoeDeriv(GaugeField &mat,const FermionField &U,const this->DhopDerivOE(mat,U,Din,dag); } else { // U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call - MeooeDag5D(U,Din); - this->DhopDerivOE(mat,Din,V,dag); + Meooe5D(U,Din); + this->DhopDerivOE(mat,Din,V,dag); } }; template @@ -350,7 +350,7 @@ void CayleyFermion5D::MeoDeriv(GaugeField &mat,const FermionField &U,const this->DhopDerivEO(mat,U,Din,dag); } else { // U d/du [D_w D5]^dag V = U D5^dag d/du DW^dag Y // implicit adj on U in call - MeooeDag5D(U,Din); + Meooe5D(U,Din); this->DhopDerivEO(mat,Din,V,dag); } }; @@ -380,6 +380,8 @@ void CayleyFermion5D::SetCoefficientsInternal(RealD zolo_hi,std::vector::SetCoefficientsInternal(RealD zolo_hi,std::vector0.0); bs[i] = 0.5*(bpc/omega[i] + bmc); cs[i] = 0.5*(bpc/omega[i] - bmc); } - + //////////////////////////////////////////////////////// // Constants for the preconditioned matrix Cayley form //////////////////////////////////////////////////////// @@ -425,12 +428,12 @@ void CayleyFermion5D::SetCoefficientsInternal(RealD zolo_hi,std::vectorM5) +1.0); + bee[i]=as[i]*(bs[i]*(4.0-this->M5) +1.0); + // assert(fabs(bee[i])>0.0); cee[i]=as[i]*(1.0-cs[i]*(4.0-this->M5)); beo[i]=as[i]*bs[i]; ceo[i]=-as[i]*cs[i]; } - aee.resize(Ls); aeo.resize(Ls); for(int i=0;i::SetCoefficientsInternal(RealD zolo_hi,std::vector0.0); + delta_d *= cee[j]/bee[j]; + } dee[Ls-1] += delta_d; } int inv=1; this->MooeeInternalCompute(0,inv,MatpInv,MatmInv); this->MooeeInternalCompute(1,inv,MatpInvDag,MatmInvDag); - } @@ -495,7 +500,9 @@ void CayleyFermion5D::MooeeInternalCompute(int dag, int inv, GridBase *grid = this->FermionRedBlackGrid(); int LLs = grid->_rdimensions[0]; - if ( LLs == Ls ) return; // Not vectorised in 5th direction + if ( LLs == Ls ) { + return; // Not vectorised in 5th direction + } Eigen::MatrixXcd Pplus = Eigen::MatrixXcd::Zero(Ls,Ls); Eigen::MatrixXcd Pminus = Eigen::MatrixXcd::Zero(Ls,Ls); diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h index ad4bf87e..72ce8f42 100644 --- a/lib/qcd/action/fermion/DomainWallFermion.h +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -68,7 +68,7 @@ namespace Grid { Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham assert(zdata->n==this->Ls); - // std::cout<" + template + inline std::istream & operator>>(std::istream &is, std::pair &buf) + { + T1 buf1; + T2 buf2; + char c; + + // Search for "pair" delimiters. + do + { + is.get(c); + } while (c != '<' && !is.eof()); + if (c == '<') + { + int start = is.tellg(); + do + { + is.get(c); + } while (c != '>' && !is.eof()); + if (c == '>') + { + int end = is.tellg(); + int psize = end - start - 1; + + // Only read data between pair limiters. + is.seekg(start); + std::string tmpstr(psize, ' '); + is.read(&tmpstr[0], psize); + std::istringstream temp(tmpstr); + temp >> buf1 >> buf2; + buf = std::make_pair(buf1, buf2); + is.seekg(end); + } + } + is.peek(); + return is; + } + + // output to streams for pairs + template + inline std::ostream & operator<<(std::ostream &os, const std::pair &p) + { + os << "<" << p.first << " " << p.second << ">"; + return os; + } + // Abstract writer/reader classes //////////////////////////////////////////// // static polymorphism implemented using CRTP idiom class Serializable; diff --git a/scripts/configure.commit b/scripts/configure.commit index afc662eb..4fa6e822 100755 --- a/scripts/configure.commit +++ b/scripts/configure.commit @@ -16,4 +16,4 @@ else # configure.ac is modified printf 'commit: %s-dirty\n' `git rev-parse --short HEAD` printf 'branch: %s\n' `git rev-parse --abbrev-ref HEAD` printf 'date : %s\n' `git log -1 --date=short --pretty=format:%cd` - fi \ No newline at end of file +fi diff --git a/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 212b653f..384a001f 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -115,6 +115,7 @@ int main(int argc,char **argv) // test serializable class writing myclass obj(1234); // non-trivial constructor std::vector vec; + std::pair pair; std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl; write(WR,"obj",obj); @@ -122,6 +123,8 @@ int main(int argc,char **argv) vec.push_back(myclass(1234)); vec.push_back(myclass(5678)); vec.push_back(myclass(3838)); + pair = std::make_pair(myenum::red, myenum::blue); + write(WR, "objvec", vec); std::cout << "-- serialisable class writing to std::cout:" << std::endl; std::cout << obj << std::endl; @@ -129,21 +132,30 @@ int main(int argc,char **argv) std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl; std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl; + write(WR, "objpair", pair); + std::cout << "-- pair writing to std::cout:" << std::endl; + std::cout << pair << std::endl; + // read tests std::cout << "\n==== IO self-consistency tests" << std::endl; //// XML ioTest("iotest.xml", obj, "XML (object) "); ioTest("iotest.xml", vec, "XML (vector of objects)"); + ioTest("iotest.xml", pair, "XML (pair of objects)"); //// binary ioTest("iotest.bin", obj, "binary (object) "); ioTest("iotest.bin", vec, "binary (vector of objects)"); + ioTest("iotest.bin", pair, "binary (pair of objects)"); //// text ioTest("iotest.dat", obj, "text (object) "); ioTest("iotest.dat", vec, "text (vector of objects)"); + ioTest("iotest.dat", pair, "text (pair of objects)"); //// HDF5 +#undef HAVE_HDF5 #ifdef HAVE_HDF5 ioTest("iotest.h5", obj, "HDF5 (object) "); ioTest("iotest.h5", vec, "HDF5 (vector of objects)"); + ioTest("iotest.h5", pair, "HDF5 (pair of objects)"); #endif std::cout << "\n==== vector flattening/reconstruction" << std::endl; diff --git a/tests/core/Test_mobius_even_odd.cc b/tests/core/Test_mobius_even_odd.cc new file mode 100644 index 00000000..3c9ac520 --- /dev/null +++ b/tests/core/Test_mobius_even_odd.cc @@ -0,0 +1,287 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_even_odd.cc + + Copyright (C) 2015 + +Author: Peter Boyle +Author: paboyle + + 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 + *************************************************************************************/ + /* END LEGAL */ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +template +struct scal { + d internal; +}; + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout< seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + + LatticeFermion src (FGrid); random(RNG5,src); + LatticeFermion phi (FGrid); random(RNG5,phi); + LatticeFermion chi (FGrid); random(RNG5,chi); + LatticeFermion result(FGrid); result=zero; + LatticeFermion ref(FGrid); ref=zero; + LatticeFermion tmp(FGrid); tmp=zero; + LatticeFermion err(FGrid); tmp=zero; + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + std::vector U(4,UGrid); + + // Only one non-zero (y) + Umu=zero; + for(int nn=0;nn0 ) + U[nn]=zero; + PokeIndex(Umu,U[nn],nn); + } + + RealD mass=0.1; + RealD M5 =1.8; + std::vector < std::complex > omegas; +#if 0 + for(int i=0;i temp (0.25+0.01*i, imag*0.01); + omegas.push_back(temp); + } +#else + omegas.push_back( std::complex(1.45806438985048,-0) ); + omegas.push_back( std::complex(1.18231318389348,-0) ); + omegas.push_back( std::complex(0.830951166685955,-0) ); + omegas.push_back( std::complex(0.542352409156791,-0) ); + omegas.push_back( std::complex(0.341985020453729,-0) ); + omegas.push_back( std::complex(0.21137902619029,-0) ); + omegas.push_back( std::complex(0.126074299502912,-0) ); + omegas.push_back( std::complex(0.0990136651962626,-0) ); + omegas.push_back( std::complex(0.0686324988446592,0.0550658530827402) ); + omegas.push_back( std::complex(0.0686324988446592,-0.0550658530827402) ); +#endif + + MobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, 0.5,0.5); + // DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + + LatticeFermion src_e (FrbGrid); + LatticeFermion src_o (FrbGrid); + LatticeFermion r_e (FrbGrid); + LatticeFermion r_o (FrbGrid); + LatticeFermion r_eo (FGrid); + LatticeFermion r_eeoo(FGrid); + + std::cout< * = < chi | Deo^dag| phi> "< * = < chi | Deo^dag| phi> "< HermOpEO(Ddwf); + HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); + HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); + + HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2); + HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2); + + pDce = innerProduct(phi_e,dchi_e); + pDco = innerProduct(phi_o,dchi_o); + cDpe = innerProduct(chi_e,dphi_e); + cDpo = innerProduct(chi_o,dphi_o); + + std::cout< dist(1, 100); - - auto gen = std::bind(dist, mersenne_engine); - std::vector seeds4(4); - generate(begin(seeds4), end(seeds4), gen); - - //std::vector seeds4({1,2,3,5}); - std::vector seeds5({5,6,7,1}); + std::vector seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); @@ -76,15 +65,13 @@ int main (int argc, char ** argv) //////////////////////////////////// // Unmodified matrix element //////////////////////////////////// - RealD mass = 0.01; - RealD M5 = 1.8; - RealD b = 2.0; - RealD c = 1.0; + RealD mass=0.01; + RealD M5=1.8; + RealD b=0.5; + RealD c=0.5; MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); Ddwf.M (phi,Mphi); - - ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p // get the deriv of phidag MdagM phi with respect to "U" diff --git a/tests/forces/Test_wilson_force_phiMdagMphi.cc b/tests/forces/Test_wilson_force_phiMdagMphi.cc deleted file mode 100644 index 2a5814fe..00000000 --- a/tests/forces/Test_wilson_force_phiMdagMphi.cc +++ /dev/null @@ -1,167 +0,0 @@ - /************************************************************************************* - - Grid physics library, www.github.com/paboyle/Grid - - Source file: ./tests/Test_wilson_force_phiMdagMphi.cc - - Copyright (C) 2015 - -Author: Peter Boyle - - 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 - *************************************************************************************/ - /* END LEGAL */ -#include - -using namespace std; -using namespace Grid; -using namespace Grid::QCD; - - - -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); - - std::vector latt_size = GridDefaultLatt(); - std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); - - GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); - - int threads = GridThread::GetThreads(); - std::cout< seeds({1,2,3,4}); - - GridParallelRNG pRNG(&Grid); - pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); - - LatticeFermion phi (&Grid); gaussian(pRNG,phi); - LatticeFermion Mphi (&Grid); - LatticeFermion Mdagphi (&Grid); - LatticeFermion MphiPrime (&Grid); - LatticeFermion MdagphiPrime (&Grid); - LatticeFermion dMphi (&Grid); - - LatticeGaugeField U(&Grid); - - - SU3::HotConfiguration(pRNG,U); - // SU3::ColdConfiguration(pRNG,U); - - //////////////////////////////////// - // Unmodified matrix element - //////////////////////////////////// - RealD mass=-4.0; //kills the diagonal term - WilsonFermionR Dw (U, Grid,RBGrid,mass); - Dw.M (phi,Mphi); - Dw.Mdag(phi,Mdagphi); - - ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p - ComplexD Sdag = innerProduct(Mdagphi,Mdagphi); // pdag MMdag p - - // get the deriv of phidag MdagM phi with respect to "U" - LatticeGaugeField UdSdU(&Grid); - LatticeGaugeField UdSdUdag(&Grid); - LatticeGaugeField tmp(&Grid); - - Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp; - - Dw.MDeriv(tmp , Mdagphi, phi,DaggerYes ); UdSdUdag=tmp; - - - LatticeFermion dMdagphi (&Grid); dMdagphi=zero; - LatticeFermion Ftmp (&Grid); - - - // Dw.MDeriv(UdSdU,Mdagphi, phi,DaggerYes );// UdSdU =UdSdU +tmp; - - //////////////////////////////////// - // Modify the gauge field a little in one dir - //////////////////////////////////// - RealD dt = 1.0e-3; - - LatticeColourMatrix mommu(&Grid); - LatticeGaugeField mom(&Grid); - LatticeGaugeField Uprime(&Grid); - - for(int mu=0;mu(mom,mommu,mu); - - parallel_for(auto i=mom.begin();i - - 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 - *************************************************************************************/ - /* END LEGAL */ -#include - -using namespace std; -using namespace Grid; -using namespace Grid::QCD; - - - -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); - - std::vector latt_size = GridDefaultLatt(); - std::vector simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); - - GridCartesian Grid(latt_size,simd_layout,mpi_layout); - GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout); - - int threads = GridThread::GetThreads(); - std::cout< seeds({1,2,3,4}); - - GridParallelRNG pRNG(&Grid); - pRNG.SeedFixedIntegers(std::vector({45,12,81,9})); - - LatticeFermion phi (&Grid); gaussian(pRNG,phi); - LatticeFermion Mphi (&Grid); - LatticeFermion MphiPrime (&Grid); - LatticeFermion dMphi (&Grid); - - LatticeGaugeField U(&Grid); - - - SU3::HotConfiguration(pRNG,U); - - //////////////////////////////////// - // Unmodified matrix element - //////////////////////////////////// - RealD mass=-4.0; //kills the diagonal term - WilsonFermionR Dw (U, Grid,RBGrid,mass); - Dw.M(phi,Mphi); - - ComplexD S = innerProduct(phi,Mphi); - - // get the deriv - LatticeGaugeField UdSdU(&Grid); - Dw.MDeriv(UdSdU,phi, phi,DaggerNo ); - - - //////////////////////////////////// - // Modify the gauge field a little in one dir - //////////////////////////////////// - RealD dt = 1.0e-3; - Complex Complex_i(0,1); - - LatticeColourMatrix Umu(&Grid); - LatticeColourMatrix Umu_save(&Grid); - LatticeColourMatrix dU (&Grid); - LatticeColourMatrix mom(&Grid); - SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mom); // Traceless antihermitian momentum; gaussian in lie alg - - - // check mom is as i expect - LatticeColourMatrix tmpmom(&Grid); - tmpmom = mom+adj(mom); - std::cout << GridLogMessage << "mom anti-herm check "<< norm2(tmpmom)<(U,mu); - Umu_save=Umu; - dU = mom * Umu * dt; - Umu= Umu+dU; - PokeIndex(Dw.Umu,Umu,mu); - - Dw.M(phi,MphiPrime); - - ComplexD Sprime = innerProduct(phi,MphiPrime); - - std::cout << GridLogMessage << " S "<(Dw.Umu,dU,mu); - Dw.M(phi,dMphi); - - - ComplexD deltaS = innerProduct(phi,dMphi); - std::cout << GridLogMessage << "deltaS "<(Dw.Umu,Umu_save,mu); - Dw.Mdir(phi,dMphi,mu,1); - dMphi = dt*mom*dMphi; - - deltaS = innerProduct(phi,dMphi); - std::cout << GridLogMessage << "deltaS from inner prod of mom* M[u] "< + + 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; + +/******************************************************************************* + * Macros to reduce code duplication. + ******************************************************************************/ +// Useful definitions +#define ZERO_MOM "0. 0. 0. 0." +#define INIT_INDEX(s, n) (std::string(s) + "_" + std::to_string(n)) +#define ADD_INDEX(s, n) (s + "_" + std::to_string(n)) +#define LABEL_3PT(s, t1, t2) ADD_INDEX(INIT_INDEX(s, t1), t2) +#define LABEL_4PT(s, t1, t2, t3) ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3) +#define LABEL_4PT_NOISE(s, t1, t2, t3, nn) ADD_INDEX(ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3), nn) + +// Wall source/sink macros +#define NAME_3MOM_WALL_SOURCE(t, mom) ("wall_" + std::to_string(t) + "_" + mom) +#define NAME_WALL_SOURCE(t) NAME_3MOM_WALL_SOURCE(t, ZERO_MOM) +#define NAME_POINT_SOURCE(pos) ("point_" + pos) + +#define MAKE_3MOM_WALL_PROP(tW, mom, propName, solver)\ +{\ + std::string srcName = NAME_3MOM_WALL_SOURCE(tW, mom);\ + makeWallSource(application, srcName, tW, mom);\ + makePropagator(application, propName, srcName, solver);\ +} + +#define MAKE_WALL_PROP(tW, propName, solver)\ + MAKE_3MOM_WALL_PROP(tW, ZERO_MOM, propName, solver) + +// Sequential source macros +#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, propName, solver)\ +{\ + std::string srcName = ADD_INDEX(qSrc + "_seq", tS);\ + makeSequentialSource(application, srcName, qSrc, tS, mom);\ + makePropagator(application, propName, srcName, solver);\ +} + +// Point source macros +#define MAKE_POINT_PROP(pos, propName, solver)\ +{\ + std::string srcName = NAME_POINT_SOURCE(pos);\ + makePointSource(application, srcName, pos);\ + makePropagator(application, propName, srcName, solver);\ +} + +/******************************************************************************* + * Functions for propagator construction. + ******************************************************************************/ + +/******************************************************************************* + * Name: makePointSource + * Purpose: Construct point source and add to application module. + * Parameters: application - main application that stores modules. + * srcName - name of source module to create. + * pos - Position of point source. + * Returns: None. + ******************************************************************************/ +inline void makePointSource(Application &application, std::string srcName, + std::string pos) +{ + // If the source already exists, don't make the module again. + if (!(Environment::getInstance().hasModule(srcName))) + { + MSource::Point::Par pointPar; + pointPar.position = pos; + application.createModule(srcName, pointPar); + } +} + +/******************************************************************************* + * Name: makeSequentialSource + * Purpose: Construct sequential source and add to application module. + * Parameters: application - main application that stores modules. + * srcName - name of source module to create. + * qSrc - Input quark for sequential inversion. + * tS - sequential source timeslice. + * mom - momentum insertion (default is zero). + * Returns: None. + ******************************************************************************/ +inline void makeSequentialSource(Application &application, std::string srcName, + std::string qSrc, unsigned int tS, + std::string mom = ZERO_MOM) +{ + // If the source already exists, don't make the module again. + if (!(Environment::getInstance().hasModule(srcName))) + { + MSource::SeqGamma::Par seqPar; + seqPar.q = qSrc; + seqPar.tA = tS; + seqPar.tB = tS; + seqPar.mom = mom; + seqPar.gamma = Gamma::Algebra::GammaT; + application.createModule(srcName, seqPar); + } +} + +/******************************************************************************* + * Name: makeWallSource + * Purpose: Construct wall source and add to application module. + * Parameters: application - main application that stores modules. + * srcName - name of source module to create. + * tW - wall source timeslice. + * mom - momentum insertion (default is zero). + * Returns: None. + ******************************************************************************/ +inline void makeWallSource(Application &application, std::string srcName, + unsigned int tW, std::string mom = ZERO_MOM) +{ + // If the source already exists, don't make the module again. + if (!(Environment::getInstance().hasModule(srcName))) + { + MSource::Wall::Par wallPar; + wallPar.tW = tW; + wallPar.mom = mom; + application.createModule(srcName, wallPar); + } +} + +/******************************************************************************* + * Name: makeWallSink + * Purpose: Wall sink smearing of a propagator. + * Parameters: application - main application that stores modules. + * propName - name of input propagator. + * wallName - name of smeared propagator. + * mom - momentum insertion (default is zero). + * Returns: None. + ******************************************************************************/ +inline void makeWallSink(Application &application, std::string propName, + std::string wallName, std::string mom = ZERO_MOM) +{ + // 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))) + { + MSink::Wall::Par wallPar; + wallPar.q = propName; + wallPar.mom = mom; + application.createModule(wallName, wallPar); + }*/ +} + +/******************************************************************************* + * Name: makePropagator + * Purpose: Construct source and propagator then add to application module. + * Parameters: application - main application that stores modules. + * propName - name of propagator module to create. + * srcName - name of source module to use. + * solver - solver to use (default is CG). + * Returns: None. + ******************************************************************************/ +inline void makePropagator(Application &application, std::string &propName, + std::string &srcName, std::string &solver) +{ + // If the propagator already exists, don't make the module again. + if (!(Environment::getInstance().hasModule(propName))) + { + Quark::Par quarkPar; + quarkPar.source = srcName; + quarkPar.solver = solver; + application.createModule(propName, quarkPar); + } +} + +/******************************************************************************* + * Name: makeLoop + * Purpose: Use noise source and inversion result to make loop propagator, then + * add to application module. + * Parameters: application - main application that stores modules. + * propName - name of propagator module to create. + * srcName - name of noise source module to use. + * resName - name of inversion result on given noise source. + * Returns: None. + ******************************************************************************/ +inline void makeLoop(Application &application, std::string &propName, + std::string &srcName, std::string &resName) +{ + // If the loop propagator already exists, don't make the module again. + if (!(Environment::getInstance().hasModule(propName))) + { + MLoop::NoiseLoop::Par loopPar; + loopPar.q = resName; + loopPar.eta = srcName; + application.createModule(propName, loopPar); + } +} + +/******************************************************************************* + * Contraction module creation. + ******************************************************************************/ + +/******************************************************************************* + * Name: mesonContraction + * Purpose: Create meson 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. + * q2 - quark propagator 2. + * label - unique label to construct module name. + * mom - momentum to project (default is zero) + * gammas - gamma insertions at source and sink. + * Returns: None. + ******************************************************************************/ +inline void mesonContraction(Application &application, unsigned int npt, + std::string &q1, std::string &q2, + std::string &label, + std::string mom = ZERO_MOM, + std::string gammas = "") +{ + std::string modName = std::to_string(npt) + "pt_" + label; + if (!(Environment::getInstance().hasModule(modName))) + { + MContraction::Meson::Par mesPar; + mesPar.output = std::to_string(npt) + "pt/" + label; + mesPar.q1 = q1; + mesPar.q2 = q2; + mesPar.mom = mom; + mesPar.gammas = gammas; + application.createModule(modName, mesPar); + } + } + +/******************************************************************************* + * Name: gamma3ptContraction + * 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. + * q2 - quark propagator 2. + * q3 - quark propagator 3. + * label - unique label to construct module name. + * 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, + Gamma::Algebra gamma = Gamma::Algebra::Identity) +{ + std::string modName = std::to_string(npt) + "pt_" + label; + if (!(Environment::getInstance().hasModule(modName))) + { + MContraction::Gamma3pt::Par gamma3ptPar; + gamma3ptPar.output = std::to_string(npt) + "pt/" + label; + gamma3ptPar.q1 = q1; + gamma3ptPar.q2 = q2; + gamma3ptPar.q3 = q3; + gamma3ptPar.gamma = gamma; + application.createModule(modName, gamma3ptPar); + } + } + +/******************************************************************************* + * Name: weakContraction[Eye,NonEye] + * Purpose: Create Weak Hamiltonian contraction module for Eye/NonEye topology + * and add to application module. + * Parameters: application - main application that stores modules. + * npt - specify n-point correlator (for labelling). + * q1 - quark propagator 1. + * q2 - quark propagator 2. + * q3 - quark propagator 3. + * q4 - quark propagator 4. + * label - unique label to construct module name. + * 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 modName = std::to_string(npt) + "pt_" + label;\ + if (!(Environment::getInstance().hasModule(modName)))\ + {\ + MContraction::WeakHamiltonian##top::Par weakPar;\ + weakPar.output = std::to_string(npt) + "pt/" + label;\ + weakPar.q1 = q1;\ + weakPar.q2 = q2;\ + weakPar.q3 = q3;\ + weakPar.q4 = q4;\ + application.createModule(modName, weakPar);\ + }\ +} +HW_CONTRACTION(Eye) // weakContractionEye +HW_CONTRACTION(NonEye) // weakContractionNonEye + +/******************************************************************************* + * Name: disc0Contraction + * Purpose: Create contraction module for 4pt Weak Hamiltonian + current + * disconnected topology for neutral mesons and add to application + * module. + * Parameters: application - main application that stores modules. + * q1 - quark propagator 1. + * q2 - quark propagator 2. + * q3 - quark propagator 3. + * q4 - quark propagator 4. + * label - unique label to construct module name. + * Returns: None. + ******************************************************************************/ +inline void disc0Contraction(Application &application, + std::string &q1, std::string &q2, + std::string &q3, std::string &q4, + std::string &label) +{ + std::string modName = "4pt_" + label; + if (!(Environment::getInstance().hasModule(modName))) + { + MContraction::WeakNeutral4ptDisc::Par disc0Par; + disc0Par.output = "4pt/" + label; + disc0Par.q1 = q1; + disc0Par.q2 = q2; + disc0Par.q3 = q3; + disc0Par.q4 = q4; + application.createModule(modName, disc0Par); + } + } + +/******************************************************************************* + * Name: discLoopContraction + * Purpose: Create contraction module for disconnected loop and add to + * application module. + * Parameters: application - main application that stores modules. + * q_loop - loop quark propagator. + * modName - unique module name. + * gamma - gamma matrix to use in contraction. + * Returns: None. + ******************************************************************************/ +inline void discLoopContraction(Application &application, + std::string &q_loop, std::string &modName, + Gamma::Algebra gamma = Gamma::Algebra::Identity) +{ + if (!(Environment::getInstance().hasModule(modName))) + { + MContraction::DiscLoop::Par discPar; + discPar.output = "disc/" + modName; + discPar.q_loop = q_loop; + discPar.gamma = gamma; + application.createModule(modName, discPar); + } + } diff --git a/tests/hadrons/Test_hadrons_meson_3pt.cc b/tests/hadrons/Test_hadrons_meson_3pt.cc index da7ddb84..efef6931 100644 --- a/tests/hadrons/Test_hadrons_meson_3pt.cc +++ b/tests/hadrons/Test_hadrons_meson_3pt.cc @@ -30,14 +30,6 @@ using namespace Grid; using namespace Hadrons; -static Gamma::Algebra gmu[4] = -{ - Gamma::Algebra::GammaX, - Gamma::Algebra::GammaY, - Gamma::Algebra::GammaZ, - Gamma::Algebra::GammaT -}; - int main(int argc, char *argv[]) { // initialization ////////////////////////////////////////////////////////// @@ -110,7 +102,7 @@ int main(int argc, char *argv[]) seqName.push_back(std::vector(Nd)); for (unsigned int mu = 0; mu < Nd; ++mu) { - seqPar.gamma = gmu[mu]; + seqPar.gamma = 0x1 << mu; seqName[i][mu] = "G" + std::to_string(seqPar.gamma) + "_" + std::to_string(seqPar.tA) + "-" + qName[i]; @@ -135,11 +127,11 @@ int main(int argc, char *argv[]) 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]; - mesPar.gammaSource = Gamma::Algebra::Gamma5; - mesPar.gammaSink = Gamma::Algebra::Gamma5; + mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j]; + mesPar.q1 = qName[i]; + mesPar.q2 = qName[j]; + mesPar.gammas = "all"; + mesPar.mom = "0. 0. 0. 0."; application.createModule("meson_Z2_" + std::to_string(t) + "_" @@ -157,6 +149,8 @@ int main(int argc, char *argv[]) + std::to_string(mu); mesPar.q1 = qName[i]; mesPar.q2 = seqName[j][mu]; + mesPar.gammas = "all"; + mesPar.mom = "0. 0. 0. 0."; application.createModule("3pt_Z2_" + std::to_string(t) + "_" diff --git a/tests/hadrons/Test_hadrons_rarekaon.cc b/tests/hadrons/Test_hadrons_rarekaon.cc new file mode 100644 index 00000000..89d7d501 --- /dev/null +++ b/tests/hadrons/Test_hadrons_rarekaon.cc @@ -0,0 +1,337 @@ +/******************************************************************************* + Grid physics library, www.github.com/paboyle/Grid + + Source file: tests/hadrons/Test_hadrons_rarekaon.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; + +enum quarks +{ + light = 0, + strange = 1, + charm = 2 +}; + +int main(int argc, char *argv[]) +{ + // parse command line ////////////////////////////////////////////////////// + std::string configStem; + + if (argc < 2) + { + std::cerr << "usage: " << argv[0] << " [Grid options]"; + std::cerr << std::endl; + std::exit(EXIT_FAILURE); + } + configStem = argv[1]; + + // 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 mass = {.01, .04, .2}; + std::vector flavour = {"l", "s", "c"}; + std::vector solvers = {"CG_l", "CG_s", "CG_c"}; + std::string kmom = "0. 0. 0. 0."; + std::string pmom = "1. 0. 0. 0."; + std::string qmom = "-1. 0. 0. 0."; + std::string mqmom = "1. 0. 0. 0."; + std::vector tKs = {0}; + unsigned int dt_pi = 16; + std::vector tJs = {8}; + unsigned int n_noise = 1; + unsigned int nt = 32; + bool do_disconnected(false); + + // 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 = 1000; + globalPar.genetic.maxCstGen = 200; + globalPar.genetic.popSize = 20; + globalPar.genetic.mutationRate = .1; + application.setPar(globalPar); + + // gauge field + if (configStem == "None") + { + application.createModule("gauge"); + } + else + { + MGauge::Load::Par gaugePar; + gaugePar.file = configStem; + application.createModule("gauge", gaugePar); + } + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 16; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + application.createModule("DWF_" + flavour[i], actionPar); + + // solvers + // RBPrecCG -> CG + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + application.createModule(solvers[i], + solverPar); + } + + // Create noise propagators for loops. + std::vector noiseSrcs; + std::vector> noiseRes; + std::vector> noiseProps; + if (n_noise > 0) + { + MSource::Z2::Par noisePar; + noisePar.tA = 0; + noisePar.tB = nt - 1; + std::string loop_stem = "loop_"; + + noiseRes.resize(flavour.size()); + noiseProps.resize(flavour.size()); + for (unsigned int nn = 0; nn < n_noise; ++nn) + { + std::string eta = INIT_INDEX("noise", nn); + application.createModule(eta, noisePar); + noiseSrcs.push_back(eta); + + for (unsigned int f = 0; f < flavour.size(); ++f) + { + std::string loop_prop = INIT_INDEX(loop_stem + flavour[f], nn); + std::string loop_res = loop_prop + "_res"; + makePropagator(application, loop_res, eta, solvers[f]); + makeLoop(application, loop_prop, eta, loop_res); + noiseRes[f].push_back(loop_res); + noiseProps[f].push_back(loop_prop); + } + } + } + + // Translate rare kaon decay across specified timeslices. + for (unsigned int i = 0; i < tKs.size(); ++i) + { + // Zero-momentum wall source propagators for kaon and pion. + unsigned int tK = tKs[i]; + unsigned int tpi = (tK + dt_pi) % nt; + std::string q_Kl_0 = INIT_INDEX("Q_l_0", tK); + std::string q_pil_0 = INIT_INDEX("Q_l_0", tpi); + MAKE_WALL_PROP(tK, q_Kl_0, solvers[light]); + MAKE_WALL_PROP(tpi, q_pil_0, solvers[light]); + + // Wall sources for kaon and pion with momentum insertion. If either + // p or k are zero, or p = k, re-use the existing name to avoid + // duplicating a propagator. + std::string q_Ks_k = INIT_INDEX("Q_Ks_k", tK); + std::string q_Ks_p = INIT_INDEX((kmom == pmom) ? "Q_Ks_k" : "Q_Ks_p", tK); + std::string q_pil_k = INIT_INDEX((kmom == ZERO_MOM) ? "Q_l_0" : "Q_l_k", tpi); + std::string q_pil_p = INIT_INDEX((pmom == kmom) ? q_pil_k : ((pmom == ZERO_MOM) ? "Q_l_0" : "Q_l_p"), tpi); + MAKE_3MOM_WALL_PROP(tK, kmom, q_Ks_k, solvers[strange]); + MAKE_3MOM_WALL_PROP(tK, pmom, q_Ks_p, solvers[strange]); + MAKE_3MOM_WALL_PROP(tpi, kmom, q_pil_k, solvers[light]); + MAKE_3MOM_WALL_PROP(tpi, pmom, q_pil_p, solvers[light]); + + /*********************************************************************** + * CONTRACTIONS: pi and K 2pt contractions with mom = p, k. + **********************************************************************/ + // Wall-Point + std::string PW_K_k = INIT_INDEX("PW_K_k", tK); + std::string PW_K_p = INIT_INDEX("PW_K_p", tK); + std::string PW_pi_k = INIT_INDEX("PW_pi_k", tpi); + std::string PW_pi_p = INIT_INDEX("PW_pi_p", tpi); + mesonContraction(application, 2, q_Kl_0, q_Ks_k, PW_K_k, kmom); + mesonContraction(application, 2, q_Kl_0, q_Ks_p, PW_K_p, pmom); + mesonContraction(application, 2, q_pil_k, q_pil_0, PW_pi_k, kmom); + mesonContraction(application, 2, q_pil_p, q_pil_0, PW_pi_p, pmom); + // Wall-Wall, to be done - requires modification of meson module. + + /*********************************************************************** + * CONTRACTIONS: 3pt Weak Hamiltonian, C & W (non-Eye type) classes. + **********************************************************************/ + std::string HW_CW_k = LABEL_3PT("HW_CW_k", tK, tpi); + std::string HW_CW_p = LABEL_3PT("HW_CW_p", tK, tpi); + weakContractionNonEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, q_pil_0, HW_CW_k); + weakContractionNonEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, q_pil_0, HW_CW_p); + + /*********************************************************************** + * CONTRACTIONS: 3pt sd insertion. + **********************************************************************/ + // Note: eventually will use wall sink smeared q_Kl_0 instead. + std::string sd_k = LABEL_3PT("sd_k", tK, tpi); + std::string sd_p = LABEL_3PT("sd_p", tK, tpi); + gamma3ptContraction(application, 3, q_Kl_0, q_Ks_k, q_pil_k, sd_k); + gamma3ptContraction(application, 3, q_Kl_0, q_Ks_p, q_pil_p, sd_p); + + for (unsigned int nn = 0; nn < n_noise; ++nn) + { + /******************************************************************* + * CONTRACTIONS: 3pt Weak Hamiltonian, S and E (Eye type) classes. + ******************************************************************/ + // Note: eventually will use wall sink smeared q_Kl_0 instead. + for (unsigned int f = 0; f < flavour.size(); ++f) + { + if ((f != strange) || do_disconnected) + { + std::string HW_SE_k = LABEL_3PT("HW_SE_k_" + flavour[f], tK, tpi); + std::string HW_SE_p = LABEL_3PT("HW_SE_p_" + flavour[f], tK, tpi); + std::string loop_q = noiseProps[f][nn]; + weakContractionEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, loop_q, HW_CW_k); + weakContractionEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, loop_q, HW_CW_p); + } + } + } + + // Perform separate contractions for each t_J position. + for (unsigned int j = 0; j < tJs.size(); ++j) + { + // Sequential sources for current insertions. Local for now, + // gamma_0 only. + unsigned int tJ = (tJs[j] + tK) % nt; + MSource::SeqGamma::Par seqPar; + std::string q_KlCl_q = LABEL_3PT("Q_KlCl_q", tK, tJ); + std::string q_KsCs_mq = LABEL_3PT("Q_KsCs_mq", tK, tJ); + std::string q_pilCl_q = LABEL_3PT("Q_pilCl_q", tpi, tJ); + std::string q_pilCl_mq = LABEL_3PT("Q_pilCl_mq", tpi, tJ); + MAKE_SEQUENTIAL_PROP(tJ, q_Kl_0, qmom, q_KlCl_q, solvers[light]); + MAKE_SEQUENTIAL_PROP(tJ, q_Ks_k, mqmom, q_KsCs_mq, solvers[strange]); + MAKE_SEQUENTIAL_PROP(tJ, q_pil_p, qmom, q_pilCl_q, solvers[light]); + MAKE_SEQUENTIAL_PROP(tJ, q_pil_0, mqmom, q_pilCl_mq, solvers[light]); + + /******************************************************************* + * CONTRACTIONS: pi and K 3pt contractions with current insertion. + ******************************************************************/ + // Wall-Point + std::string C_PW_Kl = LABEL_3PT("C_PW_Kl", tK, tJ); + std::string C_PW_Ksb = LABEL_3PT("C_PW_Ksb", tK, tJ); + std::string C_PW_pilb = LABEL_3PT("C_PW_pilb", tK, tJ); + std::string C_PW_pil = LABEL_3PT("C_PW_pil", tK, tJ); + mesonContraction(application, 3, q_KlCl_q, q_Ks_k, C_PW_Kl, pmom); + mesonContraction(application, 3, q_Kl_0, q_KsCs_mq, C_PW_Ksb, pmom); + mesonContraction(application, 3, q_pil_0, q_pilCl_q, C_PW_pilb, kmom); + mesonContraction(application, 3, q_pilCl_mq, q_pil_p, C_PW_pil, kmom); + // Wall-Wall, to be done. + + /******************************************************************* + * CONTRACTIONS: 4pt contractions, C & W classes. + ******************************************************************/ + std::string CW_Kl = LABEL_4PT("CW_Kl", tK, tJ, tpi); + std::string CW_Ksb = LABEL_4PT("CW_Ksb", tK, tJ, tpi); + std::string CW_pilb = LABEL_4PT("CW_pilb", tK, tJ, tpi); + std::string CW_pil = LABEL_4PT("CW_pil", tK, tJ, tpi); + weakContractionNonEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, q_pil_0, CW_Kl); + weakContractionNonEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, q_pil_0, CW_Ksb); + weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, q_pil_0, CW_pilb); + weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, q_pilCl_mq, CW_pil); + + /******************************************************************* + * CONTRACTIONS: 4pt contractions, sd insertions. + ******************************************************************/ + // Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead. + std::string sd_Kl = LABEL_4PT("sd_Kl", tK, tJ, tpi); + std::string sd_Ksb = LABEL_4PT("sd_Ksb", tK, tJ, tpi); + std::string sd_pilb = LABEL_4PT("sd_pilb", tK, tJ, tpi); + gamma3ptContraction(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, sd_Kl); + gamma3ptContraction(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, sd_Ksb); + gamma3ptContraction(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, sd_pilb); + + // Sequential sources for each noise propagator. + for (unsigned int nn = 0; nn < n_noise; ++nn) + { + std::string loop_stem = "loop_"; + + // Contraction required for each quark flavour - alternatively + // drop the strange loop if not performing disconnected + // contractions or neglecting H_W operators Q_3 -> Q_10. + for (unsigned int f = 0; f < flavour.size(); ++f) + { + if ((f != strange) || do_disconnected) + { + std::string eta = noiseSrcs[nn]; + std::string loop_q = noiseProps[f][nn]; + std::string loop_qCq = LABEL_3PT(loop_stem + flavour[f], tJ, nn); + std::string loop_qCq_res = loop_qCq + "_res"; + MAKE_SEQUENTIAL_PROP(tJ, noiseRes[f][nn], qmom, + loop_qCq_res, solvers[f]); + makeLoop(application, loop_qCq, eta, loop_qCq_res); + + /******************************************************* + * CONTRACTIONS: 4pt contractions, S & E classes. + ******************************************************/ + // Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead. + std::string SE_Kl = LABEL_4PT_NOISE("SE_Kl", tK, tJ, tpi, nn); + std::string SE_Ksb = LABEL_4PT_NOISE("SE_Ksb", tK, tJ, tpi, nn); + std::string SE_pilb = LABEL_4PT_NOISE("SE_pilb", tK, tJ, tpi, nn); + std::string SE_loop = LABEL_4PT_NOISE("SE_loop", tK, tJ, tpi, nn); + weakContractionEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, loop_q, SE_Kl); + weakContractionEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, loop_q, SE_Ksb); + weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, SE_pilb); + weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, loop_qCq, SE_loop); + + /******************************************************* + * CONTRACTIONS: 4pt contractions, pi0 disconnected + * loop. + ******************************************************/ + std::string disc0 = LABEL_4PT_NOISE("disc0", tK, tJ, tpi, nn); + disc0Contraction(application, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, disc0); + + /******************************************************* + * CONTRACTIONS: Disconnected loop. + ******************************************************/ + std::string discLoop = "disc_" + loop_qCq; + discLoopContraction(application, loop_qCq, discLoop); + } + } + } + } + } + // execution + std::string par_file_name = "rarekaon_000_100_tK0_tpi16_tJ8_noloop_mc0.2.xml"; + application.saveParameterFile(par_file_name); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +} diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index ea6e4479..2d731ff4 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -96,12 +96,16 @@ 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"; + mesPar.mom = "0. 0. 0. 0."; application.createModule("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"; + mesPar.mom = "0. 0. 0. 0."; application.createModule("meson_Z2_" + flavour[i] + flavour[j], mesPar); diff --git a/tests/qdpxx/Test_qdpxx_loops_staples.cc b/tests/qdpxx/Test_qdpxx_loops_staples.cc index 3bed9601..db384fcc 100644 --- a/tests/qdpxx/Test_qdpxx_loops_staples.cc +++ b/tests/qdpxx/Test_qdpxx_loops_staples.cc @@ -317,6 +317,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu) Grid::QCD::LatticeComplex rect(UGrid); Grid::QCD::TComplex trect; Grid::QCD::Complex crect; + Grid::RealD rrect; Grid::RealD vol = UGrid->gSites(); for(int mu=0;mu