mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
GitHub
commit
202f025fc7
252
Grid/qcd/utils/BaryonUtils.h
Normal file
252
Grid/qcd/utils/BaryonUtils.h
Normal file
@ -0,0 +1,252 @@
|
||||
/*************************************************************************************
|
||||
|
||||
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>
|
||||
|
||||
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);
|
||||
|
||||
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 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 <class mobj, class robj>
|
||||
static void baryon_site(const mobj &D1,
|
||||
const mobj &D2,
|
||||
const mobj &D3,
|
||||
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_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 <class mobj, class robj>
|
||||
static void ContractBaryons_Sliced(const mobj &D1,
|
||||
const mobj &D2,
|
||||
const mobj &D3,
|
||||
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);
|
||||
};
|
||||
|
||||
template <class FImpl>
|
||||
constexpr int BaryonUtils<FImpl>::epsilon[6][3];
|
||||
template <class FImpl>
|
||||
constexpr Complex BaryonUtils<FImpl>::epsilon_sgn[6];
|
||||
|
||||
template <class FImpl>
|
||||
template <class mobj, class robj>
|
||||
void BaryonUtils<FImpl>::baryon_site(const mobj &D1,
|
||||
const mobj &D2,
|
||||
const mobj &D3,
|
||||
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)
|
||||
{
|
||||
|
||||
Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4)
|
||||
|
||||
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_right * D3;
|
||||
|
||||
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_left;
|
||||
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
|
||||
for (int beta_left=0; beta_left<Ns; beta_left++){
|
||||
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
|
||||
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,beta_left)(b_right,b_left);
|
||||
}}}
|
||||
}
|
||||
//This is the \delta_{456}^{231} part
|
||||
if (wick_contraction[1]){
|
||||
auto pD1g = pD1 * GammaB_left;
|
||||
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
|
||||
for (int beta_left=0; beta_left<Ns; beta_left++){
|
||||
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
|
||||
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3()(alpha_right,gamma_left)(b_right,c_left);
|
||||
}}}
|
||||
}
|
||||
//This is the \delta_{456}^{312} part
|
||||
if (wick_contraction[2]){
|
||||
auto gD3g = gD3 * GammaB_left;
|
||||
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
|
||||
for (int beta_left=0; beta_left<Ns; beta_left++){
|
||||
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
|
||||
result()()() += epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,beta_left)(c_right,b_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3g()(alpha_right,beta_left)(b_right,a_left);
|
||||
}}}
|
||||
}
|
||||
//This is the \delta_{456}^{132} part
|
||||
if (wick_contraction[3]){
|
||||
auto gD3g = gD3 * GammaB_left;
|
||||
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
|
||||
for (int beta_left=0; beta_left<Ns; beta_left++){
|
||||
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
|
||||
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,gamma_left)(c_right,c_left)*D2()(alpha_right,beta_left)(a_right,b_left)*gD3g()(alpha_right,beta_left)(b_right,a_left);
|
||||
}}}
|
||||
}
|
||||
//This is the \delta_{456}^{321} part
|
||||
if (wick_contraction[4]){
|
||||
auto D2g = D2 * GammaB_left;
|
||||
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
|
||||
for (int beta_left=0; beta_left<Ns; beta_left++){
|
||||
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
|
||||
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1()(gamma_left,beta_left)(c_right,b_left)*D2g()(alpha_right,beta_left)(a_right,a_left)*gD3()(alpha_right,gamma_left)(b_right,c_left);
|
||||
}}}
|
||||
}
|
||||
//This is the \delta_{456}^{213} part
|
||||
if (wick_contraction[5]){
|
||||
auto pD1g = pD1 * GammaB_left;
|
||||
for (int alpha_right=0; alpha_right<Ns; alpha_right++){
|
||||
for (int beta_left=0; beta_left<Ns; beta_left++){
|
||||
for (int gamma_left=0; gamma_left<Ns; gamma_left++){
|
||||
result()()() -= epsilon_sgn[ie_left] * epsilon_sgn[ie_right] * pD1g()(gamma_left,beta_left)(c_right,a_left)*D2()(alpha_right,gamma_left)(a_right,c_left)*gD3()(alpha_right,beta_left)(b_right,b_left);
|
||||
}}}
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
template<class FImpl>
|
||||
void BaryonUtils<FImpl>::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_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_left.Grid();
|
||||
|
||||
int wick_contraction[6];
|
||||
for (int ie=0; ie < 6 ; ie++)
|
||||
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_left.View();
|
||||
auto v2 = q2_left.View();
|
||||
auto v3 = q3_left.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_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result);
|
||||
vbaryon_corr[ss] = result;
|
||||
} );//end loop over lattice sites
|
||||
}
|
||||
template <class FImpl>
|
||||
template <class mobj, class robj>
|
||||
void BaryonUtils<FImpl>::ContractBaryons_Sliced(const mobj &D1,
|
||||
const mobj &D2,
|
||||
const mobj &D3,
|
||||
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_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_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_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result);
|
||||
}
|
||||
NAMESPACE_END(Grid);
|
||||
@ -7,7 +7,7 @@ Source file: Hadrons/Modules/MContraction/Baryon.hpp
|
||||
Copyright (C) 2015-2019
|
||||
|
||||
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
|
||||
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/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
#include <Grid/qcd/utils/BaryonUtils.h>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
@ -41,6 +42,9 @@ BEGIN_HADRONS_NAMESPACE
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MContraction)
|
||||
|
||||
typedef std::pair<Gamma::Algebra, Gamma::Algebra> GammaAB;
|
||||
typedef std::pair<GammaAB, GammaAB> GammaABPair;
|
||||
|
||||
class BaryonPar: Serializable
|
||||
{
|
||||
public:
|
||||
@ -48,6 +52,11 @@ public:
|
||||
std::string, q1,
|
||||
std::string, q2,
|
||||
std::string, q3,
|
||||
std::string, gammas,
|
||||
std::string, quarks,
|
||||
std::string, prefactors,
|
||||
std::string, parity,
|
||||
std::string, sink,
|
||||
std::string, output);
|
||||
};
|
||||
|
||||
@ -58,12 +67,21 @@ public:
|
||||
FERM_TYPE_ALIASES(FImpl1, 1);
|
||||
FERM_TYPE_ALIASES(FImpl2, 2);
|
||||
FERM_TYPE_ALIASES(FImpl3, 3);
|
||||
class Result: Serializable
|
||||
BASIC_TYPE_ALIASES(ScalarImplCR, Scalar);
|
||||
SINK_TYPE_ALIASES(Scalar);
|
||||
class Metadata: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
|
||||
std::vector<std::vector<std::vector<Complex>>>, corr);
|
||||
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);
|
||||
};
|
||||
typedef Correlator<Metadata> Result;
|
||||
public:
|
||||
// constructor
|
||||
TBaryon(const std::string name);
|
||||
@ -72,11 +90,14 @@ public:
|
||||
// dependency relation
|
||||
virtual std::vector<std::string> getInput(void);
|
||||
virtual std::vector<std::string> getOutput(void);
|
||||
virtual void parseGammaString(std::vector<GammaABPair> &gammaList);
|
||||
protected:
|
||||
// setup
|
||||
virtual void setup(void);
|
||||
// execution
|
||||
virtual void execute(void);
|
||||
// Which gamma algebra was specified
|
||||
Gamma::Algebra al;
|
||||
};
|
||||
|
||||
MODULE_REGISTER_TMP(Baryon, ARG(TBaryon<FIMPL, FIMPL, FIMPL>), MContraction);
|
||||
@ -94,7 +115,7 @@ TBaryon<FImpl1, FImpl2, FImpl3>::TBaryon(const std::string name)
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
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, par().q2, par().q3, par().sink};
|
||||
|
||||
return input;
|
||||
}
|
||||
@ -107,30 +128,199 @@ std::vector<std::string> TBaryon<FImpl1, FImpl2, FImpl3>::getOutput(void)
|
||||
return out;
|
||||
}
|
||||
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
void TBaryon<FImpl1, FImpl2,FImpl3>::parseGammaString(std::vector<GammaABPair> &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 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 MinusSigmaYT)");
|
||||
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.");
|
||||
|
||||
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;
|
||||
for (std::sregex_iterator i = gamma_begin; i != gamma_end; ++i) {
|
||||
nGamma++;
|
||||
}
|
||||
*/
|
||||
gammaList.resize(nGamma/4);
|
||||
std::vector<std::string> 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<Gamma::Algebra> gS1 = strToVec<Gamma::Algebra>(gS[4*i]);
|
||||
std::vector<Gamma::Algebra> gS2 = strToVec<Gamma::Algebra>(gS[4*i+1]);
|
||||
std::vector<Gamma::Algebra> gS3 = strToVec<Gamma::Algebra>(gS[4*i+2]);
|
||||
std::vector<Gamma::Algebra> gS4 = strToVec<Gamma::Algebra>(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];
|
||||
}
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
void TBaryon<FImpl1, FImpl2, FImpl3>::setup(void)
|
||||
{
|
||||
envTmpLat(LatticeComplex, "c");
|
||||
envTmpLat(LatticeComplex, "c2");
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl1, typename FImpl2, typename FImpl3>
|
||||
void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
|
||||
{
|
||||
LOG(Message) << "Computing baryon contractions '" << getName() << "' using"
|
||||
<< " quarks '" << par().q1 << "', '" << par().q2 << "', and '"
|
||||
<< par().q3 << "'" << std::endl;
|
||||
|
||||
auto &q1 = envGet(PropagatorField1, par().q1);
|
||||
auto &q2 = envGet(PropagatorField2, par().q2);
|
||||
auto &q3 = envGet(PropagatorField3, par().q2);
|
||||
std::vector<std::string> quarks = strToVec<std::string>(par().quarks);
|
||||
std::vector<double> prefactors = strToVec<double>(par().prefactors);
|
||||
int nQ=quarks.size();
|
||||
const int parity {par().parity.size()>0 ? std::stoi(par().parity) : 1};
|
||||
|
||||
std::vector<GammaABPair> 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.");
|
||||
|
||||
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 << 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);
|
||||
Result result;
|
||||
envGetTmp(LatticeComplex, c2);
|
||||
int nt = env().getDim(Tp);
|
||||
std::vector<TComplex> buf;
|
||||
TComplex cs;
|
||||
TComplex ch;
|
||||
|
||||
// FIXME: do contractions
|
||||
std::vector<Result> 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
|
||||
envHasType(SlicedPropagator3, par().q3))
|
||||
{
|
||||
auto &q1 = envGet(SlicedPropagator1, par().q1);
|
||||
auto &q2 = envGet(SlicedPropagator2, par().q2);
|
||||
auto &q3 = envGet(SlicedPropagator3, par().q3);
|
||||
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();
|
||||
for (int iQ1 = 0; iQ1 < nQ; iQ1++){
|
||||
for (int iQ2 = 0; iQ2 < nQ; iQ2++){
|
||||
BaryonUtils<FIMPL>::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;
|
||||
}
|
||||
}
|
||||
r.corr.push_back(TensorRemove(cs));
|
||||
}
|
||||
result.push_back(r);
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
auto &q1 = envGet(PropagatorField1, par().q1);
|
||||
auto &q2 = envGet(PropagatorField2, par().q2);
|
||||
auto &q3 = envGet(PropagatorField3, par().q3);
|
||||
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);
|
||||
|
||||
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<FIMPL>::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);
|
||||
}
|
||||
else if (ns == "MSink")
|
||||
{
|
||||
c=Zero();
|
||||
for (int iQ1 = 0; iQ1 < nQ; iQ1++){
|
||||
for (int iQ2 = 0; iQ2 < nQ; iQ2++){
|
||||
BaryonUtils<FIMPL>::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);
|
||||
}
|
||||
r.corr.clear();
|
||||
for (unsigned int t = 0; t < buf.size(); ++t)
|
||||
{
|
||||
r.corr.push_back(TensorRemove(buf[t]));
|
||||
}
|
||||
result.push_back(r);
|
||||
}
|
||||
}
|
||||
|
||||
saveResult(par().output, "baryon", result);
|
||||
|
||||
// saveResult(par().output, "meson", result);
|
||||
}
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
@ -46,6 +46,7 @@ int main(int argc, char *argv[])
|
||||
// run setup ///////////////////////////////////////////////////////////////
|
||||
Application application;
|
||||
std::vector<std::string> flavour = {"l", "s", "c1", "c2", "c3"};
|
||||
std::vector<std::string> flavour_baryon = {"l", "s", "a", "b", "c"}; //needs to be a single character
|
||||
std::vector<double> mass = {.01, .04, .2 , .25 , .3 };
|
||||
|
||||
// global parameters
|
||||
@ -134,6 +135,10 @@ 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.quarks = flavour_baryon[i] + flavour_baryon[j] + flavour_baryon[k];
|
||||
barPar.prefactors = "1.0";
|
||||
barPar.sink = "sink";
|
||||
application.createModule<MContraction::Baryon>(
|
||||
"baryon_pt_" + flavour[i] + flavour[j] + flavour[k], barPar);
|
||||
}
|
||||
|
||||
604
tests/qdpxx/Test_qdpxx_baryon.cc
Normal file
604
tests/qdpxx/Test_qdpxx_baryon.cc
Normal file
@ -0,0 +1,604 @@
|
||||
/*************************************************************************************
|
||||
|
||||
Grid physics library, www.github.com/paboyle/Grid
|
||||
|
||||
Source file: ./tests/qdpxx/Test_qdpxx_wilson.cc
|
||||
|
||||
Copyright (C) 2017
|
||||
|
||||
Author: Felix Erben <felix.erben@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
|
||||
|
||||
Grid::LatticeComplex c(UGrid);
|
||||
Grid::LatticeComplex c1(UGrid);
|
||||
|
||||
if(! baryon.compare("OmegaX")){
|
||||
BaryonUtils<Grid::WilsonImplR>::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<Grid::WilsonImplR>::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<Grid::WilsonImplR>::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<Grid::WilsonImplR>::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<Grid::WilsonImplR>::ContractBaryons(qS,qU,qD,G_A,G_B,G_A,G_B,"sud","sud",1,c1); //<ud>s
|
||||
c = 4.*c1;
|
||||
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qD,qU,qS,G_A,G_B,G_A,G_B,"dus","dus",1,c1); //<us>d
|
||||
c += c1;
|
||||
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qU,qD,qS,G_A,G_B,G_A,G_B,"uds","uds",1,c1); //<ds>u
|
||||
c += c1;
|
||||
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qD,qU,qS,G_A,G_B,G_A,G_B,"dus","sud",1,c1); //(sud)
|
||||
c += 2.*c1;
|
||||
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qU,qD,qS,G_A,G_B,G_A,G_B,"uds","sud",1,c1); //(sdu)
|
||||
c -= 2.*c1;
|
||||
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qS,qU,qD,G_A,G_B,G_A,G_B,"sud","dus",1,c1); //(dus)
|
||||
c += 2.*c1;
|
||||
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qU,qD,qS,G_A,G_B,G_A,G_B,"uds","dus",1,c1); //-(dsu)
|
||||
c -= c1;
|
||||
BaryonUtils<Grid::WilsonImplR>::ContractBaryons(qS,qU,qD,G_A,G_B,G_A,G_B,"sud","uds",1,c1); //(uds)
|
||||
c -= 2.*c1;
|
||||
BaryonUtils<Grid::WilsonImplR>::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 {
|
||||
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;
|
||||
Grid::Complex res_chroma_g;
|
||||
|
||||
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++){
|
||||
res_chroma_g = Grid::Complex(toDouble(real(res_chroma[t])), toDouble(imag(res_chroma[t])));
|
||||
std::cout << " Difference "<<t<<" "<< res_chroma_g - res_grid[t] << std::endl;
|
||||
}
|
||||
|
||||
std::cout << "Finished test " << std::endl;
|
||||
|
||||
}
|
||||
Chroma::finalize();
|
||||
}
|
||||
Loading…
Reference in New Issue
Block a user