From e66d48c14246f3bf367283816c146c3c19c671f5 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 30 Jul 2019 13:46:59 +0100 Subject: [PATCH] second way to compute baryons - qdp style --- Grid/qcd/utils/BaryonUtils.h | 58 ++++++- Hadrons/Modules.hpp | 2 + Hadrons/Modules/MContraction/Baryon.hpp | 86 +--------- Hadrons/Modules/MContraction/Baryon2.cc | 35 ++++ Hadrons/Modules/MContraction/Baryon2.hpp | 206 +++++++++++++++++++++++ Hadrons/modules.inc | 4 +- 6 files changed, 306 insertions(+), 85 deletions(-) create mode 100644 Hadrons/Modules/MContraction/Baryon2.cc create mode 100644 Hadrons/Modules/MContraction/Baryon2.hpp diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index ea2034de..468c1f78 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -27,6 +27,9 @@ public: const Gamma GammaA, const Gamma GammaB, ComplexField &baryon_corr); + + static LatticeSpinColourMatrix quarkContract13(const PropagatorField &q1, + const PropagatorField &q2); }; @@ -139,19 +142,70 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1, }}} } - if (ie_src==0 && ie_snk==0){ + /*if (ie_src==0 && ie_snk==0){ baryon_corr._odata[ss] = result; } else { baryon_corr._odata[ss] += result; - } + }*/ } } + baryon_corr._odata[ss] = result; + } //end loop over lattice sites } +//QDP / CHROMA - style diquark construction +// (q_out)^{c'c}_{alpha,beta} = epsilon^{abc} epsilon^{a'b'c'} (q1)^{aa'}_{rho alpha}^* (q2)^{bb'}_{rho beta} +template +LatticeSpinColourMatrix BaryonUtils::quarkContract13(const PropagatorField &q1, + const PropagatorField &q2) +{ + GridBase *grid = q1._grid; + + + std::vector> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + std::vector epsilon_sgn = {1,1,1,-1,-1,-1}; + std::vector wick_contraction = {0,0,0,0,0,0}; + + LatticeSpinColourMatrix q_out=zero; + + parallel_for(int ss=0;ssoSites();ss++){ + + typedef typename ComplexField::vector_object vobj; + + auto D1 = q1._odata[ss]; + auto D2 = q2._odata[ss]; + //auto D_out = q_out._odata[ss]; + //D_out=zero; + + SpinColourMatrix D_out=zero; + + for (int ie_src=0; ie_src < 6 ; ie_src++){ + int a_src = epsilon[ie_src][0]; //a + int b_src = epsilon[ie_src][1]; //b + int c_src = epsilon[ie_src][2]; //c + for (int ie_snk=0; ie_snk < 6 ; ie_snk++){ + int a_snk = epsilon[ie_snk][0]; //a' + int b_snk = epsilon[ie_snk][1]; //b' + int c_snk = epsilon[ie_snk][2]; //c' + for (int alpha=0; alpha #include #include +#include #include #include #include +#include #include #include #include diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 52508864..f78759d8 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -164,13 +164,13 @@ void TBaryon::setup(void) template void TBaryon::execute(void) { - LOG(Message) << "Computing nucleon contractions '" << getName() << "' using" - << " quarks '" << par().q1 << "', '" << par().q2 << "', and '" - << par().q3 << "'" << std::endl; + LOG(Message) << "Computing baryon contractions '" << getName() << "' using" + << " quarks '" << par().q1 << "', and a diquark formed of ('" << par().q2 << "', and '" + << par().q3 << "')" << std::endl; auto &q1 = envGet(PropagatorField1, par().q1); auto &q2 = envGet(PropagatorField2, par().q2); - auto &q3 = envGet(PropagatorField3, par().q2); + auto &q3 = envGet(PropagatorField3, par().q3); envGetTmp(LatticeComplex, c); envGetTmp(LatticeComplex, diquark); Result result; @@ -178,84 +178,6 @@ void TBaryon::execute(void) result.corr.resize(nt); const std::string gamma{ par().gamma }; std::vector buf; - // C = i gamma_2 gamma_4 => C gamma_5 = - i gamma_1 gamma_3 -/* Gamma GammaA(Gamma::Algebra::Identity); //Still hardcoded 1 - Gamma GammaB(Gamma::Algebra::SigmaXZ); //Still hardcoded Cg5 - Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) - - std::vector> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; - std::vector epsilon_sgn = {1,1,1,-1,-1,-1}; - - char left[] = "uud"; - char right[] = "uud"; - std::vector wick_contraction = {0,0,0,0,0,0}; - - for (int ie=0; ie < 6 ; ie++) - if (left[0] == right[epsilon[ie][0]] && left[1] == right[epsilon[ie][1]] && left[2] == right[epsilon[ie][2]]) - wick_contraction[ie]=1; - - - int parity = 1; - - - for (int ie_src=0; ie_src < 6 ; ie_src++){ - int a_src = epsilon[ie_src][0]; //a - int b_src = epsilon[ie_src][1]; //b - int c_src = epsilon[ie_src][2]; //c - for (int ie_snk=0; ie_snk < 6 ; ie_snk++){ - int a_snk = epsilon[ie_snk][0]; //a' - int b_snk = epsilon[ie_snk][1]; //b' - int c_snk = epsilon[ie_snk][2]; //c' - auto Daa = peekColour(q2,a_snk,a_src); //D_{alpha' alpha} - auto Dbb = peekColour(q3,b_snk,b_src); //D_{beta' beta} - auto Dcc = peekColour(q1,c_snk,c_src); //D_{gamma' gamma} - auto Dab = peekColour(q2,a_snk,b_src); //D_{alpha' beta} - auto Dac = peekColour(q2,a_snk,c_src); //D_{alpha' gamma} - auto Dba = peekColour(q3,b_snk,a_src); //D_{beta' alpha} - auto Dbc = peekColour(q3,b_snk,c_src); //D_{beta' gamma} - auto Dca = peekColour(q1,c_snk,a_src); //D_{gamma' alpha} - auto Dcb = peekColour(q1,c_snk,b_src); //D_{gamma' beta} - // This is the \delta_{123}^{123} part - if (wick_contraction[0]){ - diquark = trace(GammaB * Daa * GammaB * Dbb); //1st GammaB and Daa transposed???? - auto temp = GammaA * Dcc * diquark; - auto g4_temp = GammaA * g4 * temp; - c += epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * 0.5 * trace(GammaA * temp + (double)parity * g4_temp); - } - // This is the \delta_{123}^{231} part - if (wick_contraction[1]){ - auto temp = GammaA * Dca * GammaB * Dab * GammaB * Dbc; //Dab transposed??? - auto g4_temp = GammaA * g4 * temp; - c += epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * 0.5 * trace(GammaA * temp + (double)parity * g4_temp); - } - // This is the \delta_{123}^{312} part - if (wick_contraction[2]){ - auto temp = GammaA * Dcb * GammaB * Dba * GammaB * Dac; //both GammaB and Dba transposed??? - auto g4_temp = GammaA * g4 * temp; - c += epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * 0.5 * trace(GammaA * temp + (double)parity * g4_temp); - } - // This is the \delta_{123}^{132} part - if (wick_contraction[3]){ - diquark = trace(GammaB * Dba * GammaB * Dab); //2nd GammaB and Dab transposed???? - auto temp = GammaA * Dcc * diquark; - auto g4_temp = GammaA * g4 * temp; - c -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * 0.5 * trace(GammaA * temp + (double)parity * g4_temp); - } - // This is the \delta_{123}^{321} part - if (wick_contraction[4]){ - auto temp = GammaA * Dcb * GammaB * Daa * GammaB * Dbc; //1st GammaB and Daa transposed??? - auto g4_temp = GammaA * g4 * temp; - c -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * 0.5 * trace(GammaA * temp + (double)parity * g4_temp); - } - // This is the \delta_{123}^{213} part - if (wick_contraction[5]){ - auto temp = GammaA * Dca * GammaB * Dbb * GammaB * Dac; //(Dbb*GammaB) transposed??? - auto g4_temp = GammaA * g4 * temp; - c -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * 0.5 * trace(GammaA * temp + (double)parity * g4_temp); - } - } - } -*/ const Gamma GammaA{ Gamma::Algebra::Identity }; const Gamma GammaB{ al }; diff --git a/Hadrons/Modules/MContraction/Baryon2.cc b/Hadrons/Modules/MContraction/Baryon2.cc new file mode 100644 index 00000000..b8d3bc52 --- /dev/null +++ b/Hadrons/Modules/MContraction/Baryon2.cc @@ -0,0 +1,35 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/Baryon2.cc + +Copyright (C) 2015-2019 + +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 MContraction; + +template class Grid::Hadrons::MContraction::TBaryon2; + diff --git a/Hadrons/Modules/MContraction/Baryon2.hpp b/Hadrons/Modules/MContraction/Baryon2.hpp new file mode 100644 index 00000000..df928472 --- /dev/null +++ b/Hadrons/Modules/MContraction/Baryon2.hpp @@ -0,0 +1,206 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/Baryon2.hpp + +Copyright (C) 2015-2019 + +Author: Antonin Portelli +Author: Felix Erben + +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_MContraction_Baryon2_hpp_ +#define Hadrons_MContraction_Baryon2_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * Baryon2 * + ******************************************************************************/ +BEGIN_MODULE_NAMESPACE(MContraction) + +class Baryon2Par: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Baryon2Par, + std::string, q1, + std::string, q2, + std::string, q3, + std::string, gamma, + std::string, output); +}; + +template +class TBaryon2: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl1, 1); + FERM_TYPE_ALIASES(FImpl2, 2); + FERM_TYPE_ALIASES(FImpl3, 3); + class Result: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + std::vector, corr); + }; +public: + // constructor + TBaryon2(const std::string name); + // destructor + virtual ~TBaryon2(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); +protected: + // setup + virtual void setup(void); + // execution + virtual void execute(void); + // Which gamma algebra was specified + Gamma::Algebra al; +}; + +MODULE_REGISTER_TMP(Baryon2, ARG(TBaryon2), MContraction); + +/****************************************************************************** + * TBaryon2 implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TBaryon2::TBaryon2(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TBaryon2::getInput(void) +{ + std::vector input = {par().q1, par().q2, par().q3}; + + return input; +} + +template +std::vector TBaryon2::getOutput(void) +{ + std::vector out = {}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TBaryon2::setup(void) +{ + envTmpLat(LatticeComplex, "c"); + envTmpLat(LatticeComplex, "diquark"); + // Translate the full string naming the desired gamma structure into the one we need to use + const std::string gamma{ par().gamma }; + int iGamma = 0; + do + { + const char * pGammaName = Gamma::name[iGamma]; + int iLen = 0; + while( pGammaName[iLen] && pGammaName[iLen] != ' ' ) + iLen++; + if( !gamma.compare( 0, gamma.size(), pGammaName, iLen ) ) + break; + } + while( ++iGamma < Gamma::nGamma ); + if( iGamma >= Gamma::nGamma ) { + LOG(Message) << "Unrecognised gamma structure \"" << gamma << "\"" << std::endl; + assert( 0 && "Invalid gamma structure specified" ); + } + switch( iGamma ) { + case Gamma::Algebra::GammaX: + std::cout << "using interpolator C gamma_X" << std::endl; + al = Gamma::Algebra::GammaZGamma5; //Still hardcoded CgX = i gamma_3 gamma_5 + break; + case Gamma::Algebra::GammaY: + std::cout << "using interpolator C gamma_Y" << std::endl; + al = Gamma::Algebra::GammaT; //Still hardcoded CgX = - gamma_4 + break; + case Gamma::Algebra::GammaZ: + std::cout << "using interpolator C gamma_Z" << std::endl; + al = Gamma::Algebra::GammaXGamma5; //Still hardcoded CgX = i gamma_1 gamma_5 + break; + default: + { + LOG(Message) << "Unsupported gamma structure " << gamma << " = " << iGamma << std::endl; + assert( 0 && "Unsupported gamma structure" ); + // or you could do something like + al = static_cast( iGamma ); + break; + } + } + LOG(Message) << "Gamma structure " << gamma << " = " << iGamma + << " translated to " << Gamma::name[al] << std::endl; +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TBaryon2::execute(void) +{ + LOG(Message) << "Computing baryon contractions '" << getName() << "' using" + << " quarks '" << par().q1 << "', and a diquark formed of ('" << par().q2 << "', and '" + << par().q3 << "')" << std::endl; + + auto &q1 = envGet(PropagatorField1, par().q1); + auto &q2 = envGet(PropagatorField2, par().q2); + auto &q3 = envGet(PropagatorField3, par().q3); + envGetTmp(LatticeComplex, c); + //envGetTmp(LatticeComplex, diquark); + Result result; + int nt = env().getDim(Tp); + result.corr.resize(nt); + const std::string gamma{ par().gamma }; + std::vector buf; + + const Gamma GammaA{ Gamma::Algebra::Identity }; + const Gamma GammaB{ al }; + + LatticeSpinColourMatrix diquark; + + diquark = BaryonUtils::quarkContract13(q2*GammaB,GammaB*q3); + + //result = trace(GammaA*GammaA * traceColour(q1*traceSpin(diquark))) + 2.0 * trace(GammaA*GammaA*traceColour(q1*diquark)); + result = trace(q1*diquark); + + sliceSum(c,buf,Tp); + + for (unsigned int t = 0; t < buf.size(); ++t) + { + result.corr[t] = TensorRemove(buf[t]); + } + + saveResult(par().output, "baryon", result); +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MContraction_Baryon2_hpp_ diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index bcab1232..a18afd53 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -55,6 +55,7 @@ modules_cc =\ Modules/MContraction/WeakEye3pt.cc \ Modules/MContraction/Meson.cc \ Modules/MContraction/A2AAslashField.cc \ + Modules/MContraction/Baryon2.cc \ Modules/MContraction/Baryon.cc \ Modules/MContraction/Nucleon.cc \ Modules/MContraction/WeakNonEye3pt.cc \ @@ -63,7 +64,6 @@ modules_cc =\ Modules/MContraction/A2AMesonField.cc \ Modules/MContraction/A2ALoop.cc \ Modules/MContraction/Gamma3pt.cc \ - Modules/MContraction/Nucleon.cc \ Modules/MAction/MobiusDWF.cc \ Modules/MAction/WilsonClover.cc \ Modules/MAction/Wilson.cc \ @@ -138,9 +138,11 @@ modules_hpp =\ Modules/MContraction/WeakMesonDecayKl2.hpp \ Modules/MContraction/Nucleon.hpp \ Modules/MContraction/A2AAslashField.hpp \ + Modules/MContraction/Baryon2.hpp \ Modules/MContraction/WeakEye3pt.hpp \ Modules/MContraction/WeakNonEye3pt.hpp \ Modules/MContraction/Baryon.hpp \ + Modules/MContraction/Baryon_old.hpp \ Modules/MContraction/Meson.hpp \ Modules/MContraction/A2ALoop.hpp \ Modules/MContraction/Gamma3pt.hpp \