portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-11-04 05:54:32 +00:00
Code Releases Activity

baryon field structure is now eigentensor - started on contractions for 2pt functions

Browse Source
This commit is contained in:
Felix Erben
2019-02-14 16:44:54 +00:00
parent 59c8cc1588
commit 886c895f81
1 changed files with 32 additions and 57 deletions
Download Patch File Download Diff File Expand all files Collapse all files

89
Hadrons/Modules/MDistil/BContraction.hpp
View File

@@ -18,6 +18,20 @@ BEGIN_HADRONS_NAMESPACE
******************************************************************************/
BEGIN_MODULE_NAMESPACE(MDistil)
// general baryon tensor set based on Eigen tensors and Grid-allocated memory
// Dimensions:
// 0 - ext - external field (momentum, EM field, ...)
// 1 - str - spin-color structure
// 2 - t - timeslice
// 3 - i - left distillation mode index
// 4 - j - middle distillation mode index
// 5 - k - left distillation mode index
// template <typename T>
// using BaryonTensorSet = Eigen::TensorMap<Eigen::Tensor<T, 6, Eigen::RowMajor>>;
class BContractionPar: Serializable
{
public:
@@ -141,6 +155,7 @@ void TBContraction<FImpl>::execute(void)
SpinVector tmp222;
SpinVector tmp111;
assert(parity == 1 || parity == -1);
std::vector<std::vector<int>> epsilon = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
@@ -159,61 +174,11 @@ void TBContraction<FImpl>::execute(void)
Gamma::Algebra::GammaYGamma5, // i gamma_4 C gamma_5 = i gamma_2 gamma_5
};
std::vector<Complex> factor23{{0.,-1.},{0.,1.},{0.,1.}};
if((0)){
for (int i1=0 ; i1 < N_1 ; i1++){
for (int i2=0 ; i2 < N_2 ; i2++){
for (int i3=0 ; i3 < N_3 ; i3++){
for (int imom=0 ; imom < Nmom ; imom++){
for (int t=0 ; t < Nt ; t++){
Bindex = i1 + N_1*(i2 + N_2*(i3 + N_3*(imom+Nmom*t)));
ExtractSliceLocal(tmp1,one[i1],0,t,3);
ExtractSliceLocal(tmp2,two[i2],0,t,3);
ExtractSliceLocal(tmp3,three[i3],0,t,3);
parallel_for (unsigned int sU = 0; sU < grid3d->oSites(); ++sU)
{
for (int ie=0 ; ie < 6 ; ie++){
// Why does peekColour not work????
for (int is=0 ; is < 4 ; is++){
tmp11s()(is)() = tmp11[sU]()(is)(epsilon[ie][0]);
tmp22s()(is)() = tmp22[sU]()(is)(epsilon[ie][1]);
tmp33s()(is)() = tmp33[sU]()(is)(epsilon[ie][2]);
}
tmp333 = Gamma(gamma23_[0])*tmp33s;
// this should be outerProduct??? Does not work.
for (int isl=0 ; isl < 4 ; isl++){
for (int isr=0 ; isr < 4 ; isr++){
diquark()(isl,isr)() = factor23[0]*tmp22s()(isl)()*tmp333()(isr)();
}
}
// Is there a way to compute gamma * SpinMatrix (left component)???
for (int isr=0 ; isr < 4 ; isr++){
for (int isl=0 ; isl < 4 ; isl++){
tmp222()(isl)() = diquark()(isl,isr)();
}
tmp111 = Gamma(gamma12_[0]) * tmp222;
for (int isl=0 ; isl < 4 ; isl++){
g_diquark()(isl,isr)() = tmp111()(isl)();
}
}
for (int is=0 ; is < 4 ; is++){
BField[Bindex]+=(double)epsilon_sgn[ie]*tmp11s()(is)()*g_diquark()(is,is)();
}
}
}
}
}
}
}
}
//BaryonTensorSet<Complex> BField3(storage,Nmom,4,Nt,N_1,N_2,N_3);
using BaryonTensorSet = Eigen::Tensor<Complex, 6, Eigen::RowMajor>;
BaryonTensorSet BField3(Nmom,4,Nt,N_1,N_2,N_3);
for (int t=0 ; t < Nt ; t++){
Bindex = 0 + N_1*(0 + N_2*(0 + N_3*(0+Nmom*t)));
std::cout << "BaryonField(t=" << t << ") = " << BField[Bindex] << std::endl;
}
} // end if 0
// ONLY THIS IS CORRECT?
std::vector<SpinVector> BField2(Nmom*Nt*N_1*N_2*N_3);
std::vector<SpinVector> BField2(Nmom*Nt*N_1*N_2*N_3);
Complex diquark2;
for (int i1=0 ; i1 < N_1 ; i1++){
for (int i2=0 ; i2 < N_2 ; i2++){
@@ -238,7 +203,10 @@ void TBContraction<FImpl>::execute(void)
tmp222 = g4*tmp111;
tmp111 = 0.5*(double)parity*(tmp111 + tmp222); // P_\pm * ...
diquark2 = factor23[0]*innerProduct(tmp22s,tmp333);
BField2[Bindex]+=(double)epsilon_sgn[ie]*tmp111*diquark2;
//BField2[Bindex]+=(double)epsilon_sgn[ie]*tmp111*diquark2;
for (int is=0 ; is < 4 ; is++){
BField3(imom,is,t,i1,i2,i3)+=(double)epsilon_sgn[ie]*tmp111()(is)()*diquark2;
}
}
}
}
@@ -248,10 +216,17 @@ void TBContraction<FImpl>::execute(void)
}
for (int is=0 ; is < 4 ; is++){
for (int t=0 ; t < Nt ; t++){
Bindex = 0 + N_1*(0 + N_2*(0 + N_3*(0+Nmom*t)));
std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField2[Bindex]()(is)() << std::endl;
// Bindex = 0 + N_1*(0 + N_2*(0 + N_3*(0+Nmom*t)));
// std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField2[Bindex]()(is)() << std::endl;
std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField3(0,is,t,0,0,0) << std::endl;
}
}
//Product ijk * ijk
// for ijk * jik: (4,5),(5,4),(6,6) z.b.
Eigen::array<Eigen::IndexPair<int>, 3> product_dims = { Eigen::IndexPair<int>(4, 4),Eigen::IndexPair<int>(5, 5) ,Eigen::IndexPair<int>(6, 6) };
// Whycan't I choose the dimension to be 3??? Want to sum over them, not save each element!
Eigen::Tensor<Complex,6,Eigen::RowMajor> C2 = BField3.contract(BField3,product_dims);
}
END_MODULE_NAMESPACE