diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index 4f65ca67..e34c4d72 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -53,8 +53,8 @@ public:
 					 const PropagatorField &q3_src,
 					 const Gamma GammaA,
 					 const Gamma GammaB,
-					 const char quarks_snk[],
-					 const char quarks_src[],
+					 const char * quarks_snk,
+					 const char * quarks_src,
 					 const int parity,
 					 ComplexField &baryon_corr);
  
@@ -67,12 +67,15 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
 						 const PropagatorField &q3_src,
 						 const Gamma GammaA,
 						 const Gamma GammaB,
-						 const char quarks_snk[],
-						 const char quarks_src[],
+						 const char * quarks_snk,
+						 const char * quarks_src,
 						 const int parity,
 						 ComplexField &baryon_corr)
 {
-
+  std::cout << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] <<  std::endl;
+    std::cout << "GammaA " << (GammaA.g) <<  std::endl;
+    std::cout << "GammaB " << (GammaB.g) <<  std::endl;
+ 
   assert(parity==1 || parity == -1 && "Parity must be +1 or -1");
 
   GridBase *grid = q1_src.Grid();
@@ -88,12 +91,16 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
       wick_contraction[ie]=1;
 
   typedef typename ComplexField::vector_object vobj;
-  LatticeView<vobj> vbaryon_corr{ baryon_corr };
-  accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+  auto vbaryon_corr= baryon_corr.View();
+  auto v1 = q1_src.View();
+  auto v2 = q2_src.View();
+  auto v3 = q3_src.View();
+
+      std::cout << "wick contract "  << wick_contraction <<  std::endl;
+    int count=0;
+ // accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+  thread_for(ss,grid->oSites(),{
 
-    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];
@@ -103,7 +110,11 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
     auto pD1 = 0.5* (gD1a + (double)parity * gD1b);
     auto gD3 = GammaB * D3;
 
-    vobj result{ 0 };
+    vobj result=Zero();
+    if (count<10){
+      std::cout << "outside epsilon "  << count <<  std::endl;
+      count++;
+    }
     
     for (int ie_src=0; ie_src < 6 ; ie_src++){
       int a_src = epsilon[ie_src][0]; //a
@@ -132,6 +143,7 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_src,
           }}}
         }	  
         //This is the \delta_{456}^{312} part
+       // for(int test=0;test<3;test++){
 	if (wick_contraction[2]){
           auto gD3g = gD3 * GammaB;
 	  for (int alpha_snk=0; alpha_snk<Ns; alpha_snk++){
@@ -149,6 +161,7 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_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;
diff --git a/Hadrons/Modules/MContraction/Baryon.hpp b/Hadrons/Modules/MContraction/Baryon.hpp
index fc4cc8fc..8f831c9b 100644
--- a/Hadrons/Modules/MContraction/Baryon.hpp
+++ b/Hadrons/Modules/MContraction/Baryon.hpp
@@ -34,6 +34,7 @@ See the full license in the file "LICENSE" in the top level distribution directo
 #include <Hadrons/Module.hpp>
 #include <Hadrons/ModuleFactory.hpp>
 #include <Grid/qcd/utils/BaryonUtils.h>
+#include <Grid/qcd/utils/A2Autils.h>
 
 BEGIN_HADRONS_NAMESPACE
 
@@ -134,24 +135,49 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
                  << " ) and parity " << par().parity << "." << std::endl;
     
     auto       &q1_src = envGet(PropagatorField1, par().q1_src);
+    LOG(Message) << "1" << std::endl;
     auto       &q2_src = envGet(PropagatorField2, par().q2_src);
     auto       &q3_src = envGet(PropagatorField3, par().q3_src);
     envGetTmp(LatticeComplex, c);
     Result     result;
     int nt = env().getDim(Tp);
     result.corr.resize(nt);
+    LOG(Message) << "2" << std::endl;
     std::vector<Gamma::Algebra> ggA = strToVec<Gamma::Algebra>(par().GammaA);
-    Gamma GammaA = ggA[0];
+    LOG(Message) << "3" << std::endl;
+    Gamma GammaA(ggA[0]);
+    LOG(Message) << "4" << std::endl;
     std::vector<Gamma::Algebra> ggB = strToVec<Gamma::Algebra>(par().GammaB);
-    Gamma GammaB = ggB[0];
+    Gamma GammaB(ggB[0]);
     std::vector<TComplex> buf;
     const int  parity {par().parity};
+    LOG(Message) << "5" << std::endl;
     const char * quarks_snk{par().quarks_snk.c_str()};
+    LOG(Message) << "6" << std::endl;
     const char * quarks_src{par().quarks_src.c_str()};
+    LOG(Message) << "quarks_snk " << quarks_snk[0] << quarks_snk[1] << quarks_snk[2] <<  std::endl;
+    LOG(Message) << "GammaA " << (GammaA.g) <<  std::endl;
+    LOG(Message) << "GammaB " << (GammaB.g) <<  std::endl;
 
-    BaryonUtils<FIMPL>::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c);
+    GridCartesian *Ugrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(),GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
+    LatticePropagator q(Ugrid);
+    GridParallelRNG RNG4(Ugrid);
+    gaussian(RNG4,q);
+    Gamma gA = Gamma(Gamma::Algebra::Identity);
+    Gamma gB = Gamma(Gamma::Algebra::Gamma5);
+    int p=1;
+    char * om = "sss";
+    LatticeComplex c2(Ugrid);
+
+    //BaryonUtils<FIMPL>::ContractBaryons(q1_src,q2_src,q3_src,GammaA,GammaB,quarks_snk,quarks_src,parity,c);
+    BaryonUtils<FIMPL>::ContractBaryons(q,q,q,gA,gB,om,om,p,c);
+    std::vector<Gamma> GA={gA};
+    std::vector<Gamma> GB={gB};
+    //A2Autils<FIMPL>::ContractFourQuarkColourMix(q,q,GA,GB,c,c2);
+    LOG(Message) << "survived ContractBaryons" << std::endl;
 
     sliceSum(c,buf,Tp);
+    LOG(Message) << "survived sliceSum" << std::endl;
 
     for (unsigned int t = 0; t < buf.size(); ++t)
     {