portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-04 19:25:56 +01:00
Code Releases Activity

added Baryon code

Browse Source
This commit is contained in:
Felix Erben 2019-09-16 17:20:54 +01:00
parent b473405652
commit dab8c01c3d
3 changed files with 826 additions and 16 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

176
Grid/qcd/utils/BaryonUtils.h Normal file
View File

@ -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 <felix.erben@ed.ac.uk>
Author: Michael Marshall <Michael.Marshall@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#pragma once
//#include <Grid/Hadrons/Global.hpp>
#include <Grid/Eigen/unsupported/CXX11/Tensor>
NAMESPACE_BEGIN(Grid);
#undef DELTA_F_EQ_2
template <typename FImpl>
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<class FImpl>
void BaryonUtils<FImpl>::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<std::vector<int>> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
std::vector<int> epsilon_sgn = {1,1,1,-1,-1,-1};
std::vector<int> 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<vobj> vbaryon_corr{ baryon_corr };
accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
LatticeView<pobj> v1(q1_src);
LatticeView<pobj> v2(q2_src);
LatticeView<pobj> 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<Ns; alpha_snk++){
for (int beta_src=0; beta_src<Ns; beta_src++){
for (int gamma_src=0; gamma_src<Ns; gamma_src++){
result()()() += epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1()(gamma_src,gamma_src)(c_snk,c_src)*D2g()(alpha_snk,beta_src)(a_snk,a_src)*gD3()(alpha_snk,beta_src)(b_snk,b_src);
}}}
}
//This is the \delta_{456}^{231} part
if (wick_contraction[1]){
auto pD1g = pD1 * GammaB;
for (int alpha_snk=0; alpha_snk<Ns; alpha_snk++){
for (int beta_src=0; beta_src<Ns; beta_src++){
for (int gamma_src=0; gamma_src<Ns; gamma_src++){
result()()() += epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1g()(gamma_src,beta_src)(c_snk,a_src)*D2()(alpha_snk,beta_src)(a_snk,b_src)*gD3()(alpha_snk,gamma_src)(b_snk,c_src);
}}}
}
//This is the \delta_{456}^{312} part
if (wick_contraction[2]){
auto gD3g = gD3 * GammaB;
for (int alpha_snk=0; alpha_snk<Ns; alpha_snk++){
for (int beta_src=0; beta_src<Ns; beta_src++){
for (int gamma_src=0; gamma_src<Ns; gamma_src++){
result()()() += epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1()(gamma_src,beta_src)(c_snk,b_src)*D2()(alpha_snk,gamma_src)(a_snk,c_src)*gD3g()(alpha_snk,beta_src)(b_snk,a_src);
}}}
}
//This is the \delta_{456}^{132} part
if (wick_contraction[3]){
auto gD3g = gD3 * GammaB;
for (int alpha_snk=0; alpha_snk<Ns; alpha_snk++){
for (int beta_src=0; beta_src<Ns; beta_src++){
for (int gamma_src=0; gamma_src<Ns; gamma_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;
for (int alpha_snk=0; alpha_snk<Ns; alpha_snk++){
for (int beta_src=0; beta_src<Ns; beta_src++){
for (int gamma_src=0; gamma_src<Ns; gamma_src++){
result()()() -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1()(gamma_src,beta_src)(c_snk,b_src)*D2g()(alpha_snk,beta_src)(a_snk,a_src)*gD3()(alpha_snk,gamma_src)(b_snk,c_src);
}}}
}
//This is the \delta_{456}^{213} part
if (wick_contraction[5]){
auto pD1g = pD1 * GammaB;
for (int alpha_snk=0; alpha_snk<Ns; alpha_snk++){
for (int beta_src=0; beta_src<Ns; beta_src++){
for (int gamma_src=0; gamma_src<Ns; gamma_src++){
result()()() -= epsilon_sgn[ie_src] * epsilon_sgn[ie_snk] * pD1g()(gamma_src,beta_src)(c_snk,a_src)*D2()(alpha_snk,gamma_src)(a_snk,c_src)*gD3()(alpha_snk,beta_src)(b_snk,b_src);
}}}
}
}
}
vbaryon_corr[ss] = result;
} );//end loop over lattice sites
}
NAMESPACE_END(Grid);

61
Hadrons/Modules/MContraction/Baryon.hpp
View File

@ -7,7 +7,7 @@ Source file: Hadrons/Modules/MContraction/Baryon.hpp
Copyright (C) 2015-2019 Copyright (C) 2015-2019
Author: Antonin Portelli <antonin.portelli@me.com> Author: Antonin Portelli <antonin.portelli@me.com>
Author: Lanny91 <andrew.lawson@gmail.com> Author: Felix Erben <felix.erben@ed.ac.uk>
This program is free software; you can redistribute it and/or modify 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 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 <Hadrons/Global.hpp> #include <Hadrons/Global.hpp>
#include <Hadrons/Module.hpp> #include <Hadrons/Module.hpp>
#include <Hadrons/ModuleFactory.hpp> #include <Hadrons/ModuleFactory.hpp>
#include <Grid/qcd/utils/BaryonUtils.h>
BEGIN_HADRONS_NAMESPACE BEGIN_HADRONS_NAMESPACE
@ -45,9 +46,16 @@ class BaryonPar: Serializable
{ {
public: public:
GRID_SERIALIZABLE_CLASS_MEMBERS(BaryonPar, GRID_SERIALIZABLE_CLASS_MEMBERS(BaryonPar,
std::string, q1, std::string, q1_src,
std::string, q2, std::string, q2_src,
std::string, q3, 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); std::string, output);
}; };
@ -62,7 +70,7 @@ public:
{ {
public: public:
GRID_SERIALIZABLE_CLASS_MEMBERS(Result, GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
std::vector<std::vector<std::vector<Complex>>>, corr); std::vector<Complex>, corr);
}; };
public: public:
// constructor // constructor
@ -77,6 +85,8 @@ protected:
virtual void setup(void); virtual void setup(void);
// execution // execution
virtual void execute(void); virtual void execute(void);
// Which gamma algebra was specified
Gamma::Algebra al;
}; };
MODULE_REGISTER_TMP(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction); MODULE_REGISTER_TMP(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
@ -94,7 +104,7 @@ TBaryon<FImpl1, FImpl2, FImpl3>::TBaryon(const std::string name)
template <typename FImpl1, typename FImpl2, typename FImpl3> template <typename FImpl1, typename FImpl2, typename FImpl3>
std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getInput(void) std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getInput(void)
{ {
std::vector<std::string> input = {par().q1, par().q2, par().q3}; std::vector<std::string> input = {par().q1_src, par().q2_src, par().q3_src};
return input; return input;
} }
@ -118,19 +128,38 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::setup(void)
template <typename FImpl1, typename FImpl2, typename FImpl3> template <typename FImpl1, typename FImpl2, typename FImpl3>
void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void) void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
{ {
LOG(Message) << "Computing baryon contractions '" << getName() << "' using" LOG(Message) << "Computing baryon contraction '" << getName() << "' < " << par().quarks_snk << " | " << par().quarks_src << " > using"
<< " quarks '" << par().q1 << "', '" << par().q2 << "', and '" << " quarks '" << par().q1_src << "', and a diquark formed of ('" << par().q2_src << "', and '"
<< par().q3 << "'" << std::endl; << 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 &q1_src = envGet(PropagatorField1, par().q1_src);
auto &q2 = envGet(PropagatorField2, par().q2); auto &q2_src = envGet(PropagatorField2, par().q2_src);
auto &q3 = envGet(PropagatorField3, par().q2); auto &q3_src = envGet(PropagatorField3, par().q3_src);
envGetTmp(LatticeComplex, c); envGetTmp(LatticeComplex, c);
Result result; Result result;
int nt = env().getDim(Tp);
// FIXME: do contractions result.corr.resize(nt);
std::vector<Gamma::Algebra> ggA = strToVec<Gamma::Algebra>(par().GammaA);
// saveResult(par().output, "meson", result); Gamma GammaA = ggA[0];
std::vector<Gamma::Algebra> ggB = strToVec<Gamma::Algebra>(par().GammaB);
Gamma GammaB = ggB[0];
std::vector<TComplex> 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<FIMPL>::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 END_MODULE_NAMESPACE

605
tests/qdpxx/Test_qdpxx_baryon.cc Normal file
View File

@ -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 <guido.cossu@ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <chroma.h>
#include <Grid/Grid.h>
#include <Grid/qcd/utils/BaryonUtils.h>
typedef Grid::LatticeGaugeField GaugeField;
namespace Chroma
{
class ChromaWrapper
{
public:
typedef multi1d<LatticeColorMatrix> U;
typedef LatticeFermion T4;
static void ImportGauge(GaugeField &gr,
QDP::multi1d<QDP::LatticeColorMatrix> &ch)
{
Grid::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> 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<QDP::LatticeColorMatrix> &ch)
{
Grid::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> 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<Grid::WilsonImplR>(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<Grid::WilsonImplR>(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<int> x(4); // explicit 4d fermions in Grid
QDP::multi1d<int> 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<int> x(4); // 4d fermions
QDP::multi1d<int> 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<int> 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<QDP::Complex> &res, std::string baryon)
{
QDP::multi1d<QDP::LatticeColorMatrix> 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<Grid::TComplex> 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<DComplex> 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[ <us>d + <ds>u + 4<ud>s - (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 <us>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 <ds>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 4<ud>s -- 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 <us>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 <ds>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 4<ud>s -- 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<DComplex> 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<Grid::Complex> &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<Grid::WilsonImplR>::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<Grid::WilsonImplR>::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<Grid::WilsonImplR>::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<Grid::WilsonImplR>::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<Grid::WilsonImplR>::ContractBaryons(qS,qU,qD,G_A,G_B,"sud","sud",1,c1); //<ud>s
c = 4.*c1;
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qD,qU,qS,G_A,G_B,"dus","dus",1,c1); //<us>d
c += c1;
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qU,qD,qS,G_A,G_B,"uds","uds",1,c1); //<ds>u
c += c1;
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qD,qU,qS,G_A,G_B,"sud","dus",1,c1); //(sud)
c += 2.*c1;
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qU,qD,qS,G_A,G_B,"sud","uds",1,c1); //(sdu)
c -= 2.*c1;
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qS,qU,qD,G_A,G_B,"dus","sud",1,c1); //(dus)
c += 2.*c1;
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qU,qD,qS,G_A,G_B,"dus","uds",1,c1); //-(dsu)
c -= c1;
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qS,qU,qD,G_A,G_B,"uds","sud",1,c1); //(uds)
c -= 2.*c1;
BaryonUtils<Grid::WilsonImplR>::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<Grid::TComplex> 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<int> 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<ComplexD> res_chroma;
std::vector<Grid::Complex> res_grid;
std::vector<Grid::Complex> difference;
std::vector<std::string> 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<res_chroma.size();t++){
std::cout << " Chroma baryon "<<t<<" "<< res_chroma[t] << std::endl;
}
calc_grid(Ug, up, down, strange, res_grid,baryons[iB]);
for(int t=0;t<res_chroma.size();t++){
std::cout << " Grid baryon "<<t<<" "<< res_grid[t] << std::endl;
}
/* for(int t=0;t<res_chroma.size();t++){
std::cout << " Difference "<<t<<" "<< res_chroma[t] - res_grid[t] << std::endl;
}*/
std::cout << "Finished test " << std::endl;
}
Chroma::finalize();
}