Merge branch 'feature/distil' of github.com:mmphys/Grid into feature/distil

Hadrons/Modules/MDistil/BContraction.hpp

@@ -26,6 +26,7 @@ public:
 		                    std::string, two,
 		                    std::string, three,
 		                    std::string, output,
+				    int, parity,
                                     std::vector<std::string>, mom);
 };
 
@@ -104,6 +105,8 @@ void TBContraction<FImpl>::execute(void)
     int N_2     = two.size();
     int N_3     = three.size();
 
+    int parity = par().parity;
+
     LOG(Message) << "Computing distillation baryon fields" << std::endl;
     LOG(Message) << "One: '" << par().one << "' Two: '" << par().two  << "' Three: '" << par().three << "'" << std::endl;
     LOG(Message) << "Momenta:" << std::endl;
@@ -138,9 +141,13 @@ 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}};
     std::vector<int> epsilon_sgn = {1,1,1,-1,-1,-1};
 
+    Gamma g4(Gamma::Algebra::GammaT);
+
     gamma12_ = {
        Gamma::Algebra::Identity, // I
        Gamma::Algebra::Gamma5, // gamma_5
@@ -160,6 +167,8 @@ void TBContraction<FImpl>::execute(void)
             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++){
@@ -186,7 +195,6 @@ void TBContraction<FImpl>::execute(void)
 		      g_diquark()(isl,isr)() = tmp111()(isl)();
 		    }
 		  }
-		  // Really only the trace? Should check baryons again! laph paper lists c_{alpha,beta,gamma}, gattringer-lang two gamma matrices.
                   for (int is=0 ; is < 4 ; is++){
                     BField[Bindex]+=(double)epsilon_sgn[ie]*tmp11s()(is)()*g_diquark()(is,is)();
   		  }
@@ -214,6 +222,8 @@ void TBContraction<FImpl>::execute(void)
             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++){
@@ -224,8 +234,11 @@ void TBContraction<FImpl>::execute(void)
 	            tmp33s()(is)() = tmp33[sU]()(is)(epsilon[ie][2]);
 		  }
 		  tmp333 = Gamma(gamma23_[0])*tmp33s;
+		  tmp111 = Gamma(gamma12_[0])*tmp11s;
+		  tmp222 = g4*tmp111;
+		  tmp111 = 0.5*(double)parity*(tmp111 + tmp222); // P_\pm * ...
                   diquark2 = factor23[0]*innerProduct(tmp22s,tmp333);
-                  BField2[Bindex]+=(double)epsilon_sgn[ie]*tmp11s*diquark2;
+                  BField2[Bindex]+=(double)epsilon_sgn[ie]*tmp111*diquark2;
 		}
   	      }
             }

tests/hadrons/Test_hadrons_distil.cc

@@ -225,19 +225,36 @@ void test_MesonField(Application &application)
   application.createModule<MContraction::A2AMesonField>("DistilMesonField",A2AMesonFieldPar);
 }
 /////////////////////////////////////////////////////////////
-// BaryonFields
+// BaryonFields - phiphiphi
 /////////////////////////////////////////////////////////////
 
-void test_BaryonField(Application &application)
+void test_BaryonFieldPhi(Application &application)
 {
   // DistilVectors parameters
   MDistil::BContraction::Par BContractionPar;
   BContractionPar.one="DistilVecs_phi";
   BContractionPar.two="DistilVecs_phi";
   BContractionPar.three="DistilVecs_phi";
-  BContractionPar.output="BaryonField";
+  BContractionPar.output="BaryonFieldPhi";
+  BContractionPar.parity=1;
   BContractionPar.mom={"0 0 0"};
-  application.createModule<MDistil::BContraction>("BaryonField",BContractionPar);
+  application.createModule<MDistil::BContraction>("BaryonFieldPhi",BContractionPar);
+}
+/////////////////////////////////////////////////////////////
+// BaryonFields - rhorhorho
+/////////////////////////////////////////////////////////////
+
+void test_BaryonFieldRho(Application &application)
+{
+  // DistilVectors parameters
+  MDistil::BContraction::Par BContractionPar;
+  BContractionPar.one="DistilVecs_rho";
+  BContractionPar.two="DistilVecs_rho";
+  BContractionPar.three="DistilVecs_rho";
+  BContractionPar.output="BaryonFieldRho";
+  BContractionPar.parity=1;
+  BContractionPar.mom={"0 0 0"};
+  application.createModule<MDistil::BContraction>("BaryonFieldRho",BContractionPar);
 }
 
 bool bNumber( int &ri, const char * & pstr, bool bGobbleWhiteSpace = true )
@@ -571,7 +588,8 @@ int main(int argc, char *argv[])
       test_LapEvec( application );
       test_Perambulators( application );
       test_DistilVectors( application );
-      test_BaryonField( application );
+      test_BaryonFieldPhi( application );
+      test_BaryonFieldRho( application );
       break;
   }
   LOG(Message) << "====== XML creation for test " << iTestNum << " complete ======" << std::endl;