diff --git a/extras/Hadrons/Modules.hpp b/extras/Hadrons/Modules.hpp index 77ae08b7..e8cdf2af 100644 --- a/extras/Hadrons/Modules.hpp +++ b/extras/Hadrons/Modules.hpp @@ -30,11 +30,16 @@ 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 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/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..f9b0481f --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp @@ -0,0 +1,86 @@ +/************************************************************************************* + +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); +}; + +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..2b855ce0 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp @@ -0,0 +1,83 @@ +/************************************************************************************* + +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) + +class TWeakHamiltonianEye: public Module +{ +public: + TYPE_ALIASES(FIMPL,) + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, name, + std::vector, corr); + }; +public: + // constructor + TWeakHamiltonianEye(const std::string name); + // destructor + virtual ~TWeakHamiltonianEye(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(WeakHamiltonianEye, TWeakHamiltonianEye, MContraction); + +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..b7fbdf28 --- /dev/null +++ b/extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp @@ -0,0 +1,82 @@ +/************************************************************************************* + +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) + +class TWeakHamiltonianNonEye: public Module +{ +public: + TYPE_ALIASES(FIMPL,) + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::string, name, + std::vector, corr); + }; +public: + // constructor + TWeakHamiltonianNonEye(const std::string name); + // destructor + virtual ~TWeakHamiltonianNonEye(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(WeakHamiltonianNonEye, TWeakHamiltonianNonEye, MContraction); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_WeakHamiltonianNonEye_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.inc b/extras/Hadrons/modules.inc index 4251ffa3..3e2f1020 100644 --- a/extras/Hadrons/modules.inc +++ b/extras/Hadrons/modules.inc @@ -1,4 +1,6 @@ modules_cc =\ + Modules/MContraction/WeakHamiltonianEye.cc \ + Modules/MContraction/WeakHamiltonianNonEye.cc \ Modules/MGauge/Load.cc \ Modules/MGauge/Random.cc \ Modules/MGauge/Unit.cc @@ -8,12 +10,17 @@ modules_hpp =\ Modules/MAction/Wilson.hpp \ Modules/MContraction/Baryon.hpp \ Modules/MContraction/Meson.hpp \ + Modules/MContraction/WeakHamiltonian.hpp \ + Modules/MContraction/WeakHamiltonianEye.hpp \ + Modules/MContraction/WeakHamiltonianNonEye.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/lib/serialisation/BaseIO.h b/lib/serialisation/BaseIO.h index 0357915d..36e6fd77 100644 --- a/lib/serialisation/BaseIO.h +++ b/lib/serialisation/BaseIO.h @@ -138,6 +138,54 @@ namespace Grid { unsigned int dimInd_{0}; }; + // Pair IO utilities ///////////////////////////////////////////////////////// + // helper function to parse input in the format "" + 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/tests/IO/Test_serialisation.cc b/tests/IO/Test_serialisation.cc index 8204b05b..acafc900 100644 --- a/tests/IO/Test_serialisation.cc +++ b/tests/IO/Test_serialisation.cc @@ -113,6 +113,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); @@ -120,6 +121,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; @@ -127,21 +130,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/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_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);