From 03c3d495a26facb4715d33fdcc6d2724309d8506 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Fri, 12 Oct 2018 10:59:33 +0100 Subject: [PATCH] First cut (non functional NPR code) developed by Julia Kettle --- Hadrons/Modules/MGauge/GaugeFix.cc | 93 +++++++++ Hadrons/Modules/MGauge/GaugeFix.hpp | 81 ++++++++ Hadrons/Modules/MNPR/Amputate.cc | 8 + Hadrons/Modules/MNPR/Amputate.hpp | 199 +++++++++++++++++++ Hadrons/Modules/MNPR/Bilinear.cc | 8 + Hadrons/Modules/MNPR/Bilinear.hpp | 224 ++++++++++++++++++++++ Hadrons/Modules/MNPR/FourQuark.cc | 8 + Hadrons/Modules/MNPR/FourQuark.hpp | 273 +++++++++++++++++++++++++++ Hadrons/Modules/MSource/Momentum.cc | 8 + Hadrons/Modules/MSource/Momentum.hpp | 121 ++++++++++++ 10 files changed, 1023 insertions(+) create mode 100644 Hadrons/Modules/MGauge/GaugeFix.cc create mode 100644 Hadrons/Modules/MGauge/GaugeFix.hpp create mode 100644 Hadrons/Modules/MNPR/Amputate.cc create mode 100644 Hadrons/Modules/MNPR/Amputate.hpp create mode 100644 Hadrons/Modules/MNPR/Bilinear.cc create mode 100644 Hadrons/Modules/MNPR/Bilinear.hpp create mode 100644 Hadrons/Modules/MNPR/FourQuark.cc create mode 100644 Hadrons/Modules/MNPR/FourQuark.hpp create mode 100644 Hadrons/Modules/MSource/Momentum.cc create mode 100644 Hadrons/Modules/MSource/Momentum.hpp diff --git a/Hadrons/Modules/MGauge/GaugeFix.cc b/Hadrons/Modules/MGauge/GaugeFix.cc new file mode 100644 index 00000000..bd57ed07 --- /dev/null +++ b/Hadrons/Modules/MGauge/GaugeFix.cc @@ -0,0 +1,93 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Fix.cc + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Antonin Portelli + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MGauge; + + +/****************************************************************************** +* TGaugeFix implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TGaugeFix::TGaugeFix(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TGaugeFix::getInput(void) +{ + std::vector in = {par().gauge}; + return in; +} + +template +std::vector TGaugeFix::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TGaugeFix::setup(void) +{ + envCreateLat(GaugeField, getName()); +} + + +// execution /////////////////////////////////////////////////////////////////// +template +void TGaugeFix::execute(void) +//Loads the gauge and fixes it +{ + + std::cout << "executing" << std::endl; + LOG(Message) << "Fixing the Gauge" << std::endl; + LOG(Message) << par().gauge << std::endl; + auto &U = envGet(GaugeField, par().gauge); + auto &Umu = envGet(GaugeField, getName()); + LOG(Message) << "Gauge Field fetched" << std::endl; + //do we allow maxiter etc to be user set? + Real alpha = par().alpha; + int maxiter = par().maxiter; + Real Omega_tol = par().Omega_tol; + Real Phi_tol = par().Phi_tol; + bool Fourier = par().Fourier; + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(U,alpha,maxiter,Omega_tol,Phi_tol,Fourier); + Umu = U; + LOG(Message) << "Gauge Fixed" << std::endl; + +} + +template class Grid::Hadrons::MGauge::TGaugeFix; diff --git a/Hadrons/Modules/MGauge/GaugeFix.hpp b/Hadrons/Modules/MGauge/GaugeFix.hpp new file mode 100644 index 00000000..154237a7 --- /dev/null +++ b/Hadrons/Modules/MGauge/GaugeFix.hpp @@ -0,0 +1,81 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MGauge/Fix.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Antonin Portelli + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_MGaugeFix_hpp_ +#define Hadrons_MGaugeFix_hpp_ + +#include +#include +#include +#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Fix gauge * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MGauge) + +class GaugeFixPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(GaugeFixPar, + std::string, gauge, + Real, alpha, + int, maxiter, + Real, Omega_tol, + Real, Phi_tol, + bool, Fourier); +}; + +template +class TGaugeFix: public Module +{ +public: + GAUGE_TYPE_ALIASES(GImpl,); +public: + // constructor + TGaugeFix(const std::string name); + // destructor + virtual ~TGaugeFix(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(GaugeFix, TGaugeFix, MGauge); + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MGaugeFix_hpp_ diff --git a/Hadrons/Modules/MNPR/Amputate.cc b/Hadrons/Modules/MNPR/Amputate.cc new file mode 100644 index 00000000..c4ada9d9 --- /dev/null +++ b/Hadrons/Modules/MNPR/Amputate.cc @@ -0,0 +1,8 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TAmputate; + diff --git a/Hadrons/Modules/MNPR/Amputate.hpp b/Hadrons/Modules/MNPR/Amputate.hpp new file mode 100644 index 00000000..797247b5 --- /dev/null +++ b/Hadrons/Modules/MNPR/Amputate.hpp @@ -0,0 +1,199 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MNpr/Amputate.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk + +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_Amputate_hpp_ +#define Hadrons_Amputate_hpp_ + +#include +#include +#include +#include +//#include +//#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TAmputate * + Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class AmputatePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(AmputatePar, + std::string, Sin, //need to make this a propogator type? + std::string, Sout, //same + std::string, vertex, + std::string, pin, + std::string, pout, + std::string, output, + std::string, input); +}; + +template +class TAmputate: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, Vamp, + ); + }; +public: + // constructor + TAmputate(const std::string name); + // destructor + virtual ~TAmputate(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + virtual SpinColourMatrix invertspincolmat(SpinColourMatrix &scmat); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Amputate, ARG(TAmputate), MNPR); + +/****************************************************************************** + * TAmputate implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TAmputate::TAmputate(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TAmputate::getInput(void) +{ + std::vector input = {par().Sin, par().Sout, par().vertex}; + + return input; +} + +template +std::vector TAmputate::getOutput(void) +{ + std::vector output = {getName()}; + + + return output; +} + +// Invert spin colour matrix using Eigen +template +SpinColourMatrix TAmputate::invertspincolmat(SpinColourMatrix &scmat) +{ + Eigen::MatrixXcf scmat_2d(Ns*Nc,Ns*Nc); + for(int ic=0; ic +void TAmputate::execute(void) +{ + LOG(Message) << "Computing bilinear amputations '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + BinaryWriter writer(par().output); + PropagatorField1 &Sin = *env().template getObject(par().Sin); //Do these have the phases taken into account?? Don't think so. FIX + PropagatorField2 &Sout = *env().template getObject(par().Sout); + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + std::vector latt_size(pin.begin(), pin.end()); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + LOG(Message) << "Propagators set up " << std::endl; + std::vector vertex; // Let's read from file here + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + LOG(Message) << "reading file - " << par().input << std::endl; + BinaryReader reader(par().input); + Complex Ci(0.0,1.0); + + std::string svertex; + read(reader,"vertex", vertex); + LOG(Message) << "vertex read" << std::endl; + + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + SpinColourMatrix Sin_mom = sum(Sin); + SpinColourMatrix Sout_mom = sum(Sout); + LOG(Message) << "summed over lattice" << std::endl; + + LOG(Message) << "Lattice -> spincolourmatrix conversion" << std::endl; + + SpinColourMatrix Sin_inv = invertspincolmat(Sin_mom); + SpinColourMatrix Sout_inv = invertspincolmat(Sout_mom); + LOG(Message) << "Inversions done" << std::endl; + + result.Vamp.resize(Gamma::nGamma/2); + for( int mu=0; mu < Gamma::nGamma/2; mu++){ + Gamma::Algebra gam = mu; + result.Vamp[mu] = 1/12.0*trace(adj(Gamma(mu*2+1))*g5*Sout_inv*g5*vertex[mu]*Sin_inv); + LOG(Message) << "Vamp[" << mu << "] - " << result.Vamp[mu] << std::endl; + } +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_Amputate_hpp_ diff --git a/Hadrons/Modules/MNPR/Bilinear.cc b/Hadrons/Modules/MNPR/Bilinear.cc new file mode 100644 index 00000000..01fb08c6 --- /dev/null +++ b/Hadrons/Modules/MNPR/Bilinear.cc @@ -0,0 +1,8 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TBilinear; + diff --git a/Hadrons/Modules/MNPR/Bilinear.hpp b/Hadrons/Modules/MNPR/Bilinear.hpp new file mode 100644 index 00000000..0dc621b9 --- /dev/null +++ b/Hadrons/Modules/MNPR/Bilinear.hpp @@ -0,0 +1,224 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MContraction/Bilinear.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk + +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_Bilinear_hpp_ +#define Hadrons_Bilinear_hpp_ + +#include +#include +#include +#include +//#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TBilinear * + Performs bilinear contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta in Rome-Southampton NPR +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class BilinearPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(BilinearPar, + std::string, Sin, + std::string, Sout, + std::string, pin, + std::string, pout, + std::string, output); +}; + +template +class TBilinear: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, bilinear); + }; +public: + // constructor + TBilinear(const std::string name); + // destructor + virtual ~TBilinear(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + //LatticeSpinColourMatrix PhaseProps(LatticeSpinColourMatrix S, std::vector p); + // setup + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(Bilinear, ARG(TBilinear), MNPR); + +/****************************************************************************** + * TBilinear implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TBilinear::TBilinear(const std::string name) +: Module(name) +{} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TBilinear::setup(void) +{ + //env().template registerLattice(getName()); + //env().template registerObject(getName()); +} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TBilinear::getInput(void) +{ + std::vector input = {par().Sin, par().Sout}; + + return input; +} + +template +std::vector TBilinear::getOutput(void) +{ + std::vector out = {getName()}; + + return out; +} + +/* +/////Phase propagators////////////////////////// +template +LatticeSpinColourMatrix TBilinear::PhaseProps(LatticeSpinColourMatrix S, std::vector p) +{ + GridBase *grid = S._grid; + LatticeComplex pdotx(grid), coor(grid); + std::vector latt_size = grid->_fdimensions; + Complex Ci(0.0,1.0); + pdotx=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotx = pdotx +(TwoPiL * p[mu]) * coor; + } + S = S*exp(-Ci*pdotx); + return S; +} +*/ +// execution /////////////////////////////////////////////////////////////////// +template +void TBilinear::execute(void) +{ +/************************************************************************** + +Compute the bilinear vertex needed for the NPR. +V(G) = sum_x [ g5 * adj(S'(x,p2)) * g5 * G * S'(x,p1) ]_{si,sj,ci,cj} +G is one of the 16 gamma vertices [I,gmu,g5,g5gmu,sig(mu,nu)] + + * G + / \ + p1/ \p2 + / \ + / \ + +Returns a spin-colour matrix, with indices si,sj, ci,cj + +Conventions: +p1 - incoming momenta +p2 - outgoing momenta +q = (p1-p2) +**************************************************************************/ + + LOG(Message) << "Computing bilinear contractions '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + + BinaryWriter writer(par().output); + + + // Propogators + LatticeSpinColourMatrix &Sin = *env().template getObject(par().Sin); + LatticeSpinColourMatrix &Sout = *env().template getObject(par().Sout); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + // momentum on legs + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + std::vector latt_size(pin.begin(), pin.end()); + //bilinears + LatticeSpinColourMatrix bilinear_x(env().getGrid()); + SpinColourMatrix bilinear; + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + Complex Ci(0.0,1.0); + + // + + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + ////Set up gamma vector////////////////////////// + std::vector gammavector; + for( int i=0; i + +using namespace Grid; +using namespace Hadrons; +using namespace MNPR; + +template class Grid::Hadrons::MNPR::TFourQuark; + diff --git a/Hadrons/Modules/MNPR/FourQuark.hpp b/Hadrons/Modules/MNPR/FourQuark.hpp new file mode 100644 index 00000000..d71ddf43 --- /dev/null +++ b/Hadrons/Modules/MNPR/FourQuark.hpp @@ -0,0 +1,273 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: extras/Hadrons/Modules/MNPR/FourQuark.hpp + +Copyright (C) 2015 +Copyright (C) 2016 + +Author: Julia Kettle J.R.Kettle-2@sms.ed.ac.uk + +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_FourQuark_hpp_ +#define Hadrons_FourQuark_hpp_ + +#include +#include +#include +#include +#include +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * TFourQuark * + Performs fourquark contractions of the type tr[g5*adj(Sout)*g5*G*Sin] + Suitable for non exceptional momenta +******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MNPR) + +class FourQuarkPar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(FourQuarkPar, + std::string, Sin, //need to make this a propogator type? + std::string, Sout, //same + std::string, pin, + std::string, pout, + bool, fullbasis, + std::string, output); +}; + +template +class TFourQuark: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, fourquark); + }; +public: + // constructor + TFourQuark(const std::string name); + // destructor + virtual ~TFourQuark(void) = default; + // dependencies/products + virtual std::vector getInput(void); + virtual std::vector getOutput(void); + // setup + virtual void tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b); + virtual void setup(void); + // execution + virtual void execute(void); +}; + +MODULE_REGISTER_TMP(FourQuark, ARG(TFourQuark), MNPR); + +/****************************************************************************** + * TFourQuark implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TFourQuark::TFourQuark(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TFourQuark::getInput(void) +{ + std::vector input = {par().Sin, par().Sout}; + + return input; +} + +template +std::vector TFourQuark::getOutput(void) +{ + std::vector output = {getName()}; + + return output; +} + + +template +void TFourQuark::tensorprod(LatticeSpinColourSpinColourMatrix &lret, LatticeSpinColourMatrix a, LatticeSpinColourMatrix b) +{ +#if 0 + parallel_for(auto site=lret.begin();site +void TFourQuark::setup(void) +{ + envCreateLat(LatticeSpinColourMatrix, getName()); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TFourQuark::execute(void) +{ + +/********************************************************************************* + +TFourQuark : Creates the four quark vertex required for the NPR of four-quark ops + +V_{Gamma_1,Gamma_2} = sum_x [ ( g5 * adj(S'(x,p2)) * g5 * G1 * S'(x,p1) )_ci,cj;si,sj x ( g5 * adj(S'(x,p2)) * g5 * G2 S'(x,p1) )_ck,cl;sk,cl ] + +Create a bilinear vertex for G1 and G2 the spin and colour indices are kept free. Where there are 16 potential Gs. +We then find the outer product of V1 and V2, keeping the spin and colour indices uncontracted +Then this is summed over the lattice coordinate +Result is a SpinColourSpinColourMatrix - with 4 colour and 4 spin indices. +We have up to 256 of these including the offdiag (G1 != G2). + + \ / + \p1 p1/ + \ / + \ / + G1 * * G2 + / \ + / \ + /p2 p2\ + / \ + +*********************************************************************************/ + + + + + LOG(Message) << "Computing fourquark contractions '" << getName() << "' using" + << " momentum '" << par().Sin << "' and '" << par().Sout << "'" + << std::endl; + + BinaryWriter writer(par().output); + + PropagatorField1 &Sin = *env().template getObject(par().Sin); + PropagatorField2 &Sout = *env().template getObject(par().Sout); + std::vector pin = strToVec(par().pin), pout = strToVec(par().pout); + bool fullbasis = par().fullbasis; + Gamma g5(Gamma::Algebra::Gamma5); + Result result; + std::vector latt_size(pin.begin(), pin.end()); + LatticeComplex pdotxin(env().getGrid()), pdotxout(env().getGrid()), coor(env().getGrid()); + LatticeSpinColourMatrix bilinear_mu(env().getGrid()), bilinear_nu(env().getGrid()); + LatticeSpinColourSpinColourMatrix lret(env().getGrid()); + Complex Ci(0.0,1.0); + + //Phase propagators + //Sin = Grid::QCD::PropUtils::PhaseProps(Sin,pin); + //Sout = Grid::QCD::PropUtils::PhaseProps(Sout,pout); + + //find p.x for in and out so phase can be accounted for in propagators + pdotxin=zero; + pdotxout=zero; + for (unsigned int mu = 0; mu < 4; ++mu) + { + Real TwoPiL = M_PI * 2.0/ latt_size[mu]; + LatticeCoordinate(coor,mu); + pdotxin = pdotxin +(TwoPiL * pin[mu]) * coor; + pdotxout= pdotxout +(TwoPiL * pout[mu]) * coor; + } + Sin = Sin*exp(-Ci*pdotxin); //phase corrections + Sout = Sout*exp(-Ci*pdotxout); + + + //Set up Gammas + std::vector gammavector; + for( int i=1; i + +using namespace Grid; +using namespace Hadrons; +using namespace MSource; + +template class Grid::Hadrons::MSource::TMomSource; + diff --git a/Hadrons/Modules/MSource/Momentum.hpp b/Hadrons/Modules/MSource/Momentum.hpp new file mode 100644 index 00000000..7e4b3c09 --- /dev/null +++ b/Hadrons/Modules/MSource/Momentum.hpp @@ -0,0 +1,121 @@ +#ifndef Hadrons_MomSource_hpp_ +#define Hadrons_MomSource_hpp_ + +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/* +Plane Wave source +----------------- +src_x = e^i2pi/L * p *position +*/ + +/****************************************************************************** + * TPlane * + ******************************************************************************/ + BEGIN_MODULE_NAMESPACE(MSource) + + class MomSourcePar: Serializable + { + public: +//What is meant by serializable in this context +GRID_SERIALIZABLE_CLASS_MEMBERS(MomSourcePar, +std::string, mom); +}; + + +template +class TMomSource: public Module +{ +public: +FERM_TYPE_ALIASES(FImpl,); +public: +// constructor +TMomSource(const std::string name); +// destructor +virtual ~TMomSource(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_TMP(MomSource, TMomSource, MSource); +//MODULE_REGISTER_NS(MomSource, TMomSource, MSource); + +/****************************************************************************** +* TMomSource template implementation * +******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TMomSource::TMomSource(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TMomSource::getInput(void) +{ + std::vector in; + return in; +} + +template +std::vector TMomSource::getOutput(void) +{ + std::vector out = {getName()}; + return out; +} + + +// setup /////////////////////////////////////////////////////////////////////// +template +void TMomSource::setup(void) +{ + envCreateLat(PropagatorField, getName()); +} + + +//execution////////////////////////////////////////////////////////////////// +template +void TMomSource::execute(void) +{ + LOG(Message) << "Generating planewave momentum source with momentum " << par().mom << std::endl; + //what does this env do? + PropagatorField &src = envGet(PropagatorField, getName()); + Lattice> t(env().getGrid()); + LatticeComplex C(env().getGrid()), coor(env().getGrid()); + std::vector p; + std::vector latt_size(GridDefaultLatt().begin(), GridDefaultLatt().end()); + Complex i(0.0,1.0); + + LOG(Message) << " " << std::endl; + //get the momentum from parameters + p = strToVec(par().mom); + C = zero; + LOG(Message) << "momentum converted from string - " << std::to_string(p[0]) <