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