From dab8c01c3d254fd9db1cf3b6d3948858db6a15a7 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Mon, 16 Sep 2019 17:20:54 +0100 Subject: [PATCH 01/18] added Baryon code --- Grid/qcd/utils/BaryonUtils.h | 176 +++++++ Hadrons/Modules/MContraction/Baryon.hpp | 61 ++- tests/qdpxx/Test_qdpxx_baryon.cc | 605 ++++++++++++++++++++++++ 3 files changed, 826 insertions(+), 16 deletions(-) create mode 100644 Grid/qcd/utils/BaryonUtils.h create mode 100644 tests/qdpxx/Test_qdpxx_baryon.cc diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h new file mode 100644 index 00000000..4f65ca67 --- /dev/null +++ b/Grid/qcd/utils/BaryonUtils.h @@ -0,0 +1,176 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./lib/qcd/utils/BaryonUtils.h + + Copyright (C) 2019 + + Author: Felix Erben + Author: Michael Marshall + + 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 */ +#pragma once +//#include +#include + +NAMESPACE_BEGIN(Grid); + +#undef DELTA_F_EQ_2 + +template +class BaryonUtils +{ +public: + typedef typename FImpl::ComplexField ComplexField; + typedef typename FImpl::FermionField FermionField; + typedef typename FImpl::PropagatorField PropagatorField; + + typedef typename FImpl::SitePropagator pobj; + typedef typename FImpl::SiteSpinor vobj; + typedef typename vobj::scalar_object sobj; + typedef typename vobj::scalar_type scalar_type; + typedef typename vobj::vector_type vector_type; + + static void ContractBaryons(const PropagatorField &q1_src, + const PropagatorField &q2_src, + const PropagatorField &q3_src, + const Gamma GammaA, + const Gamma GammaB, + const char quarks_snk[], + const char quarks_src[], + const int parity, + ComplexField &baryon_corr); + +}; + + +template +void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, + const PropagatorField &q2_src, + const PropagatorField &q3_src, + const Gamma GammaA, + const Gamma GammaB, + const char quarks_snk[], + const char quarks_src[], + const int parity, + ComplexField &baryon_corr) +{ + + assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); + + GridBase *grid = q1_src.Grid(); + + 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}; + std::vector wick_contraction = {0,0,0,0,0,0}; + + for (int ie=0; ie < 6 ; ie++) + if (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) + wick_contraction[ie]=1; + + typedef typename ComplexField::vector_object vobj; + LatticeView vbaryon_corr{ baryon_corr }; + accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + + LatticeView v1(q1_src); + LatticeView v2(q2_src); + LatticeView v3(q3_src); + auto D1 = v1[ss]; + auto D2 = v2[ss]; + auto D3 = v3[ss]; + + auto gD1a = GammaA * GammaA * D1; + auto gD1b = GammaA * g4 * GammaA * D1; + auto pD1 = 0.5* (gD1a + (double)parity * gD1b); + auto gD3 = GammaB * D3; + + vobj result{ 0 }; + + 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' + //This is the \delta_{456}^{123} part + if (wick_contraction[0]){ + auto D2g = D2 * GammaB; + for (int alpha_snk=0; alpha_snk -Author: Lanny91 +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 @@ -33,6 +33,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE @@ -45,9 +46,16 @@ class BaryonPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(BaryonPar, - std::string, q1, - std::string, q2, - std::string, q3, + std::string, q1_src, + std::string, q2_src, + std::string, q3_src, + std::string, GammaA, + std::string, GammaB, + // char[], quarks_snk, + // char[], quarks_src, + std::string, quarks_snk, + std::string, quarks_src, + int, parity, std::string, output); }; @@ -62,7 +70,7 @@ public: { public: GRID_SERIALIZABLE_CLASS_MEMBERS(Result, - std::vector>>, corr); + std::vector, corr); }; public: // constructor @@ -77,6 +85,8 @@ protected: virtual void setup(void); // execution virtual void execute(void); + // Which gamma algebra was specified + Gamma::Algebra al; }; MODULE_REGISTER_TMP(Baryon, ARG(TBaryon), MContraction); @@ -94,7 +104,7 @@ TBaryon::TBaryon(const std::string name) template std::vector TBaryon::getInput(void) { - std::vector input = {par().q1, par().q2, par().q3}; + std::vector input = {par().q1_src, par().q2_src, par().q3_src}; return input; } @@ -118,19 +128,38 @@ void TBaryon::setup(void) template void TBaryon::execute(void) { - LOG(Message) << "Computing baryon contractions '" << getName() << "' using" - << " quarks '" << par().q1 << "', '" << par().q2 << "', and '" - << par().q3 << "'" << std::endl; + LOG(Message) << "Computing baryon contraction '" << getName() << "' < " << par().quarks_snk << " | " << par().quarks_src << " > using" + << " quarks '" << par().q1_src << "', and a diquark formed of ('" << par().q2_src << "', and '" + << par().q3_src << "') at the source and (Gamma^A,Gamma^B) = ( " << par().GammaA << " , " << par().GammaB + << " ) and parity " << par().parity << "." << std::endl; - auto &q1 = envGet(PropagatorField1, par().q1); - auto &q2 = envGet(PropagatorField2, par().q2); - auto &q3 = envGet(PropagatorField3, par().q2); + auto &q1_src = envGet(PropagatorField1, par().q1_src); + auto &q2_src = envGet(PropagatorField2, par().q2_src); + auto &q3_src = envGet(PropagatorField3, par().q3_src); envGetTmp(LatticeComplex, c); Result result; - - // FIXME: do contractions - - // saveResult(par().output, "meson", result); + int nt = env().getDim(Tp); + result.corr.resize(nt); + std::vector ggA = strToVec(par().GammaA); + Gamma GammaA = ggA[0]; + std::vector ggB = strToVec(par().GammaB); + Gamma GammaB = ggB[0]; + std::vector buf; + const int parity {par().parity}; + const char quarks_snk[]{par().quarks_snk.c_str()}; + const char quarks_src[]{par().quarks_src.c_str()}; + + BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); + + 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 diff --git a/tests/qdpxx/Test_qdpxx_baryon.cc b/tests/qdpxx/Test_qdpxx_baryon.cc new file mode 100644 index 00000000..9b5fbee8 --- /dev/null +++ b/tests/qdpxx/Test_qdpxx_baryon.cc @@ -0,0 +1,605 @@ +/************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/qdpxx/Test_qdpxx_wilson.cc + + Copyright (C) 2017 + + Author: Guido Cossu + + 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 +#include +#include + +typedef Grid::LatticeGaugeField GaugeField; + +namespace Chroma +{ + +class ChromaWrapper +{ +public: + typedef multi1d U; + typedef LatticeFermion T4; + + static void ImportGauge(GaugeField &gr, + QDP::multi1d &ch) + { + Grid::LorentzColourMatrix LCM; + Grid::Complex cc; + QDP::ColorMatrix cm; + QDP::Complex c; + + std::vector x(4); + QDP::multi1d cx(4); + Grid::Coordinate gd = gr.Grid()->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + Grid::peekSite(LCM, gr, x); + + for (int mu = 0; mu < 4; mu++) + { + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + cc = LCM(mu)()(i, j); + c = QDP::cmplx(QDP::Real(real(cc)), QDP::Real(imag(cc))); + QDP::pokeColor(cm, c, i, j); + } + } + QDP::pokeSite(ch[mu], cm, cx); + } + } + } + } + } + } + + static void ExportGauge(GaugeField &gr, + QDP::multi1d &ch) + { + Grid::LorentzColourMatrix LCM; + Grid::Complex cc; + QDP::ColorMatrix cm; + QDP::Complex c; + + std::vector x(4); + QDP::multi1d cx(4); + Grid::Coordinate gd = gr.Grid()->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + for (int mu = 0; mu < 4; mu++) + { + for (int i = 0; i < 3; i++) + { + for (int j = 0; j < 3; j++) + { + cm = QDP::peekSite(ch[mu], cx); + c = QDP::peekColor(cm, i, j); + cc = Grid::Complex(toDouble(real(c)), toDouble(imag(c))); + LCM(mu) + ()(i, j) = cc; + } + } + } + Grid::pokeSite(LCM, gr, x); + } + } + } + } + } + + // Specific for Wilson Fermions + static void ImportPropagator(Grid::LatticePropagator &gr, + QDP::LatticePropagator &ch) + { + Grid::LatticeSpinColourVector LF(gr.Grid()); + QDP::LatticeFermion cLF; + + int Nspin=4; + int Ncolour=3; + + for (int is = 0; is < Nspin; is++){ + for (int ic = 0; ic < Ncolour; ic++){ + Grid::PropToFerm(LF,gr,is,ic); + ImportFermion(LF,cLF); + Chroma::FermToProp(cLF,ch,ic,is); + } + } + } + + static void ExportPropagator(Grid::LatticePropagator &gr, + QDP::LatticePropagator &ch) + { + Grid::LatticeSpinColourVector LF(gr.Grid()); + QDP::LatticeFermion cLF; + + int Nspin=4; + int Ncolour=3; + + for (int is = 0; is < Nspin; is++){ + for (int ic = 0; ic < Ncolour; ic++){ + Chroma::PropToFerm(ch,cLF,ic,is); + ExportFermion(LF,cLF); + Grid::FermToProp(gr,LF,is,ic); + } + } + } + + // Specific for Wilson Fermions + static void ImportFermion(Grid::LatticeFermion &gr, + QDP::LatticeFermion &ch) + { + Grid::SpinColourVector F; + Grid::Complex c; + + QDP::Fermion cF; + QDP::SpinVector cS; + QDP::Complex cc; + + std::vector x(4); // explicit 4d fermions in Grid + QDP::multi1d cx(4); + Grid::Coordinate gd = gr.Grid()->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + Grid::peekSite(F, gr, x); + + for (int j = 0; j < 3; j++) + { + for (int sp = 0; sp < 4; sp++) + { + + c = F()(sp)(j); + + cc = QDP::cmplx(QDP::Real(real(c)), QDP::Real(imag(c))); + + QDP::pokeSpin(cS, cc, sp); + } + QDP::pokeColor(cF, cS, j); + } + QDP::pokeSite(ch, cF, cx); + } + } + } + } + } + + // Specific for 4d Wilson fermions + static void ExportFermion(Grid::LatticeFermion &gr, + QDP::LatticeFermion &ch) + { + Grid::SpinColourVector F; + Grid::Complex c; + + QDP::Fermion cF; + QDP::SpinVector cS; + QDP::Complex cc; + + std::vector x(4); // 4d fermions + QDP::multi1d cx(4); + Grid::Coordinate gd = gr.Grid()->GlobalDimensions(); + + for (x[0] = 0; x[0] < gd[0]; x[0]++) + { + for (x[1] = 0; x[1] < gd[1]; x[1]++) + { + for (x[2] = 0; x[2] < gd[2]; x[2]++) + { + for (x[3] = 0; x[3] < gd[3]; x[3]++) + { + cx[0] = x[0]; + cx[1] = x[1]; + cx[2] = x[2]; + cx[3] = x[3]; + + cF = QDP::peekSite(ch, cx); + for (int sp = 0; sp < 4; sp++) + { + for (int j = 0; j < 3; j++) + { + cS = QDP::peekColor(cF, j); + cc = QDP::peekSpin(cS, sp); + c = Grid::Complex(QDP::toDouble(QDP::real(cc)), + QDP::toDouble(QDP::imag(cc))); + F() + (sp)(j) = c; + } + } + Grid::pokeSite(F, gr, x); + } + } + } + } + } + +}; +} // namespace Chroma + +void make_gauge(GaugeField &Umu, Grid::LatticePropagator &q1,Grid::LatticePropagator &q2,Grid::LatticePropagator &q3) +{ + using namespace Grid; + using namespace Grid::QCD; + + std::vector seeds4({1, 2, 3, 4}); + + Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu.Grid(); + Grid::GridParallelRNG RNG4(UGrid); + RNG4.SeedFixedIntegers(seeds4); + Grid::SU3::HotConfiguration(RNG4, Umu); + + // Propagator + Grid::gaussian(RNG4, q1); + Grid::gaussian(RNG4, q2); + Grid::gaussian(RNG4, q3); +} + +void calc_chroma(GaugeField &lat, Grid::LatticePropagator &qU,Grid::LatticePropagator &qD,Grid::LatticePropagator &qS, std::vector &res, std::string baryon) +{ + QDP::multi1d u(4); + Chroma::ChromaWrapper::ImportGauge(lat, u); + + QDP::LatticePropagator check; + QDP::LatticePropagator result; + QDP::LatticePropagator psiU; + QDP::LatticePropagator psiD; + QDP::LatticePropagator psiS; + + + Chroma::ChromaWrapper::ImportPropagator(qU, psiU); + Chroma::ChromaWrapper::ImportPropagator(qD, psiD); + Chroma::ChromaWrapper::ImportPropagator(qS, psiS); + + if(0){ + std::cout << "Testing ImportPropagator(): " << std::endl; + Grid::GridCartesian *UGrid = (Grid::GridCartesian *)lat.Grid(); + std::vector buf; + Grid::LatticeComplex tmp(UGrid); + tmp = Grid::trace(qU); + Grid::sliceSum(tmp,buf,Grid::Nd-1); + for (unsigned int t = 0; t < buf.size(); ++t) + { + std::cout << "Grid qU " << t << " " << Grid::TensorRemove(buf[t]) << std::endl; + } + QDP::LatticeComplex ctmp; + ctmp = QDP::trace(psiU); + Chroma::SftMom phases0(0,true,3); //How do I circumvent this? sliceSum equivalent? + QDP::multi2d hsum0; + hsum0 = phases0.sft(ctmp); + for(int t = 0; t < phases0.numSubsets(); ++t){ + std::cout << "Chroma qU " << t << " " << hsum0[0][t] << std::endl; + } + } + + SpinMatrix C; + SpinMatrix C_5; + SpinMatrix C_4_5; + SpinMatrix CG_1; + SpinMatrix CG_2; + SpinMatrix CG_3; + SpinMatrix CG_4; + + + SpinMatrix g_one = 1.0; + //C = \gamma_2\gamma_4 + C = (Gamma(10)*g_one); + + //C_5 = C*gamma_5 + C_5 = (Gamma(5)*g_one); + + //C_4_5 = C*gamma_4*gamma_5 + C_4_5 = (Gamma(13)*g_one); + + //CG_1 = C*gamma_1 + CG_1 = (Gamma(11)*g_one); + + //CG_2 = C*gamma_2 + CG_2 = (Gamma(8)*g_one); + + //CG_3 = C*gamma_3 + CG_3 = (Gamma(14)*g_one); + + //CG_4 = C*gamma_4 + CG_4 = (Gamma(2)*g_one); + + // S_proj_unpol = (1/2)(1 + gamma_4) + SpinMatrix S_proj_unpol = 0.5 * (g_one + (g_one * Gamma(8))); + + QDP::LatticeComplex b_prop; + QDP::LatticePropagator di_quark; + + if(! baryon.compare("OmegaX")){ + // Omega_x - this esentially is degenerate (s C\gamma_1 s)s + // C gamma_1 = Gamma(10) * Gamma(1) = Gamma(11) + di_quark = QDP::quarkContract13(psiS * CG_1, CG_1 * psiS); + b_prop = QDP::trace(S_proj_unpol * QDP::traceColor(psiS * QDP::traceSpin(di_quark))) + + 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * di_quark)); + } else if (! baryon.compare("OmegaY")){ + // Omega_x - this esentially is degenerate (s C\gamma_3 s)s + // C gamma_1 = Gamma(10) * Gamma(2) = Gamma(8) + di_quark = QDP::quarkContract13(psiS * CG_2, CG_2 * psiS); + b_prop = QDP::trace(S_proj_unpol * QDP::traceColor(psiS * QDP::traceSpin(di_quark))) + + 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * di_quark)); + } else if (! baryon.compare("OmegaZ")){ + // Omega_x - this esentially is degenerate (s C\gamma_3 s)s + // C gamma_1 = Gamma(10) * Gamma(4) = Gamma(14) + di_quark = QDP::quarkContract13(psiS * CG_3, CG_3 * psiS); + b_prop = QDP::trace(S_proj_unpol * QDP::traceColor(psiS * QDP::traceSpin(di_quark))) + + 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * di_quark)); + } else if (! baryon.compare("Proton")){ + // Proton - this esentially is degenerate (d C\gamma_5 u)u +// This is how the UKHadron code is written - diquarks are swapped when compared to coment above code. + //di_quark = QDP::quarkContract13(psiU * C_5, C_5 * psiD); + di_quark = QDP::quarkContract13(psiD * C_5, C_5 * psiU); + b_prop = QDP::trace(S_proj_unpol * QDP::traceColor(psiU * QDP::traceSpin(di_quark))) + + QDP::trace(S_proj_unpol * QDP::traceColor(psiU * di_quark)); + } else if (! baryon.compare("Lambda")){ + // Lambda (octet) - This is the totally antisymmetric + // one from the middle of the octet + // Lambda - (d C\gamma_5 s)u - (u C\gamma_5 s)d + // This is given by: + // 1/3[ d + u + 4s - (usd) - (dsu) + 2(sud) + 2(sdu) + 2(uds) + 2(dus) ] + +/* This is how the UKHadron code is written - diquarks are swapped when compared to coments above code. + // This gives d - (usd) -- yes + di_quark = QDP::quarkContract13(psiU * C_5, C_5 * psiS); + b_prop = QDP::trace(S_proj_unpol * QDP::traceColor(psiD * QDP::traceSpin(di_quark))) + - QDP::trace(S_proj_unpol * QDP::traceColor(psiD * di_quark)); + + // This gives u - (dsu) -- yes + di_quark = quarkContract13(psiD * C_5,C_5 * psiS); + b_prop += QDP::trace(S_proj_unpol * QDP::traceColor(psiU * QDP::traceSpin(di_quark))) + - QDP::trace(S_proj_unpol * QDP::traceColor(psiU * di_quark)); + + // This gives 4s -- yes + di_quark = quarkContract13(psiU * C_5,C_5 * psiD); + b_prop += 4.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * QDP::traceSpin(di_quark))); + + //This gives 2(sud) -- yes + di_quark = quarkContract13(psiS * C_5,C_5 * psiU); + b_prop += 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiD * di_quark)); + + // This gives 2(sdu) -- yes + di_quark = quarkContract13(psiS * C_5,C_5 * psiD); + b_prop += 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiU * di_quark)); + + // This gives 2(uds) -- yes + di_quark = quarkContract13(psiU * C_5,C_5 * psiD); + b_prop += 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * di_quark)); + + // This gives 2(dus) -- yes + di_quark = quarkContract13(psiD * C_5,C_5 * psiU); + b_prop += 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * di_quark));*/ + + // This gives d - (usd) -- yes + di_quark = QDP::quarkContract13(psiS * C_5, C_5 * psiU); + b_prop = QDP::trace(S_proj_unpol * QDP::traceColor(psiD * QDP::traceSpin(di_quark))) + - QDP::trace(S_proj_unpol * QDP::traceColor(psiD * di_quark)); + + // This gives u - (dsu) -- yes + di_quark = quarkContract13(psiS * C_5,C_5 * psiD); + b_prop += QDP::trace(S_proj_unpol * QDP::traceColor(psiU * QDP::traceSpin(di_quark))) + - QDP::trace(S_proj_unpol * QDP::traceColor(psiU * di_quark)); + + // This gives 4s -- yes + di_quark = quarkContract13(psiD * C_5,C_5 * psiU); + b_prop += 4.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * QDP::traceSpin(di_quark))); + + //This gives 2(sud) -- yes + di_quark = quarkContract13(psiU * C_5,C_5 * psiS); + b_prop += 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiD * di_quark)); + + // This gives 2(sdu) -- yes + di_quark = quarkContract13(psiD * C_5,C_5 * psiS); + b_prop += 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiU * di_quark)); + + // This gives 2(uds) -- yes + di_quark = quarkContract13(psiD * C_5,C_5 * psiU); + b_prop += 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * di_quark)); + + // This gives 2(dus) -- yes + di_quark = quarkContract13(psiU * C_5,C_5 * psiD); + b_prop += 2.0 * QDP::trace(S_proj_unpol * QDP::traceColor(psiS * di_quark)); + } else { + std::cout << "baryon not part of test " << std::endl; + return; + } + std::cout<< "Chroma computing " << baryon << std::endl; + + Chroma::SftMom phases(0,true,3); //How do I circumvent this? sliceSum equivalent? + QDP::multi2d hsum; + hsum = phases.sft(b_prop); + int length = phases.numSubsets(); + res.resize(length); + for(int t = 0; t < length; ++t){ + res[t] = hsum[0][t]; //Should I test momentum? + } + +} + + +void calc_grid(Grid::LatticeGaugeField &Umu, Grid::LatticePropagator &qU, Grid::LatticePropagator &qD, Grid::LatticePropagator &qS, std::vector &res, std::string baryon) +{ + using namespace Grid; + using namespace Grid::QCD; + + Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu.Grid(); + + + Grid::Gamma G_A = Grid::Gamma(Grid::Gamma::Algebra::Identity); + Grid::Gamma G_B = Grid::Gamma(Grid::Gamma::Algebra::GammaZGamma5); // OmegaX: C*GammaX = i* GammaZ*Gamma5 + char quarks[] = "sss"; + + Grid::LatticeComplex c(UGrid); + Grid::LatticeComplex c1(UGrid); + + if(! baryon.compare("OmegaX")){ + BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,quarks,quarks,1,c); + c*=0.5; + std::cout << "Grid-Omega factor 2 larger than Chroma-Omega!!!" << std::endl; + } else if (! baryon.compare("OmegaY")){ + G_B = Grid::Gamma(Grid::Gamma::Algebra::GammaT); + BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,quarks,quarks,1,c); + c*=0.5; + std::cout << "Grid-Omega factor 2 larger than Chroma-Omega!!!" << std::endl; + } else if (! baryon.compare("OmegaZ")){ + G_B = Grid::Gamma(Grid::Gamma::Algebra::GammaXGamma5); + BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,quarks,quarks,1,c); + c*=0.5; + std::cout << "Grid-Omega factor 2 larger than Chroma-Omega!!!" << std::endl; + } else if (! baryon.compare("Proton")){ + G_B = Grid::Gamma(Grid::Gamma::Algebra::SigmaXZ); + BaryonUtils::ContractBaryons(qU,qD,qU,G_A,G_B,"udu","udu",1,c); + std::cout << "UKHadron-Proton has flipped diquarks in original code." << std::endl; + } else if (! baryon.compare("Lambda")){ + G_B = Grid::Gamma(Grid::Gamma::Algebra::SigmaXZ); + BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,"sud","sud",1,c1); //s + c = 4.*c1; + BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,"dus","dus",1,c1); //d + c += c1; + BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,"uds","uds",1,c1); //u + c += c1; + BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,"sud","dus",1,c1); //(sud) + c += 2.*c1; + BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,"sud","uds",1,c1); //(sdu) + c -= 2.*c1; + BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,"dus","sud",1,c1); //(dus) + c += 2.*c1; + BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,"dus","uds",1,c1); //-(dsu) + c -= c1; + BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,"uds","sud",1,c1); //(uds) + c -= 2.*c1; + BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,"uds","dus",1,c1); //-(usd) + c -= c1; + std::cout << "UKHadron-Lambda has flipped diquarks in original code." << std::endl; + } else { + std::cout << "baryon not part of test " << std::endl; + return; + } + std::cout<< "Grid computing " << baryon << std::endl; + + std::vector buf; + Grid::sliceSum(c,buf,Grid::Nd-1); + res.resize(buf.size()); + for (unsigned int t = 0; t < buf.size(); ++t) + { + res[t]=Grid::TensorRemove(buf[t]); + } + +} + +int main(int argc, char **argv) +{ + + /******************************************************** + * Setup QDP + *********************************************************/ + Chroma::initialize(&argc, &argv); + Chroma::WilsonTypeFermActs4DEnv::registerAll(); + + /******************************************************** + * Setup Grid + *********************************************************/ + Grid::Grid_init(&argc, &argv); + Grid::GridCartesian *UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(), + Grid::GridDefaultSimd(Grid::Nd, Grid::vComplex::Nsimd()), + Grid::GridDefaultMpi()); + + Grid::Coordinate gd = UGrid->GlobalDimensions(); + QDP::multi1d nrow(QDP::Nd); + for (int mu = 0; mu < 4; mu++) + nrow[mu] = gd[mu]; + + QDP::Layout::setLattSize(nrow); + QDP::Layout::create(); + + GaugeField Ug(UGrid); + typedef Grid::LatticePropagator PropagatorField; + PropagatorField up(UGrid); + PropagatorField down(UGrid); + PropagatorField strange(UGrid); + std::vector res_chroma; + std::vector res_grid; + std::vector difference; + + std::vector baryons({"OmegaX","OmegaY","OmegaZ","Proton","Lambda"}); + int nBaryon=baryons.size(); + + for (int iB = 0; iB < nBaryon; iB++) + { + make_gauge(Ug, up, down, strange); // fills the gauge field and the propagator with random numbers + + calc_chroma(Ug, up, down, strange, res_chroma,baryons[iB]); + + for(int t=0;t Date: Mon, 16 Sep 2019 17:55:58 +0100 Subject: [PATCH 02/18] bugfix --- Hadrons/Modules/MContraction/Baryon.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 84ea76eb..fc4cc8fc 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -146,8 +146,8 @@ void TBaryon::execute(void) Gamma GammaB = ggB[0]; std::vector buf; const int parity {par().parity}; - const char quarks_snk[]{par().quarks_snk.c_str()}; - const char quarks_src[]{par().quarks_src.c_str()}; + const char * quarks_snk{par().quarks_snk.c_str()}; + const char * quarks_src{par().quarks_src.c_str()}; BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); From 76c93aa44e66ccb1d3d6f1d4a4017c9b8d062aac Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 17 Sep 2019 14:36:26 +0100 Subject: [PATCH 03/18] weird bug... --- Grid/qcd/utils/BaryonUtils.h | 35 +++++++++++++++++-------- Hadrons/Modules/MContraction/Baryon.hpp | 32 +++++++++++++++++++--- 2 files changed, 53 insertions(+), 14 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 4f65ca67..e34c4d72 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -53,8 +53,8 @@ public: const PropagatorField &q3_src, const Gamma GammaA, const Gamma GammaB, - const char quarks_snk[], - const char quarks_src[], + const char * quarks_snk, + const char * quarks_src, const int parity, ComplexField &baryon_corr); @@ -67,12 +67,15 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, const PropagatorField &q3_src, const Gamma GammaA, const Gamma GammaB, - const char quarks_snk[], - const char quarks_src[], + const char * quarks_snk, + const char * quarks_src, const int parity, ComplexField &baryon_corr) { - + std::cout << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << std::endl; + std::cout << "GammaA " << (GammaA.g) << std::endl; + std::cout << "GammaB " << (GammaB.g) << std::endl; + assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); GridBase *grid = q1_src.Grid(); @@ -88,12 +91,16 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, wick_contraction[ie]=1; typedef typename ComplexField::vector_object vobj; - LatticeView vbaryon_corr{ baryon_corr }; - accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + auto vbaryon_corr= baryon_corr.View(); + auto v1 = q1_src.View(); + auto v2 = q2_src.View(); + auto v3 = q3_src.View(); + + std::cout << "wick contract " << wick_contraction << std::endl; + int count=0; + // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + thread_for(ss,grid->oSites(),{ - LatticeView v1(q1_src); - LatticeView v2(q2_src); - LatticeView v3(q3_src); auto D1 = v1[ss]; auto D2 = v2[ss]; auto D3 = v3[ss]; @@ -103,7 +110,11 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, auto pD1 = 0.5* (gD1a + (double)parity * gD1b); auto gD3 = GammaB * D3; - vobj result{ 0 }; + vobj result=Zero(); + if (count<10){ + std::cout << "outside epsilon " << count << std::endl; + count++; + } for (int ie_src=0; ie_src < 6 ; ie_src++){ int a_src = epsilon[ie_src][0]; //a @@ -132,6 +143,7 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, }}} } //This is the \delta_{456}^{312} part + // for(int test=0;test<3;test++){ if (wick_contraction[2]){ auto gD3g = gD3 * GammaB; for (int alpha_snk=0; alpha_snk::ContractBaryons(const PropagatorField &q1_src, result()()() -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1()(gamma_src,gamma_src)(c_snk,c_src)*D2()(alpha_snk,beta_src)(a_snk,b_src)*gD3g()(alpha_snk,beta_src)(b_snk,a_src); }}} } + // } //This is the \delta_{456}^{321} part if (wick_contraction[4]){ auto D2g = D2 * GammaB; diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index fc4cc8fc..8f831c9b 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -34,6 +34,7 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include +#include BEGIN_HADRONS_NAMESPACE @@ -134,24 +135,49 @@ void TBaryon::execute(void) << " ) and parity " << par().parity << "." << std::endl; auto &q1_src = envGet(PropagatorField1, par().q1_src); + LOG(Message) << "1" << std::endl; auto &q2_src = envGet(PropagatorField2, par().q2_src); auto &q3_src = envGet(PropagatorField3, par().q3_src); envGetTmp(LatticeComplex, c); Result result; int nt = env().getDim(Tp); result.corr.resize(nt); + LOG(Message) << "2" << std::endl; std::vector ggA = strToVec(par().GammaA); - Gamma GammaA = ggA[0]; + LOG(Message) << "3" << std::endl; + Gamma GammaA(ggA[0]); + LOG(Message) << "4" << std::endl; std::vector ggB = strToVec(par().GammaB); - Gamma GammaB = ggB[0]; + Gamma GammaB(ggB[0]); std::vector buf; const int parity {par().parity}; + LOG(Message) << "5" << std::endl; const char * quarks_snk{par().quarks_snk.c_str()}; + LOG(Message) << "6" << std::endl; const char * quarks_src{par().quarks_src.c_str()}; + LOG(Message) << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << std::endl; + LOG(Message) << "GammaA " << (GammaA.g) << std::endl; + LOG(Message) << "GammaB " << (GammaB.g) << std::endl; - BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); + GridCartesian *Ugrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); + LatticePropagator q(Ugrid); + GridParallelRNG RNG4(Ugrid); + gaussian(RNG4,q); + Gamma gA = Gamma(Gamma::Algebra::Identity); + Gamma gB = Gamma(Gamma::Algebra::Gamma5); + int p=1; + char * om = "sss"; + LatticeComplex c2(Ugrid); + + //BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); + BaryonUtils::ContractBaryons(q,q,q,gA,gB,om,om,p,c); + std::vector GA={gA}; + std::vector GB={gB}; + //A2Autils::ContractFourQuarkColourMix(q,q,GA,GB,c,c2); + LOG(Message) << "survived ContractBaryons" << std::endl; sliceSum(c,buf,Tp); + LOG(Message) << "survived sliceSum" << std::endl; for (unsigned int t = 0; t < buf.size(); ++t) { From 8415e23fc647cddceb7e616513c6fb019bf6b795 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 26 Sep 2019 11:09:09 +0100 Subject: [PATCH 04/18] still bugfix --- Grid/qcd/utils/BaryonUtils.h | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index e34c4d72..f34f2926 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -111,10 +111,10 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, auto gD3 = GammaB * D3; vobj result=Zero(); - if (count<10){ + /* if (count<10){ std::cout << "outside epsilon " << count << std::endl; count++; - } + }*/ for (int ie_src=0; ie_src < 6 ; ie_src++){ int a_src = epsilon[ie_src][0]; //a @@ -143,7 +143,6 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, }}} } //This is the \delta_{456}^{312} part - // for(int test=0;test<3;test++){ if (wick_contraction[2]){ auto gD3g = gD3 * GammaB; for (int alpha_snk=0; alpha_snk::ContractBaryons(const PropagatorField &q1_src, result()()() -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1()(gamma_src,gamma_src)(c_snk,c_src)*D2()(alpha_snk,beta_src)(a_snk,b_src)*gD3g()(alpha_snk,beta_src)(b_snk,a_src); }}} } - // } //This is the \delta_{456}^{321} part if (wick_contraction[4]){ auto D2g = D2 * GammaB; From bf62ec163d480f3c4a0508d3f6f6e7f0e8783221 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 26 Sep 2019 13:33:49 +0100 Subject: [PATCH 05/18] thread_for caused the problems - slow for loop for now --- Grid/qcd/utils/BaryonUtils.h | 11 ++----- Hadrons/Modules/MContraction/Baryon.hpp | 38 ++++++------------------- 2 files changed, 11 insertions(+), 38 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index f34f2926..2a7fc9d0 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -96,10 +96,9 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, auto v2 = q2_src.View(); auto v3 = q3_src.View(); - std::cout << "wick contract " << wick_contraction << std::endl; - int count=0; // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - thread_for(ss,grid->oSites(),{ + //thread_for(ss,grid->oSites(),{ + for(int ss=0; ss < grid->oSites(); ss++){ auto D1 = v1[ss]; auto D2 = v2[ss]; @@ -111,10 +110,6 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, auto gD3 = GammaB * D3; vobj result=Zero(); - /* if (count<10){ - std::cout << "outside epsilon " << count << std::endl; - count++; - }*/ for (int ie_src=0; ie_src < 6 ; ie_src++){ int a_src = epsilon[ie_src][0]; //a @@ -182,6 +177,6 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, } vbaryon_corr[ss] = result; - } );//end loop over lattice sites + } // );//end loop over lattice sites } NAMESPACE_END(Grid); diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 8f831c9b..4a3a39e4 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -52,11 +52,10 @@ public: std::string, q3_src, std::string, GammaA, std::string, GammaB, - // char[], quarks_snk, - // char[], quarks_src, std::string, quarks_snk, std::string, quarks_src, int, parity, + std::string, sink, std::string, output); }; @@ -67,6 +66,8 @@ public: FERM_TYPE_ALIASES(FImpl1, 1); FERM_TYPE_ALIASES(FImpl2, 2); FERM_TYPE_ALIASES(FImpl3, 3); + BASIC_TYPE_ALIASES(ScalarImplCR, Scalar); + SINK_TYPE_ALIASES(Scalar); class Result: Serializable { public: @@ -105,7 +106,7 @@ TBaryon::TBaryon(const std::string name) template std::vector TBaryon::getInput(void) { - std::vector input = {par().q1_src, par().q2_src, par().q3_src}; + std::vector input = {par().q1_src, par().q2_src, par().q3_src, par().sink}; return input; } @@ -135,49 +136,26 @@ void TBaryon::execute(void) << " ) and parity " << par().parity << "." << std::endl; auto &q1_src = envGet(PropagatorField1, par().q1_src); - LOG(Message) << "1" << std::endl; auto &q2_src = envGet(PropagatorField2, par().q2_src); auto &q3_src = envGet(PropagatorField3, par().q3_src); envGetTmp(LatticeComplex, c); Result result; int nt = env().getDim(Tp); result.corr.resize(nt); - LOG(Message) << "2" << std::endl; std::vector ggA = strToVec(par().GammaA); - LOG(Message) << "3" << std::endl; Gamma GammaA(ggA[0]); - LOG(Message) << "4" << std::endl; std::vector ggB = strToVec(par().GammaB); Gamma GammaB(ggB[0]); std::vector buf; const int parity {par().parity}; - LOG(Message) << "5" << std::endl; const char * quarks_snk{par().quarks_snk.c_str()}; - LOG(Message) << "6" << std::endl; const char * quarks_src{par().quarks_src.c_str()}; - LOG(Message) << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << std::endl; - LOG(Message) << "GammaA " << (GammaA.g) << std::endl; - LOG(Message) << "GammaB " << (GammaB.g) << std::endl; - GridCartesian *Ugrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); - LatticePropagator q(Ugrid); - GridParallelRNG RNG4(Ugrid); - gaussian(RNG4,q); - Gamma gA = Gamma(Gamma::Algebra::Identity); - Gamma gB = Gamma(Gamma::Algebra::Gamma5); - int p=1; - char * om = "sss"; - LatticeComplex c2(Ugrid); + BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); - //BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); - BaryonUtils::ContractBaryons(q,q,q,gA,gB,om,om,p,c); - std::vector GA={gA}; - std::vector GB={gB}; - //A2Autils::ContractFourQuarkColourMix(q,q,GA,GB,c,c2); - LOG(Message) << "survived ContractBaryons" << std::endl; - - sliceSum(c,buf,Tp); - LOG(Message) << "survived sliceSum" << std::endl; + //sliceSum(c,buf,Tp); + SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); + buf = sink(c); for (unsigned int t = 0; t < buf.size(); ++t) { From 94b9a9474c0c3bc85f77536e75857c9ba99a04d6 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Fri, 27 Sep 2019 15:08:56 +0100 Subject: [PATCH 06/18] Baryons module works in 1 of 3 cases - still need SlicedProp and Msource part!! --- Grid/qcd/utils/BaryonUtils.h | 164 +++++++++++++++++------- Hadrons/Modules/MContraction/Baryon.hpp | 56 ++++++-- 2 files changed, 164 insertions(+), 56 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 2a7fc9d0..5b7404c4 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -43,74 +43,67 @@ public: typedef typename FImpl::PropagatorField PropagatorField; typedef typename FImpl::SitePropagator pobj; - typedef typename FImpl::SiteSpinor vobj; + typedef typename ComplexField::vector_object vobj; + /* typedef typename FImpl::SiteSpinor vobj; typedef typename vobj::scalar_object sobj; typedef typename vobj::scalar_type scalar_type; typedef typename vobj::vector_type vector_type; - +*/ + private: + static void baryon_site(const pobj &D1, + const pobj &D2, + const pobj &D3, + const Gamma GammaA, + const Gamma GammaB, + const int parity, + const std::vector wick_contractions, + vobj &result); + public: static void ContractBaryons(const PropagatorField &q1_src, - const PropagatorField &q2_src, - const PropagatorField &q3_src, - const Gamma GammaA, - const Gamma GammaB, - const char * quarks_snk, - const char * quarks_src, - const int parity, - ComplexField &baryon_corr); - + const PropagatorField &q2_src, + const PropagatorField &q3_src, + const Gamma GammaA, + const Gamma GammaB, + const char * quarks_snk, + const char * quarks_src, + const int parity, + ComplexField &baryon_corr); + /* template + static void ContractBaryons_Sliced(const T1 &D1, + const T2 &D2, + const T3 &D3,*/ + static void ContractBaryons_Sliced(const pobj &D1, + const pobj &D2, + const pobj &D3, + const Gamma GammaA, + const Gamma GammaB, + const char * quarks_snk, + const char * quarks_src, + const int parity, + vobj &result); }; - template -void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, - const PropagatorField &q2_src, - const PropagatorField &q3_src, +void BaryonUtils::baryon_site(const pobj &D1, + const pobj &D2, + const pobj &D3, const Gamma GammaA, const Gamma GammaB, - const char * quarks_snk, - const char * quarks_src, const int parity, - ComplexField &baryon_corr) + const std::vector wick_contraction, + vobj &result) { - std::cout << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << std::endl; - std::cout << "GammaA " << (GammaA.g) << std::endl; - std::cout << "GammaB " << (GammaB.g) << std::endl; - - assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); - - GridBase *grid = q1_src.Grid(); 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}; - std::vector wick_contraction = {0,0,0,0,0,0}; - - for (int ie=0; ie < 6 ; ie++) - if (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) - wick_contraction[ie]=1; - - typedef typename ComplexField::vector_object vobj; - auto vbaryon_corr= baryon_corr.View(); - auto v1 = q1_src.View(); - auto v2 = q2_src.View(); - auto v3 = q3_src.View(); - - // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - //thread_for(ss,grid->oSites(),{ - for(int ss=0; ss < grid->oSites(); ss++){ - - auto D1 = v1[ss]; - auto D2 = v2[ss]; - auto D3 = v3[ss]; auto gD1a = GammaA * GammaA * D1; auto gD1b = GammaA * g4 * GammaA * D1; auto pD1 = 0.5* (gD1a + (double)parity * gD1b); auto gD3 = GammaB * D3; - vobj result=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 @@ -175,8 +168,85 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, } } } - vbaryon_corr[ss] = result; +} +template +void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, + const PropagatorField &q2_src, + const PropagatorField &q3_src, + const Gamma GammaA, + const Gamma GammaB, + const char * quarks_snk, + const char * quarks_src, + const int parity, + ComplexField &baryon_corr) +{ + std::cout << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << std::endl; + std::cout << "GammaA " << (GammaA.g) << std::endl; + std::cout << "GammaB " << (GammaB.g) << std::endl; + + assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); + + GridBase *grid = q1_src.Grid(); + + 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}; + std::vector wick_contraction = {0,0,0,0,0,0}; + + for (int ie=0; ie < 6 ; ie++) + if (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) + wick_contraction[ie]=1; + +// typedef typename ComplexField::vector_object vobj; + auto vbaryon_corr= baryon_corr.View(); + auto v1 = q1_src.View(); + auto v2 = q2_src.View(); + auto v3 = q3_src.View(); + + // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + //thread_for(ss,grid->oSites(),{ + for(int ss=0; ss < grid->oSites(); ss++){ + + auto D1 = v1[ss]; + auto D2 = v2[ss]; + auto D3 = v3[ss]; + + vobj result=Zero(); + baryon_site(D1,D2,D3,GammaA,GammaB,parity,wick_contraction,result); + vbaryon_corr[ss] = result; } // );//end loop over lattice sites } +/*template +void BaryonUtils::ContractBaryons_Sliced(const T1 &D1, + const T2 &D2, + const T3 &D3,*/ +template +void BaryonUtils::ContractBaryons_Sliced(const pobj &D1, + const pobj &D2, + const pobj &D3, + const Gamma GammaA, + const Gamma GammaB, + const char * quarks_snk, + const char * quarks_src, + const int parity, + vobj &result) +{ + + assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); + + 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}; + std::vector wick_contraction = {0,0,0,0,0,0}; + + for (int ie=0; ie < 6 ; ie++) + if (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) + wick_contraction[ie]=1; + + result=Zero(); + baryon_site(D1,D2,D3,GammaA,GammaB,parity,wick_contraction,result); +} NAMESPACE_END(Grid); diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 4a3a39e4..51751f0a 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -34,7 +34,6 @@ See the full license in the file "LICENSE" in the top level distribution directo #include #include #include -#include BEGIN_HADRONS_NAMESPACE @@ -147,19 +146,58 @@ void TBaryon::execute(void) std::vector ggB = strToVec(par().GammaB); Gamma GammaB(ggB[0]); std::vector buf; + vTComplex cs; const int parity {par().parity}; const char * quarks_snk{par().quarks_snk.c_str()}; const char * quarks_src{par().quarks_src.c_str()}; - BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); - - //sliceSum(c,buf,Tp); - SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); - buf = sink(c); - - for (unsigned int t = 0; t < buf.size(); ++t) + if (envHasType(SlicedPropagator1, par().q1_src) and + envHasType(SlicedPropagator2, par().q2_src) and + envHasType(SlicedPropagator3, par().q3_src)) { - result.corr[t] = TensorRemove(buf[t]); + auto &q1_src = envGet(SlicedPropagator1, par().q1_src); + auto &q2_src = envGet(SlicedPropagator2, par().q2_src); + auto &q3_src = envGet(SlicedPropagator3, par().q3_src); + + LOG(Message) << "(propagator already sinked)" << std::endl; + for (unsigned int t = 0; t < buf.size(); ++t) + { + //TODO: Get this to compile without the casts. Templates? + //BaryonUtils::ContractBaryons_Sliced(*reinterpret_cast, 4>>*>(&q1_src[t]),*reinterpret_cast, 4>>*>(&q2_src[t]),*reinterpret_cast, 4>>*>(&q3_src[t]),GammaA,GammaB,quarks_snk,quarks_src,parity,cs); + //result.corr[t] = TensorRemove(*reinterpret_cast(&cs)); + // BaryonUtils::ContractBaryons_Sliced(q1_src[t],q2_src[t],q3_src[t],GammaA,GammaB,quarks_snk,quarks_src,parity,cs); + // result.corr[t] = TensorRemove(cs); + } + } + else + { + std::string ns; + + ns = vm().getModuleNamespace(env().getObjectModule(par().sink)); + if (ns == "MSource") + { + //TODO: Understand what this is and then get it to compile. Hopefully no new function needed. The following lines are from the Meson.hpp module. + /* PropagatorField1 &sink = envGet(PropagatorField1, par().sink); + + c = trace(mesonConnected(q1, q2, gSnk, gSrc)*sink); + sliceSum(c, buf, Tp); */ +// My attempt at some code, which doesn't work. I also don't know whether anything like this is what we want here. + /* BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); + PropagatorField1 &sink = envGet(PropagatorField1, par().sink); + auto test = trace(c*sink); + sliceSum(test, buf, Tp); */ + } + else if (ns == "MSink") + { + BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); + + SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); + buf = sink(c); + } + for (unsigned int t = 0; t < buf.size(); ++t) + { + result.corr[t] = TensorRemove(buf[t]); + } } saveResult(par().output, "baryon", result); From e5d7910fa77316cf4e82caf264b02871c62ac513 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Mon, 30 Sep 2019 13:55:26 +0100 Subject: [PATCH 07/18] code compiling now - still need to test --- Grid/qcd/utils/BaryonUtils.h | 67 ++++++++++++------------- Hadrons/Modules/MContraction/Baryon.hpp | 22 +++----- 2 files changed, 38 insertions(+), 51 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 5b7404c4..6c00a414 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -43,21 +43,18 @@ public: typedef typename FImpl::PropagatorField PropagatorField; typedef typename FImpl::SitePropagator pobj; - typedef typename ComplexField::vector_object vobj; - /* typedef typename FImpl::SiteSpinor vobj; - typedef typename vobj::scalar_object sobj; - typedef typename vobj::scalar_type scalar_type; - typedef typename vobj::vector_type vector_type; -*/ + typedef typename ComplexField::vector_object vobj; + private: - static void baryon_site(const pobj &D1, - const pobj &D2, - const pobj &D3, + template + static void baryon_site(const mobj &D1, + const mobj &D2, + const mobj &D3, const Gamma GammaA, const Gamma GammaB, const int parity, - const std::vector wick_contractions, - vobj &result); + const std::vector &wick_contractions, + robj &result); public: static void ContractBaryons(const PropagatorField &q1_src, const PropagatorField &q2_src, @@ -68,36 +65,35 @@ public: const char * quarks_src, const int parity, ComplexField &baryon_corr); - /* template - static void ContractBaryons_Sliced(const T1 &D1, - const T2 &D2, - const T3 &D3,*/ - static void ContractBaryons_Sliced(const pobj &D1, - const pobj &D2, - const pobj &D3, + template + static void ContractBaryons_Sliced(const mobj &D1, + const mobj &D2, + const mobj &D3, const Gamma GammaA, const Gamma GammaB, const char * quarks_snk, const char * quarks_src, const int parity, - vobj &result); + robj &result); }; -template -void BaryonUtils::baryon_site(const pobj &D1, - const pobj &D2, - const pobj &D3, +template +template +void BaryonUtils::baryon_site(const mobj &D1, + const mobj &D2, + const mobj &D3, const Gamma GammaA, const Gamma GammaB, const int parity, - const std::vector wick_contraction, - vobj &result) + const std::vector &wick_contraction, + robj &result) { 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}; + + static const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + static const int epsilon_sgn[6]= {1,1,1,-1,-1,-1}; auto gD1a = GammaA * GammaA * D1; auto gD1b = GammaA * g4 * GammaA * D1; @@ -206,8 +202,8 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, auto v3 = q3_src.View(); // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { - //thread_for(ss,grid->oSites(),{ - for(int ss=0; ss < grid->oSites(); ss++){ + thread_for(ss,grid->oSites(),{ + //for(int ss=0; ss < grid->oSites(); ss++){ auto D1 = v1[ss]; auto D2 = v2[ss]; @@ -216,22 +212,23 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, vobj result=Zero(); baryon_site(D1,D2,D3,GammaA,GammaB,parity,wick_contraction,result); vbaryon_corr[ss] = result; - } // );//end loop over lattice sites + } );//end loop over lattice sites } /*template void BaryonUtils::ContractBaryons_Sliced(const T1 &D1, const T2 &D2, const T3 &D3,*/ -template -void BaryonUtils::ContractBaryons_Sliced(const pobj &D1, - const pobj &D2, - const pobj &D3, +template +template +void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, + const mobj &D2, + const mobj &D3, const Gamma GammaA, const Gamma GammaB, const char * quarks_snk, const char * quarks_src, const int parity, - vobj &result) + robj &result) { assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 51751f0a..a228bfeb 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -146,7 +146,7 @@ void TBaryon::execute(void) std::vector ggB = strToVec(par().GammaB); Gamma GammaB(ggB[0]); std::vector buf; - vTComplex cs; + TComplex cs; const int parity {par().parity}; const char * quarks_snk{par().quarks_snk.c_str()}; const char * quarks_src{par().quarks_src.c_str()}; @@ -162,11 +162,8 @@ void TBaryon::execute(void) LOG(Message) << "(propagator already sinked)" << std::endl; for (unsigned int t = 0; t < buf.size(); ++t) { - //TODO: Get this to compile without the casts. Templates? - //BaryonUtils::ContractBaryons_Sliced(*reinterpret_cast, 4>>*>(&q1_src[t]),*reinterpret_cast, 4>>*>(&q2_src[t]),*reinterpret_cast, 4>>*>(&q3_src[t]),GammaA,GammaB,quarks_snk,quarks_src,parity,cs); - //result.corr[t] = TensorRemove(*reinterpret_cast(&cs)); - // BaryonUtils::ContractBaryons_Sliced(q1_src[t],q2_src[t],q3_src[t],GammaA,GammaB,quarks_snk,quarks_src,parity,cs); - // result.corr[t] = TensorRemove(cs); + BaryonUtils::ContractBaryons_Sliced(q1_src[t],q2_src[t],q3_src[t],GammaA,GammaB,quarks_snk,quarks_src,parity,cs); + result.corr[t] = TensorRemove(cs); } } else @@ -176,21 +173,14 @@ void TBaryon::execute(void) ns = vm().getModuleNamespace(env().getObjectModule(par().sink)); if (ns == "MSource") { - //TODO: Understand what this is and then get it to compile. Hopefully no new function needed. The following lines are from the Meson.hpp module. - /* PropagatorField1 &sink = envGet(PropagatorField1, par().sink); - - c = trace(mesonConnected(q1, q2, gSnk, gSrc)*sink); - sliceSum(c, buf, Tp); */ -// My attempt at some code, which doesn't work. I also don't know whether anything like this is what we want here. - /* BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); + BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); PropagatorField1 &sink = envGet(PropagatorField1, par().sink); - auto test = trace(c*sink); - sliceSum(test, buf, Tp); */ + auto test = closure(trace(sink*c)); + sliceSum(test, buf, Tp); } else if (ns == "MSink") { BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); - SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); buf = sink(c); } From 155bcd4ff3636255f6c56e0b007ff4ef9207312d Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Mon, 30 Sep 2019 16:58:20 +0100 Subject: [PATCH 08/18] ready for pull request --- Grid/qcd/utils/BaryonUtils.h | 12 ++++-------- Hadrons/Modules/MContraction/Baryon.hpp | 6 +++--- 2 files changed, 7 insertions(+), 11 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 6c00a414..f049565b 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -187,8 +187,8 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, 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}; + static const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + static const int epsilon_sgn[6]= {1,1,1,-1,-1,-1}; std::vector wick_contraction = {0,0,0,0,0,0}; for (int ie=0; ie < 6 ; ie++) @@ -214,10 +214,6 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, vbaryon_corr[ss] = result; } );//end loop over lattice sites } -/*template -void BaryonUtils::ContractBaryons_Sliced(const T1 &D1, - const T2 &D2, - const T3 &D3,*/ template template void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, @@ -235,8 +231,8 @@ void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, 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}; + static const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + static const int epsilon_sgn[6]= {1,1,1,-1,-1,-1}; std::vector wick_contraction = {0,0,0,0,0,0}; for (int ie=0; ie < 6 ; ie++) diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index a228bfeb..6366f7d2 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -134,9 +134,6 @@ void TBaryon::execute(void) << par().q3_src << "') at the source and (Gamma^A,Gamma^B) = ( " << par().GammaA << " , " << par().GammaB << " ) and parity " << par().parity << "." << std::endl; - auto &q1_src = envGet(PropagatorField1, par().q1_src); - auto &q2_src = envGet(PropagatorField2, par().q2_src); - auto &q3_src = envGet(PropagatorField3, par().q3_src); envGetTmp(LatticeComplex, c); Result result; int nt = env().getDim(Tp); @@ -168,6 +165,9 @@ void TBaryon::execute(void) } else { + auto &q1_src = envGet(PropagatorField1, par().q1_src); + auto &q2_src = envGet(PropagatorField2, par().q2_src); + auto &q3_src = envGet(PropagatorField3, par().q3_src); std::string ns; ns = vm().getModuleNamespace(env().getObjectModule(par().sink)); From b88fd436e75f199a0c2d1abe101e282262404feb Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Mon, 30 Sep 2019 17:07:46 +0100 Subject: [PATCH 09/18] added author information --- Grid/qcd/utils/BaryonUtils.h | 1 - tests/qdpxx/Test_qdpxx_baryon.cc | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index f049565b..9e8d53b6 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -7,7 +7,6 @@ Copyright (C) 2019 Author: Felix Erben - Author: Michael Marshall 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 diff --git a/tests/qdpxx/Test_qdpxx_baryon.cc b/tests/qdpxx/Test_qdpxx_baryon.cc index 9b5fbee8..b19864fb 100644 --- a/tests/qdpxx/Test_qdpxx_baryon.cc +++ b/tests/qdpxx/Test_qdpxx_baryon.cc @@ -6,7 +6,7 @@ Copyright (C) 2017 - Author: Guido Cossu + 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 From c8fc0b3e0c78ee1679a07c886ad9530e95ec9044 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Wed, 2 Oct 2019 11:36:39 +0100 Subject: [PATCH 10/18] changed baryon interface --- Hadrons/Modules/MContraction/Baryon.hpp | 52 +++++++++++++++++++++++-- 1 file changed, 48 insertions(+), 4 deletions(-) diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 6366f7d2..ec01b48e 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -51,9 +51,11 @@ public: std::string, q3_src, std::string, GammaA, std::string, GammaB, +/* std::vector, quarks, + std::vector, prefactors,*/ std::string, quarks_snk, std::string, quarks_src, - int, parity, + int, parity, //=1 std::string, sink, std::string, output); }; @@ -123,18 +125,37 @@ template void TBaryon::setup(void) { envTmpLat(LatticeComplex, "c"); + envTmpLat(LatticeComplex, "c2"); } // execution /////////////////////////////////////////////////////////////////// template void TBaryon::execute(void) { + + //std::vector quarks = {"sud","dus"}; + //std::vector prefactors = {2.0, 1.0}; + std::vector quarks; + std::string qq = "sud dus"; + quarks = strToVec(qq); + std::vector prefactors; + std::string pp = "2.0 -1.0"; + prefactors = strToVec(pp); + + int nQ=quarks.size(); + for (int iQ1 = 0; iQ1 < nQ; iQ1++){ + for (int iQ2 = 0; iQ2 < nQ; iQ2++){ + LOG(Message) << prefactors[iQ1]*prefactors[iQ2] << "*<" << quarks[iQ1] << "|" << quarks[iQ2] << ">" << std::endl; + } + } + LOG(Message) << "Computing baryon contraction '" << getName() << "' < " << par().quarks_snk << " | " << par().quarks_src << " > using" << " quarks '" << par().q1_src << "', and a diquark formed of ('" << par().q2_src << "', and '" << par().q3_src << "') at the source and (Gamma^A,Gamma^B) = ( " << par().GammaA << " , " << par().GammaB << " ) and parity " << par().parity << "." << std::endl; envGetTmp(LatticeComplex, c); + envGetTmp(LatticeComplex, c2); Result result; int nt = env().getDim(Tp); result.corr.resize(nt); @@ -144,6 +165,7 @@ void TBaryon::execute(void) Gamma GammaB(ggB[0]); std::vector buf; TComplex cs; + TComplex ch; const int parity {par().parity}; const char * quarks_snk{par().quarks_snk.c_str()}; const char * quarks_src{par().quarks_src.c_str()}; @@ -159,7 +181,15 @@ void TBaryon::execute(void) LOG(Message) << "(propagator already sinked)" << std::endl; for (unsigned int t = 0; t < buf.size(); ++t) { - BaryonUtils::ContractBaryons_Sliced(q1_src[t],q2_src[t],q3_src[t],GammaA,GammaB,quarks_snk,quarks_src,parity,cs); + cs = Zero(); + for (int iQ1 = 0; iQ1 < nQ; iQ1++){ + for (int iQ2 = 0; iQ2 < nQ; iQ2++){ + BaryonUtils::ContractBaryons_Sliced(q1_src[t],q2_src[t],q3_src[t],GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,ch); + cs += prefactors[iQ1]*prefactors[iQ2]*ch; + } + } + + /* BaryonUtils::ContractBaryons_Sliced(q1_src[t],q2_src[t],q3_src[t],GammaA,GammaB,quarks_snk,quarks_src,parity,cs);*/ result.corr[t] = TensorRemove(cs); } } @@ -173,14 +203,28 @@ void TBaryon::execute(void) ns = vm().getModuleNamespace(env().getObjectModule(par().sink)); if (ns == "MSource") { - BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); + c=Zero(); + for (int iQ1 = 0; iQ1 < nQ; iQ1++){ + for (int iQ2 = 0; iQ2 < nQ; iQ2++){ + BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); + c+=prefactors[iQ1]*prefactors[iQ2]*c2; + } + } + // BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); PropagatorField1 &sink = envGet(PropagatorField1, par().sink); auto test = closure(trace(sink*c)); sliceSum(test, buf, Tp); } else if (ns == "MSink") { - BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); + c=Zero(); + for (int iQ1 = 0; iQ1 < nQ; iQ1++){ + for (int iQ2 = 0; iQ2 < nQ; iQ2++){ + BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); + c+=prefactors[iQ1]*prefactors[iQ2]*c2; + } + } + //BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); buf = sink(c); } From 89ef2b7dc28cf0f516dd56f7b6bc1c278da72296 Mon Sep 17 00:00:00 2001 From: Michael Marshall <43034299+mmphys@users.noreply.github.com> Date: Wed, 2 Oct 2019 13:20:07 +0100 Subject: [PATCH 11/18] Should compile everywhere now --- Grid/qcd/utils/BaryonUtils.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 9e8d53b6..5d462990 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -92,7 +92,7 @@ void BaryonUtils::baryon_site(const mobj &D1, static const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; - static const int epsilon_sgn[6]= {1,1,1,-1,-1,-1}; + static const Complex epsilon_sgn[6]= {1,1,1,-1,-1,-1}; auto gD1a = GammaA * GammaA * D1; auto gD1b = GammaA * g4 * GammaA * D1; From e280ec6b0bb000f45f837198fe257fb7a300ad06 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Wed, 2 Oct 2019 16:14:06 +0100 Subject: [PATCH 12/18] changed input parameters for easier use --- Hadrons/Modules/MContraction/Baryon.hpp | 79 +++++++++++-------------- 1 file changed, 33 insertions(+), 46 deletions(-) diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index ec01b48e..7d3f61ec 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -46,16 +46,14 @@ class BaryonPar: Serializable { public: GRID_SERIALIZABLE_CLASS_MEMBERS(BaryonPar, - std::string, q1_src, - std::string, q2_src, - std::string, q3_src, + std::string, q1, + std::string, q2, + std::string, q3, std::string, GammaA, std::string, GammaB, -/* std::vector, quarks, - std::vector, prefactors,*/ - std::string, quarks_snk, - std::string, quarks_src, - int, parity, //=1 + std::string, quarks, + std::string, prefactors, + std::string, parity, std::string, sink, std::string, output); }; @@ -107,7 +105,7 @@ TBaryon::TBaryon(const std::string name) template std::vector TBaryon::getInput(void) { - std::vector input = {par().q1_src, par().q2_src, par().q3_src, par().sink}; + std::vector input = {par().q1, par().q2, par().q3, par().sink}; return input; } @@ -133,26 +131,22 @@ template void TBaryon::execute(void) { - //std::vector quarks = {"sud","dus"}; - //std::vector prefactors = {2.0, 1.0}; - std::vector quarks; - std::string qq = "sud dus"; - quarks = strToVec(qq); - std::vector prefactors; - std::string pp = "2.0 -1.0"; - prefactors = strToVec(pp); - + std::vector quarks = strToVec(par().quarks); + std::vector prefactors = strToVec(par().prefactors); int nQ=quarks.size(); - for (int iQ1 = 0; iQ1 < nQ; iQ1++){ - for (int iQ2 = 0; iQ2 < nQ; iQ2++){ - LOG(Message) << prefactors[iQ1]*prefactors[iQ2] << "*<" << quarks[iQ1] << "|" << quarks[iQ2] << ">" << std::endl; - } - } + const int parity {par().parity.size()>0 ? std::stoi(par().parity) : 1}; - LOG(Message) << "Computing baryon contraction '" << getName() << "' < " << par().quarks_snk << " | " << par().quarks_src << " > using" - << " quarks '" << par().q1_src << "', and a diquark formed of ('" << par().q2_src << "', and '" - << par().q3_src << "') at the source and (Gamma^A,Gamma^B) = ( " << par().GammaA << " , " << par().GammaB - << " ) and parity " << par().parity << "." << std::endl; + assert(prefactors.size()==nQ && "number of prefactors needs to match number of quark-structures."); + for (int iQ = 0; iQ < nQ; iQ++) + assert(quarks[iQ].size()==3 && "quark-structures must consist of 3 quarks each."); + + LOG(Message) << "Computing baryon contractions '" << getName() << "'" << std::endl; + for (int iQ1 = 0; iQ1 < nQ; iQ1++) + for (int iQ2 = 0; iQ2 < nQ; iQ2++) + LOG(Message) << prefactors[iQ1]*prefactors[iQ2] << "*<" << quarks[iQ1] << "|" << quarks[iQ2] << ">" << std::endl; + LOG(Message) << " using quarks" << par().q1 << "', " << par().q2 << "', and '" + << par().q3 << "' and (Gamma^A,Gamma^B) = ( " << par().GammaA << " , " << par().GammaB + << " ) and parity " << parity << " using sink " << par().sink << "." << std::endl; envGetTmp(LatticeComplex, c); envGetTmp(LatticeComplex, c2); @@ -166,17 +160,14 @@ void TBaryon::execute(void) std::vector buf; TComplex cs; TComplex ch; - const int parity {par().parity}; - const char * quarks_snk{par().quarks_snk.c_str()}; - const char * quarks_src{par().quarks_src.c_str()}; - if (envHasType(SlicedPropagator1, par().q1_src) and - envHasType(SlicedPropagator2, par().q2_src) and - envHasType(SlicedPropagator3, par().q3_src)) + if (envHasType(SlicedPropagator1, par().q1) and + envHasType(SlicedPropagator2, par().q2) and + envHasType(SlicedPropagator3, par().q3)) { - auto &q1_src = envGet(SlicedPropagator1, par().q1_src); - auto &q2_src = envGet(SlicedPropagator2, par().q2_src); - auto &q3_src = envGet(SlicedPropagator3, par().q3_src); + auto &q1 = envGet(SlicedPropagator1, par().q1); + auto &q2 = envGet(SlicedPropagator2, par().q2); + auto &q3 = envGet(SlicedPropagator3, par().q3); LOG(Message) << "(propagator already sinked)" << std::endl; for (unsigned int t = 0; t < buf.size(); ++t) @@ -184,20 +175,18 @@ void TBaryon::execute(void) cs = Zero(); for (int iQ1 = 0; iQ1 < nQ; iQ1++){ for (int iQ2 = 0; iQ2 < nQ; iQ2++){ - BaryonUtils::ContractBaryons_Sliced(q1_src[t],q2_src[t],q3_src[t],GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,ch); + BaryonUtils::ContractBaryons_Sliced(q1[t],q2[t],q3[t],GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,ch); cs += prefactors[iQ1]*prefactors[iQ2]*ch; } } - - /* BaryonUtils::ContractBaryons_Sliced(q1_src[t],q2_src[t],q3_src[t],GammaA,GammaB,quarks_snk,quarks_src,parity,cs);*/ result.corr[t] = TensorRemove(cs); } } else { - auto &q1_src = envGet(PropagatorField1, par().q1_src); - auto &q2_src = envGet(PropagatorField2, par().q2_src); - auto &q3_src = envGet(PropagatorField3, par().q3_src); + auto &q1 = envGet(PropagatorField1, par().q1); + auto &q2 = envGet(PropagatorField2, par().q2); + auto &q3 = envGet(PropagatorField3, par().q3); std::string ns; ns = vm().getModuleNamespace(env().getObjectModule(par().sink)); @@ -206,11 +195,10 @@ void TBaryon::execute(void) c=Zero(); for (int iQ1 = 0; iQ1 < nQ; iQ1++){ for (int iQ2 = 0; iQ2 < nQ; iQ2++){ - BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); + BaryonUtils::ContractBaryons(q1,q2,q3,GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); c+=prefactors[iQ1]*prefactors[iQ2]*c2; } } - // BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); PropagatorField1 &sink = envGet(PropagatorField1, par().sink); auto test = closure(trace(sink*c)); sliceSum(test, buf, Tp); @@ -220,11 +208,10 @@ void TBaryon::execute(void) c=Zero(); for (int iQ1 = 0; iQ1 < nQ; iQ1++){ for (int iQ2 = 0; iQ2 < nQ; iQ2++){ - BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); + BaryonUtils::ContractBaryons(q1,q2,q3,GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); c+=prefactors[iQ1]*prefactors[iQ2]*c2; } } - //BaryonUtils::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c); SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); buf = sink(c); } From 2ce7f2b4d8578a26a259ef3a8634110a98300101 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 8 Oct 2019 13:19:47 +0100 Subject: [PATCH 13/18] suggested changes for 1st pull request implemented --- Grid/qcd/utils/BaryonUtils.h | 44 ++++++++++--------------- Hadrons/Modules/MContraction/Baryon.hpp | 18 ++++++++-- tests/hadrons/Test_hadrons_spectrum.cc | 6 ++++ 3 files changed, 39 insertions(+), 29 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 5d462990..b68055e0 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -31,8 +31,6 @@ NAMESPACE_BEGIN(Grid); -#undef DELTA_F_EQ_2 - template class BaryonUtils { @@ -43,6 +41,9 @@ public: typedef typename FImpl::SitePropagator pobj; typedef typename ComplexField::vector_object vobj; + + static constexpr int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; + static constexpr Complex epsilon_sgn[6]= {1,1,1,-1,-1,-1}; private: template @@ -52,7 +53,7 @@ public: const Gamma GammaA, const Gamma GammaB, const int parity, - const std::vector &wick_contractions, + const int * wick_contractions, robj &result); public: static void ContractBaryons(const PropagatorField &q1_src, @@ -76,6 +77,11 @@ public: robj &result); }; +template +constexpr int BaryonUtils::epsilon[6][3]; +template +constexpr Complex BaryonUtils::epsilon_sgn[6]; + template template void BaryonUtils::baryon_site(const mobj &D1, @@ -84,16 +90,12 @@ void BaryonUtils::baryon_site(const mobj &D1, const Gamma GammaA, const Gamma GammaB, const int parity, - const std::vector &wick_contraction, + const int * wick_contraction, robj &result) { Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) - - static const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; - static const Complex epsilon_sgn[6]= {1,1,1,-1,-1,-1}; - auto gD1a = GammaA * GammaA * D1; auto gD1b = GammaA * g4 * GammaA * D1; auto pD1 = 0.5* (gD1a + (double)parity * gD1b); @@ -176,7 +178,7 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, const int parity, ComplexField &baryon_corr) { - std::cout << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << std::endl; + std::cout << "Contraction <" << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << "|" << quarks_src[0] << quarks_src[1] << quarks_src[2] << ">" << std::endl; std::cout << "GammaA " << (GammaA.g) << std::endl; std::cout << "GammaB " << (GammaB.g) << std::endl; @@ -184,17 +186,10 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, GridBase *grid = q1_src.Grid(); - Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) - - static const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; - static const int epsilon_sgn[6]= {1,1,1,-1,-1,-1}; - std::vector wick_contraction = {0,0,0,0,0,0}; - + int wick_contraction[6]; for (int ie=0; ie < 6 ; ie++) - if (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) - wick_contraction[ie]=1; + wick_contraction[ie] = (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) ? 1 : 0; -// typedef typename ComplexField::vector_object vobj; auto vbaryon_corr= baryon_corr.View(); auto v1 = q1_src.View(); auto v2 = q2_src.View(); @@ -225,18 +220,15 @@ void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, const int parity, robj &result) { + std::cout << "Contraction <" << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << "|" << quarks_src[0] << quarks_src[1] << quarks_src[2] << ">" << std::endl; + std::cout << "GammaA " << (GammaA.g) << std::endl; + std::cout << "GammaB " << (GammaB.g) << std::endl; assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); - Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) - - static const int epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}}; - static const int epsilon_sgn[6]= {1,1,1,-1,-1,-1}; - std::vector wick_contraction = {0,0,0,0,0,0}; - + int wick_contraction[6]; for (int ie=0; ie < 6 ; ie++) - if (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) - wick_contraction[ie]=1; + wick_contraction[ie] = (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) ? 1 : 0; result=Zero(); baryon_site(D1,D2,D3,GammaA,GammaB,parity,wick_contraction,result); diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 7d3f61ec..097b5c71 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -71,6 +71,11 @@ public: { public: GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + Gamma::Algebra, gammaA, + Gamma::Algebra, gammaB, + std::string, quarks, + std::string, prefactors, + int, parity, std::vector, corr); }; public: @@ -144,15 +149,14 @@ void TBaryon::execute(void) for (int iQ1 = 0; iQ1 < nQ; iQ1++) for (int iQ2 = 0; iQ2 < nQ; iQ2++) LOG(Message) << prefactors[iQ1]*prefactors[iQ2] << "*<" << quarks[iQ1] << "|" << quarks[iQ2] << ">" << std::endl; - LOG(Message) << " using quarks" << par().q1 << "', " << par().q2 << "', and '" + LOG(Message) << " using quarks " << par().q1 << "', " << par().q2 << "', and '" << par().q3 << "' and (Gamma^A,Gamma^B) = ( " << par().GammaA << " , " << par().GammaB << " ) and parity " << parity << " using sink " << par().sink << "." << std::endl; + envGetTmp(LatticeComplex, c); envGetTmp(LatticeComplex, c2); - Result result; int nt = env().getDim(Tp); - result.corr.resize(nt); std::vector ggA = strToVec(par().GammaA); Gamma GammaA(ggA[0]); std::vector ggB = strToVec(par().GammaB); @@ -161,6 +165,14 @@ void TBaryon::execute(void) TComplex cs; TComplex ch; + Result result; + result.corr.resize(nt); + result.gammaA = ggA[0]; + result.gammaB = ggB[0]; + result.parity = parity; + result.quarks = par().quarks; + result.prefactors = par().prefactors; + if (envHasType(SlicedPropagator1, par().q1) and envHasType(SlicedPropagator2, par().q2) and envHasType(SlicedPropagator3, par().q3)) diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index af1dccd7..6c6194c9 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -46,6 +46,7 @@ int main(int argc, char *argv[]) // run setup /////////////////////////////////////////////////////////////// Application application; std::vector flavour = {"l", "s", "c1", "c2", "c3"}; + std::vector flavour_baryon = {"l", "s", "a", "b", "c"}; //needs to be a single character std::vector mass = {.01, .04, .2 , .25 , .3 }; // global parameters @@ -134,6 +135,11 @@ int main(int argc, char *argv[]) barPar.q1 = "Qpt_" + flavour[i]; barPar.q2 = "Qpt_" + flavour[j]; barPar.q3 = "Qpt_" + flavour[k]; + barPar.GammaA = "Identity"; + barPar.GammaB = "GammaZGamma5"; //C*GammaX + barPar.quarks = flavour_baryon[i] + flavour_baryon[j] + flavour_baryon[k]; + barPar.prefactors = "1.0"; + barPar.sink = "sink"; application.createModule( "baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar); } From 548b3bf43c6f381674d0c95901efbfc3a2cb9a1f Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Wed, 9 Oct 2019 14:52:33 +0100 Subject: [PATCH 14/18] second update to pull request --- Grid/qcd/utils/BaryonUtils.h | 178 +++++++++++++----------- Hadrons/Modules/MContraction/Baryon.hpp | 166 ++++++++++++++-------- tests/hadrons/Test_hadrons_spectrum.cc | 3 +- tests/qdpxx/Test_qdpxx_baryon.cc | 28 ++-- 4 files changed, 220 insertions(+), 155 deletions(-) diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index b68055e0..73d41422 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -50,29 +50,35 @@ public: static void baryon_site(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA, - const Gamma GammaB, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, const int parity, const int * wick_contractions, robj &result); public: - static void ContractBaryons(const PropagatorField &q1_src, - const PropagatorField &q2_src, - const PropagatorField &q3_src, - const Gamma GammaA, - const Gamma GammaB, - const char * quarks_snk, - const char * quarks_src, + static void ContractBaryons(const PropagatorField &q1_left, + const PropagatorField &q2_left, + const PropagatorField &q3_left, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const char * quarks_left, + const char * quarks_right, const int parity, ComplexField &baryon_corr); template static void ContractBaryons_Sliced(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA, - const Gamma GammaB, - const char * quarks_snk, - const char * quarks_src, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const char * quarks_left, + const char * quarks_right, const int parity, robj &result); }; @@ -87,8 +93,10 @@ template void BaryonUtils::baryon_site(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA, - const Gamma GammaB, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, const int parity, const int * wick_contraction, robj &result) @@ -96,71 +104,71 @@ void BaryonUtils::baryon_site(const mobj &D1, Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4) - auto gD1a = GammaA * GammaA * D1; - auto gD1b = GammaA * g4 * GammaA * D1; + auto gD1a = GammaA_left * GammaA_right * D1; + auto gD1b = GammaA_left * g4 * GammaA_right * D1; auto pD1 = 0.5* (gD1a + (double)parity * gD1b); - auto gD3 = GammaB * D3; + auto gD3 = GammaB_right * D3; - 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 ie_left=0; ie_left < 6 ; ie_left++){ + int a_left = epsilon[ie_left][0]; //a + int b_left = epsilon[ie_left][1]; //b + int c_left = epsilon[ie_left][2]; //c + for (int ie_right=0; ie_right < 6 ; ie_right++){ + int a_right = epsilon[ie_right][0]; //a' + int b_right = epsilon[ie_right][1]; //b' + int c_right = epsilon[ie_right][2]; //c' //This is the \delta_{456}^{123} part if (wick_contraction[0]){ - auto D2g = D2 * GammaB; - for (int alpha_snk=0; alpha_snk::baryon_site(const mobj &D1, } template -void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, - const PropagatorField &q2_src, - const PropagatorField &q3_src, - const Gamma GammaA, - const Gamma GammaB, - const char * quarks_snk, - const char * quarks_src, +void BaryonUtils::ContractBaryons(const PropagatorField &q1_left, + const PropagatorField &q2_left, + const PropagatorField &q3_left, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const char * quarks_left, + const char * quarks_right, const int parity, ComplexField &baryon_corr) { - std::cout << "Contraction <" << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << "|" << quarks_src[0] << quarks_src[1] << quarks_src[2] << ">" << std::endl; - std::cout << "GammaA " << (GammaA.g) << std::endl; - std::cout << "GammaB " << (GammaB.g) << std::endl; + std::cout << "Contraction <" << quarks_right[0] << quarks_right[1] << quarks_right[2] << "|" << quarks_left[0] << quarks_left[1] << quarks_left[2] << ">" << std::endl; + std::cout << "GammaA (left) " << (GammaA_left.g) << std::endl; + std::cout << "GammaB (left) " << (GammaB_left.g) << std::endl; + std::cout << "GammaA (right) " << (GammaA_right.g) << std::endl; + std::cout << "GammaB (right) " << (GammaB_right.g) << std::endl; assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); - GridBase *grid = q1_src.Grid(); + GridBase *grid = q1_left.Grid(); int wick_contraction[6]; for (int ie=0; ie < 6 ; ie++) - wick_contraction[ie] = (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) ? 1 : 0; + wick_contraction[ie] = (quarks_left[0] == quarks_right[epsilon[ie][0]] && quarks_left[1] == quarks_right[epsilon[ie][1]] && quarks_left[2] == quarks_right[epsilon[ie][2]]) ? 1 : 0; auto vbaryon_corr= baryon_corr.View(); - auto v1 = q1_src.View(); - auto v2 = q2_src.View(); - auto v3 = q3_src.View(); + auto v1 = q1_left.View(); + auto v2 = q2_left.View(); + auto v3 = q3_left.View(); // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { thread_for(ss,grid->oSites(),{ @@ -204,7 +216,7 @@ void BaryonUtils::ContractBaryons(const PropagatorField &q1_src, auto D3 = v3[ss]; vobj result=Zero(); - baryon_site(D1,D2,D3,GammaA,GammaB,parity,wick_contraction,result); + baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result); vbaryon_corr[ss] = result; } );//end loop over lattice sites } @@ -213,24 +225,28 @@ template void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, const mobj &D2, const mobj &D3, - const Gamma GammaA, - const Gamma GammaB, - const char * quarks_snk, - const char * quarks_src, + const Gamma GammaA_left, + const Gamma GammaB_left, + const Gamma GammaA_right, + const Gamma GammaB_right, + const char * quarks_left, + const char * quarks_right, const int parity, robj &result) { - std::cout << "Contraction <" << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] << "|" << quarks_src[0] << quarks_src[1] << quarks_src[2] << ">" << std::endl; - std::cout << "GammaA " << (GammaA.g) << std::endl; - std::cout << "GammaB " << (GammaB.g) << std::endl; + std::cout << "Contraction <" << quarks_right[0] << quarks_right[1] << quarks_right[2] << "|" << quarks_left[0] << quarks_left[1] << quarks_left[2] << ">" << std::endl; + std::cout << "GammaA (left) " << (GammaA_left.g) << std::endl; + std::cout << "GammaB (left) " << (GammaB_left.g) << std::endl; + std::cout << "GammaA (right) " << (GammaA_right.g) << std::endl; + std::cout << "GammaB (right) " << (GammaB_right.g) << std::endl; assert(parity==1 || parity == -1 && "Parity must be +1 or -1"); int wick_contraction[6]; for (int ie=0; ie < 6 ; ie++) - wick_contraction[ie] = (quarks_src[0] == quarks_snk[epsilon[ie][0]] && quarks_src[1] == quarks_snk[epsilon[ie][1]] && quarks_src[2] == quarks_snk[epsilon[ie][2]]) ? 1 : 0; + wick_contraction[ie] = (quarks_left[0] == quarks_right[epsilon[ie][0]] && quarks_left[1] == quarks_right[epsilon[ie][1]] && quarks_left[2] == quarks_right[epsilon[ie][2]]) ? 1 : 0; result=Zero(); - baryon_site(D1,D2,D3,GammaA,GammaB,parity,wick_contraction,result); + baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result); } NAMESPACE_END(Grid); diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 097b5c71..895db57b 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -42,6 +42,9 @@ BEGIN_HADRONS_NAMESPACE ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MContraction) +typedef std::pair GammaAB; +typedef std::pair GammaABPair; + class BaryonPar: Serializable { public: @@ -49,8 +52,7 @@ public: std::string, q1, std::string, q2, std::string, q3, - std::string, GammaA, - std::string, GammaB, + std::string, gammas, std::string, quarks, std::string, prefactors, std::string, parity, @@ -71,8 +73,10 @@ public: { public: GRID_SERIALIZABLE_CLASS_MEMBERS(Result, - Gamma::Algebra, gammaA, - Gamma::Algebra, gammaB, + Gamma::Algebra, gammaA_left, + Gamma::Algebra, gammaB_left, + Gamma::Algebra, gammaA_right, + Gamma::Algebra, gammaB_right, std::string, quarks, std::string, prefactors, int, parity, @@ -86,6 +90,7 @@ public: // dependency relation virtual std::vector getInput(void); virtual std::vector getOutput(void); + virtual void parseGammaString(std::vector &gammaList); protected: // setup virtual void setup(void); @@ -123,6 +128,33 @@ std::vector TBaryon::getOutput(void) return out; } +template +void TBaryon::parseGammaString(std::vector &gammaList) +{ + gammaList.clear(); + + std::string gammaString = par().gammas; + //Shorthands for standard baryon operators + gammaString = regex_replace(gammaString, std::regex("j12"),"Identity SigmaXZ"); + gammaString = regex_replace(gammaString, std::regex("j32X"),"Identity GammaZGamma5"); + gammaString = regex_replace(gammaString, std::regex("j32Y"),"Identity GammaT"); + gammaString = regex_replace(gammaString, std::regex("j32Z"),"Identity GammaXGamma5"); + //Shorthands for less common baryon operators + gammaString = regex_replace(gammaString, std::regex("j12_alt1"),"Gamma5 SigmaYT"); + gammaString = regex_replace(gammaString, std::regex("j12_alt2"),"Identity GammaYGamma5"); + + std::vector gamma_help; + gamma_help = strToVec(gammaString); + LOG(Message) << gamma_help.size() << " " << gamma_help.size()%2 << std::endl; + assert(gamma_help.size()%2==0 && "need even number of gamma-pairs."); + gammaList.resize(gamma_help.size()/2); + + for (int i = 0; i < gamma_help.size()/2; i++){ + gammaList[i].first=gamma_help[2*i]; + gammaList[i].second=gamma_help[2*i+1]; + } +} + // setup /////////////////////////////////////////////////////////////////////// template void TBaryon::setup(void) @@ -141,6 +173,9 @@ void TBaryon::execute(void) int nQ=quarks.size(); const int parity {par().parity.size()>0 ? std::stoi(par().parity) : 1}; + std::vector gammaList; + parseGammaString(gammaList); + assert(prefactors.size()==nQ && "number of prefactors needs to match number of quark-structures."); for (int iQ = 0; iQ < nQ; iQ++) assert(quarks[iQ].size()==3 && "quark-structures must consist of 3 quarks each."); @@ -149,29 +184,31 @@ void TBaryon::execute(void) for (int iQ1 = 0; iQ1 < nQ; iQ1++) for (int iQ2 = 0; iQ2 < nQ; iQ2++) LOG(Message) << prefactors[iQ1]*prefactors[iQ2] << "*<" << quarks[iQ1] << "|" << quarks[iQ2] << ">" << std::endl; - LOG(Message) << " using quarks " << par().q1 << "', " << par().q2 << "', and '" - << par().q3 << "' and (Gamma^A,Gamma^B) = ( " << par().GammaA << " , " << par().GammaB - << " ) and parity " << parity << " using sink " << par().sink << "." << std::endl; + LOG(Message) << " using quarks " << par().q1 << "', " << par().q2 << "', and '" << par().q3 << std::endl; + for (int iG = 0; iG < gammaList.size(); iG++) + LOG(Message) << "' with (Gamma^A,Gamma^B)_left = ( " << gammaList[iG].first.first << " , " << gammaList[iG].first.second << "') and (Gamma^A,Gamma^B)_right = ( " << gammaList[iG].second.first << " , " << gammaList[iG].second.second << ")" << std::endl; + LOG(Message) << "and parity " << parity << " using sink " << par().sink << "." << std::endl; - envGetTmp(LatticeComplex, c); envGetTmp(LatticeComplex, c2); int nt = env().getDim(Tp); - std::vector ggA = strToVec(par().GammaA); - Gamma GammaA(ggA[0]); - std::vector ggB = strToVec(par().GammaB); - Gamma GammaB(ggB[0]); std::vector buf; TComplex cs; TComplex ch; - Result result; - result.corr.resize(nt); - result.gammaA = ggA[0]; - result.gammaB = ggB[0]; - result.parity = parity; - result.quarks = par().quarks; - result.prefactors = par().prefactors; + std::vector result; + result.resize(gammaList.size()); + for (unsigned int i = 0; i < result.size(); ++i) + { + result[i].gammaA_left = gammaList[i].first.first; + result[i].gammaA_left = gammaList[i].first.first; + result[i].gammaB_right = gammaList[i].second.second; + result[i].gammaB_right = gammaList[i].second.second; + result[i].corr.resize(nt); + result[i].parity = parity; + result[i].quarks = par().quarks; + result[i].prefactors = par().prefactors; + } if (envHasType(SlicedPropagator1, par().q1) and envHasType(SlicedPropagator2, par().q2) and @@ -180,56 +217,71 @@ void TBaryon::execute(void) auto &q1 = envGet(SlicedPropagator1, par().q1); auto &q2 = envGet(SlicedPropagator2, par().q2); auto &q3 = envGet(SlicedPropagator3, par().q3); - - LOG(Message) << "(propagator already sinked)" << std::endl; - for (unsigned int t = 0; t < buf.size(); ++t) + for (unsigned int i = 0; i < result.size(); ++i) { - cs = Zero(); - for (int iQ1 = 0; iQ1 < nQ; iQ1++){ - for (int iQ2 = 0; iQ2 < nQ; iQ2++){ - BaryonUtils::ContractBaryons_Sliced(q1[t],q2[t],q3[t],GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,ch); - cs += prefactors[iQ1]*prefactors[iQ2]*ch; + Gamma gAl(gammaList[i].first.first); + Gamma gBl(gammaList[i].first.second); + Gamma gAr(gammaList[i].second.first); + Gamma gBr(gammaList[i].second.second); + + LOG(Message) << "(propagator already sinked)" << std::endl; + for (unsigned int t = 0; t < buf.size(); ++t) + { + cs = Zero(); + for (int iQ1 = 0; iQ1 < nQ; iQ1++){ + for (int iQ2 = 0; iQ2 < nQ; iQ2++){ + BaryonUtils::ContractBaryons_Sliced(q1[t],q2[t],q3[t],gAl,gBl,gAr,gBr,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,ch); + cs += prefactors[iQ1]*prefactors[iQ2]*ch; + } } + result[i].corr[t] = TensorRemove(cs); } - result.corr[t] = TensorRemove(cs); } } else { - auto &q1 = envGet(PropagatorField1, par().q1); - auto &q2 = envGet(PropagatorField2, par().q2); - auto &q3 = envGet(PropagatorField3, par().q3); - std::string ns; + auto &q1 = envGet(PropagatorField1, par().q1); + auto &q2 = envGet(PropagatorField2, par().q2); + auto &q3 = envGet(PropagatorField3, par().q3); + for (unsigned int i = 0; i < result.size(); ++i) + { + Gamma gAl(gammaList[i].first.first); + Gamma gBl(gammaList[i].first.second); + Gamma gAr(gammaList[i].second.first); + Gamma gBr(gammaList[i].second.second); + + std::string ns; - ns = vm().getModuleNamespace(env().getObjectModule(par().sink)); - if (ns == "MSource") - { - c=Zero(); - for (int iQ1 = 0; iQ1 < nQ; iQ1++){ - for (int iQ2 = 0; iQ2 < nQ; iQ2++){ - BaryonUtils::ContractBaryons(q1,q2,q3,GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); - c+=prefactors[iQ1]*prefactors[iQ2]*c2; + ns = vm().getModuleNamespace(env().getObjectModule(par().sink)); + if (ns == "MSource") + { + c=Zero(); + for (int iQ1 = 0; iQ1 < nQ; iQ1++){ + for (int iQ2 = 0; iQ2 < nQ; iQ2++){ + BaryonUtils::ContractBaryons(q1,q2,q3,gAl,gBl,gAr,gBr,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); + c+=prefactors[iQ1]*prefactors[iQ2]*c2; + } } + PropagatorField1 &sink = envGet(PropagatorField1, par().sink); + auto test = closure(trace(sink*c)); + sliceSum(test, buf, Tp); } - PropagatorField1 &sink = envGet(PropagatorField1, par().sink); - auto test = closure(trace(sink*c)); - sliceSum(test, buf, Tp); - } - else if (ns == "MSink") - { - c=Zero(); - for (int iQ1 = 0; iQ1 < nQ; iQ1++){ - for (int iQ2 = 0; iQ2 < nQ; iQ2++){ - BaryonUtils::ContractBaryons(q1,q2,q3,GammaA,GammaB,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); - c+=prefactors[iQ1]*prefactors[iQ2]*c2; + else if (ns == "MSink") + { + c=Zero(); + for (int iQ1 = 0; iQ1 < nQ; iQ1++){ + for (int iQ2 = 0; iQ2 < nQ; iQ2++){ + BaryonUtils::ContractBaryons(q1,q2,q3,gAl,gBl,gAr,gBr,quarks[iQ1].c_str(),quarks[iQ2].c_str(),parity,c2); + c+=prefactors[iQ1]*prefactors[iQ2]*c2; + } } + SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); + buf = sink(c); + } + for (unsigned int t = 0; t < buf.size(); ++t) + { + result[i].corr[t] = TensorRemove(buf[t]); } - SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); - buf = sink(c); - } - for (unsigned int t = 0; t < buf.size(); ++t) - { - result.corr[t] = TensorRemove(buf[t]); } } diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index 6c6194c9..810ef940 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -135,8 +135,7 @@ int main(int argc, char *argv[]) barPar.q1 = "Qpt_" + flavour[i]; barPar.q2 = "Qpt_" + flavour[j]; barPar.q3 = "Qpt_" + flavour[k]; - barPar.GammaA = "Identity"; - barPar.GammaB = "GammaZGamma5"; //C*GammaX + barPar.gammas = "(j12) (j12) (j32X) (j32Y)"; barPar.quarks = flavour_baryon[i] + flavour_baryon[j] + flavour_baryon[k]; barPar.prefactors = "1.0"; barPar.sink = "sink"; diff --git a/tests/qdpxx/Test_qdpxx_baryon.cc b/tests/qdpxx/Test_qdpxx_baryon.cc index b19864fb..ad51acfc 100644 --- a/tests/qdpxx/Test_qdpxx_baryon.cc +++ b/tests/qdpxx/Test_qdpxx_baryon.cc @@ -479,51 +479,49 @@ void calc_grid(Grid::LatticeGaugeField &Umu, Grid::LatticePropagator &qU, Grid:: Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu.Grid(); - Grid::Gamma G_A = Grid::Gamma(Grid::Gamma::Algebra::Identity); Grid::Gamma G_B = Grid::Gamma(Grid::Gamma::Algebra::GammaZGamma5); // OmegaX: C*GammaX = i* GammaZ*Gamma5 - char quarks[] = "sss"; Grid::LatticeComplex c(UGrid); Grid::LatticeComplex c1(UGrid); if(! baryon.compare("OmegaX")){ - BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,quarks,quarks,1,c); + BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,G_A,G_B,"sss","sss",1,c); c*=0.5; std::cout << "Grid-Omega factor 2 larger than Chroma-Omega!!!" << std::endl; } else if (! baryon.compare("OmegaY")){ G_B = Grid::Gamma(Grid::Gamma::Algebra::GammaT); - BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,quarks,quarks,1,c); + BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,G_A,G_B,"sss","sss",1,c); c*=0.5; std::cout << "Grid-Omega factor 2 larger than Chroma-Omega!!!" << std::endl; } else if (! baryon.compare("OmegaZ")){ G_B = Grid::Gamma(Grid::Gamma::Algebra::GammaXGamma5); - BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,quarks,quarks,1,c); + BaryonUtils::ContractBaryons(qS,qS,qS,G_A,G_B,G_A,G_B,"sss","sss",1,c); c*=0.5; std::cout << "Grid-Omega factor 2 larger than Chroma-Omega!!!" << std::endl; } else if (! baryon.compare("Proton")){ G_B = Grid::Gamma(Grid::Gamma::Algebra::SigmaXZ); - BaryonUtils::ContractBaryons(qU,qD,qU,G_A,G_B,"udu","udu",1,c); + BaryonUtils::ContractBaryons(qU,qD,qU,G_A,G_B,G_A,G_B,"udu","udu",1,c); std::cout << "UKHadron-Proton has flipped diquarks in original code." << std::endl; } else if (! baryon.compare("Lambda")){ G_B = Grid::Gamma(Grid::Gamma::Algebra::SigmaXZ); - BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,"sud","sud",1,c1); //s + BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,G_A,G_B,"sud","sud",1,c1); //s c = 4.*c1; - BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,"dus","dus",1,c1); //d + BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,G_A,G_B,"dus","dus",1,c1); //d c += c1; - BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,"uds","uds",1,c1); //u + BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,G_A,G_B,"uds","uds",1,c1); //u c += c1; - BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,"sud","dus",1,c1); //(sud) + BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,G_A,G_B,"dus","sud",1,c1); //(sud) c += 2.*c1; - BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,"sud","uds",1,c1); //(sdu) + BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,G_A,G_B,"uds","sud",1,c1); //(sdu) c -= 2.*c1; - BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,"dus","sud",1,c1); //(dus) + BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,G_A,G_B,"sud","dus",1,c1); //(dus) c += 2.*c1; - BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,"dus","uds",1,c1); //-(dsu) + BaryonUtils::ContractBaryons(qU,qD,qS,G_A,G_B,G_A,G_B,"uds","dus",1,c1); //-(dsu) c -= c1; - BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,"uds","sud",1,c1); //(uds) + BaryonUtils::ContractBaryons(qS,qU,qD,G_A,G_B,G_A,G_B,"sud","uds",1,c1); //(uds) c -= 2.*c1; - BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,"uds","dus",1,c1); //-(usd) + BaryonUtils::ContractBaryons(qD,qU,qS,G_A,G_B,G_A,G_B,"dus","uds",1,c1); //-(usd) c -= c1; std::cout << "UKHadron-Lambda has flipped diquarks in original code." << std::endl; } else { From 2dee4791dba5302a2b7dc0583b285d2d4823528a Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Wed, 9 Oct 2019 17:56:09 +0100 Subject: [PATCH 15/18] baryon input strings are now pairs of pairs of gammas - still ugly!! --- Hadrons/Modules/MContraction/Baryon.hpp | 60 ++++++++++++++++++------- 1 file changed, 43 insertions(+), 17 deletions(-) diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 895db57b..87c16285 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -135,23 +135,49 @@ void TBaryon::parseGammaString(std::vector & std::string gammaString = par().gammas; //Shorthands for standard baryon operators - gammaString = regex_replace(gammaString, std::regex("j12"),"Identity SigmaXZ"); - gammaString = regex_replace(gammaString, std::regex("j32X"),"Identity GammaZGamma5"); - gammaString = regex_replace(gammaString, std::regex("j32Y"),"Identity GammaT"); - gammaString = regex_replace(gammaString, std::regex("j32Z"),"Identity GammaXGamma5"); + gammaString = regex_replace(gammaString, std::regex("j12"),"(Identity SigmaXZ)"); + gammaString = regex_replace(gammaString, std::regex("j32X"),"(Identity GammaZGamma5)"); + gammaString = regex_replace(gammaString, std::regex("j32Y"),"(Identity GammaT)"); + gammaString = regex_replace(gammaString, std::regex("j32Z"),"(Identity GammaXGamma5)"); //Shorthands for less common baryon operators - gammaString = regex_replace(gammaString, std::regex("j12_alt1"),"Gamma5 SigmaYT"); - gammaString = regex_replace(gammaString, std::regex("j12_alt2"),"Identity GammaYGamma5"); - - std::vector gamma_help; - gamma_help = strToVec(gammaString); - LOG(Message) << gamma_help.size() << " " << gamma_help.size()%2 << std::endl; - assert(gamma_help.size()%2==0 && "need even number of gamma-pairs."); - gammaList.resize(gamma_help.size()/2); + gammaString = regex_replace(gammaString, std::regex("j12_alt1"),"(Gamma5 SigmaYT)"); + gammaString = regex_replace(gammaString, std::regex("j12_alt2"),"(Identity GammaYGamma5)"); + + //A single gamma matrix + std::regex rex_g("([0-9a-zA-Z]+)"); + //The full string we expect + std::regex rex("( *\\(( *\\(([0-9a-zA-Z]+) +([0-9a-zA-Z]+) *\\)){2} *\\) *)+"); + std::smatch sm; + std::regex_match(gammaString, sm, rex); + assert(sm[0].matched && "invalid gamma structure."); - for (int i = 0; i < gamma_help.size()/2; i++){ - gammaList[i].first=gamma_help[2*i]; - gammaList[i].second=gamma_help[2*i+1]; + auto gamma_begin = std::sregex_iterator(gammaString.begin(), gammaString.end(), rex_g); + auto gamma_end = std::sregex_iterator(); + + //couldn't find out how to count the size in the iterator, other than looping through it... + int nGamma=0; + for (std::sregex_iterator i = gamma_begin; i != gamma_end; ++i) { + nGamma++; + } + gammaList.resize(nGamma/4); + std::vector gS; + gS.resize(nGamma); + //even more ugly workarounds here... + int iG=0; + for (std::sregex_iterator i = gamma_begin; i != gamma_end; ++i) { + std::smatch match = *i; + gS[iG] = match.str(); + iG++; + } + for (int i = 0; i < gammaList.size(); i++){ + std::vector gS1 = strToVec(gS[4*i]); + std::vector gS2 = strToVec(gS[4*i+1]); + std::vector gS3 = strToVec(gS[4*i+2]); + std::vector gS4 = strToVec(gS[4*i+3]); + gammaList[i].first.first=gS1[0]; + gammaList[i].first.second=gS2[0]; + gammaList[i].second.first=gS3[0]; + gammaList[i].second.second=gS4[0]; } } @@ -201,8 +227,8 @@ void TBaryon::execute(void) for (unsigned int i = 0; i < result.size(); ++i) { result[i].gammaA_left = gammaList[i].first.first; - result[i].gammaA_left = gammaList[i].first.first; - result[i].gammaB_right = gammaList[i].second.second; + result[i].gammaB_left = gammaList[i].first.second; + result[i].gammaA_right = gammaList[i].second.first; result[i].gammaB_right = gammaList[i].second.second; result[i].corr.resize(nt); result[i].parity = parity; From aa62ca90467d9806878dcc4bbf1b57cd6228491b Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 10 Oct 2019 11:07:20 +0100 Subject: [PATCH 16/18] chroma-regression test now prints difference correctly --- tests/qdpxx/Test_qdpxx_baryon.cc | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/tests/qdpxx/Test_qdpxx_baryon.cc b/tests/qdpxx/Test_qdpxx_baryon.cc index ad51acfc..a1d8f738 100644 --- a/tests/qdpxx/Test_qdpxx_baryon.cc +++ b/tests/qdpxx/Test_qdpxx_baryon.cc @@ -572,7 +572,7 @@ int main(int argc, char **argv) PropagatorField strange(UGrid); std::vector res_chroma; std::vector res_grid; - std::vector difference; + Grid::Complex res_chroma_g; std::vector baryons({"OmegaX","OmegaY","OmegaZ","Proton","Lambda"}); int nBaryon=baryons.size(); @@ -592,9 +592,10 @@ int main(int argc, char **argv) for(int t=0;t Date: Mon, 14 Oct 2019 13:41:08 +0100 Subject: [PATCH 17/18] updated syntac in Test_hadrons_spectrum --- tests/hadrons/Test_hadrons_spectrum.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/hadrons/Test_hadrons_spectrum.cc b/tests/hadrons/Test_hadrons_spectrum.cc index 810ef940..7650397b 100644 --- a/tests/hadrons/Test_hadrons_spectrum.cc +++ b/tests/hadrons/Test_hadrons_spectrum.cc @@ -135,7 +135,7 @@ int main(int argc, char *argv[]) barPar.q1 = "Qpt_" + flavour[i]; barPar.q2 = "Qpt_" + flavour[j]; barPar.q3 = "Qpt_" + flavour[k]; - barPar.gammas = "(j12) (j12) (j32X) (j32Y)"; + barPar.gammas = "(j12 j12) (j32X j32Y)"; barPar.quarks = flavour_baryon[i] + flavour_baryon[j] + flavour_baryon[k]; barPar.prefactors = "1.0"; barPar.sink = "sink"; From 3c702b510bce5430653353da363eac0eff9fe1d1 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Tue, 15 Oct 2019 18:48:51 +0100 Subject: [PATCH 18/18] result layout standardised, iterator size more elegant --- Hadrons/Modules/MContraction/Baryon.hpp | 58 ++++++++++++++----------- 1 file changed, 33 insertions(+), 25 deletions(-) diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp index 87c16285..9b48ae0a 100644 --- a/Hadrons/Modules/MContraction/Baryon.hpp +++ b/Hadrons/Modules/MContraction/Baryon.hpp @@ -69,19 +69,19 @@ public: FERM_TYPE_ALIASES(FImpl3, 3); BASIC_TYPE_ALIASES(ScalarImplCR, Scalar); SINK_TYPE_ALIASES(Scalar); - class Result: Serializable + class Metadata: Serializable { public: - GRID_SERIALIZABLE_CLASS_MEMBERS(Result, + GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata, Gamma::Algebra, gammaA_left, Gamma::Algebra, gammaB_left, Gamma::Algebra, gammaA_right, Gamma::Algebra, gammaB_right, std::string, quarks, std::string, prefactors, - int, parity, - std::vector, corr); + int, parity); }; + typedef Correlator Result; public: // constructor TBaryon(const std::string name); @@ -136,11 +136,11 @@ void TBaryon::parseGammaString(std::vector & std::string gammaString = par().gammas; //Shorthands for standard baryon operators gammaString = regex_replace(gammaString, std::regex("j12"),"(Identity SigmaXZ)"); - gammaString = regex_replace(gammaString, std::regex("j32X"),"(Identity GammaZGamma5)"); + gammaString = regex_replace(gammaString, std::regex("j32X"),"(Identity MinusGammaZGamma5)"); gammaString = regex_replace(gammaString, std::regex("j32Y"),"(Identity GammaT)"); gammaString = regex_replace(gammaString, std::regex("j32Z"),"(Identity GammaXGamma5)"); //Shorthands for less common baryon operators - gammaString = regex_replace(gammaString, std::regex("j12_alt1"),"(Gamma5 SigmaYT)"); + gammaString = regex_replace(gammaString, std::regex("j12_alt1"),"(Gamma5 MinusSigmaYT)"); gammaString = regex_replace(gammaString, std::regex("j12_alt2"),"(Identity GammaYGamma5)"); //A single gamma matrix @@ -154,11 +154,13 @@ void TBaryon::parseGammaString(std::vector & auto gamma_begin = std::sregex_iterator(gammaString.begin(), gammaString.end(), rex_g); auto gamma_end = std::sregex_iterator(); + int nGamma = std::distance(gamma_begin, gamma_end); //couldn't find out how to count the size in the iterator, other than looping through it... - int nGamma=0; + /* int nGamma=0; for (std::sregex_iterator i = gamma_begin; i != gamma_end; ++i) { nGamma++; - } + } +*/ gammaList.resize(nGamma/4); std::vector gS; gS.resize(nGamma); @@ -222,19 +224,11 @@ void TBaryon::execute(void) TComplex cs; TComplex ch; - std::vector result; - result.resize(gammaList.size()); - for (unsigned int i = 0; i < result.size(); ++i) - { - result[i].gammaA_left = gammaList[i].first.first; - result[i].gammaB_left = gammaList[i].first.second; - result[i].gammaA_right = gammaList[i].second.first; - result[i].gammaB_right = gammaList[i].second.second; - result[i].corr.resize(nt); - result[i].parity = parity; - result[i].quarks = par().quarks; - result[i].prefactors = par().prefactors; - } + std::vector result; + Result r; + r.info.parity = parity; + r.info.quarks = par().quarks; + r.info.prefactors = par().prefactors; if (envHasType(SlicedPropagator1, par().q1) and envHasType(SlicedPropagator2, par().q2) and @@ -243,14 +237,20 @@ void TBaryon::execute(void) auto &q1 = envGet(SlicedPropagator1, par().q1); auto &q2 = envGet(SlicedPropagator2, par().q2); auto &q3 = envGet(SlicedPropagator3, par().q3); - for (unsigned int i = 0; i < result.size(); ++i) + for (unsigned int i = 0; i < gammaList.size(); ++i) { + r.info.gammaA_left = gammaList[i].first.first; + r.info.gammaB_left = gammaList[i].first.second; + r.info.gammaA_right = gammaList[i].second.first; + r.info.gammaB_right = gammaList[i].second.second; + Gamma gAl(gammaList[i].first.first); Gamma gBl(gammaList[i].first.second); Gamma gAr(gammaList[i].second.first); Gamma gBr(gammaList[i].second.second); LOG(Message) << "(propagator already sinked)" << std::endl; + r.corr.clear(); for (unsigned int t = 0; t < buf.size(); ++t) { cs = Zero(); @@ -260,8 +260,9 @@ void TBaryon::execute(void) cs += prefactors[iQ1]*prefactors[iQ2]*ch; } } - result[i].corr[t] = TensorRemove(cs); + r.corr.push_back(TensorRemove(cs)); } + result.push_back(r); } } else @@ -269,8 +270,13 @@ void TBaryon::execute(void) auto &q1 = envGet(PropagatorField1, par().q1); auto &q2 = envGet(PropagatorField2, par().q2); auto &q3 = envGet(PropagatorField3, par().q3); - for (unsigned int i = 0; i < result.size(); ++i) + for (unsigned int i = 0; i < gammaList.size(); ++i) { + r.info.gammaA_left = gammaList[i].first.first; + r.info.gammaB_left = gammaList[i].first.second; + r.info.gammaA_right = gammaList[i].second.first; + r.info.gammaB_right = gammaList[i].second.second; + Gamma gAl(gammaList[i].first.first); Gamma gBl(gammaList[i].first.second); Gamma gAr(gammaList[i].second.first); @@ -304,10 +310,12 @@ void TBaryon::execute(void) SinkFnScalar &sink = envGet(SinkFnScalar, par().sink); buf = sink(c); } + r.corr.clear(); for (unsigned int t = 0; t < buf.size(); ++t) { - result[i].corr[t] = TensorRemove(buf[t]); + r.corr.push_back(TensorRemove(buf[t])); } + result.push_back(r); } }