From 886c895f81f4ebacf7317289a4cb0a43ff2d3f39 Mon Sep 17 00:00:00 2001 From: Felix Erben Date: Thu, 14 Feb 2019 16:44:54 +0000 Subject: [PATCH] baryon field structure is now eigentensor - started on contractions for 2pt functions --- Hadrons/Modules/MDistil/BContraction.hpp | 89 +++++++++--------------- 1 file changed, 32 insertions(+), 57 deletions(-) diff --git a/Hadrons/Modules/MDistil/BContraction.hpp b/Hadrons/Modules/MDistil/BContraction.hpp index 4401e8cd..71f6474c 100644 --- a/Hadrons/Modules/MDistil/BContraction.hpp +++ b/Hadrons/Modules/MDistil/BContraction.hpp @@ -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 + // using BaryonTensorSet = Eigen::TensorMap>; + + class BContractionPar: Serializable { public: @@ -141,6 +155,7 @@ void TBContraction::execute(void) SpinVector tmp222; SpinVector tmp111; + assert(parity == 1 || parity == -1); std::vector> 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::execute(void) Gamma::Algebra::GammaYGamma5, // i gamma_4 C gamma_5 = i gamma_2 gamma_5 }; std::vector 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 BField3(storage,Nmom,4,Nt,N_1,N_2,N_3); + using BaryonTensorSet = Eigen::Tensor; + 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 BField2(Nmom*Nt*N_1*N_2*N_3); + std::vector 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::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::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, 3> product_dims = { Eigen::IndexPair(4, 4),Eigen::IndexPair(5, 5) ,Eigen::IndexPair(6, 6) }; + // Whycan't I choose the dimension to be 3??? Want to sum over them, not save each element! + Eigen::Tensor C2 = BField3.contract(BField3,product_dims); } END_MODULE_NAMESPACE