From 435653490e39ba50d5bcd15bc4426551013039af Mon Sep 17 00:00:00 2001
From: Felix Erben <dc-erbe1@tesseract-login1.ib0.sgi.cluster.dirac.ed.ac.uk>
Date: Fri, 17 May 2019 10:50:15 +0100
Subject: [PATCH] fixed contraction issue

---
 Grid/qcd/utils/A2Autils.h             | 10 +++++++++-
 Hadrons/Modules/MDistil/Baryon2pt.hpp | 24 ++++++++++--------------
 2 files changed, 19 insertions(+), 15 deletions(-)

diff --git a/Grid/qcd/utils/A2Autils.h b/Grid/qcd/utils/A2Autils.h
index 9e2bd7e2..81c53800 100644
--- a/Grid/qcd/utils/A2Autils.h
+++ b/Grid/qcd/utils/A2Autils.h
@@ -183,7 +183,8 @@ void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat,
 
        	  for(int j=0;j<twoBlock;j++){
 
-	    auto v2 = conjugate(two[j]._odata[ss]);
+	    //auto v2 = conjugate(two[j]._odata[ss]);
+	    auto v2 = two[j]._odata[ss];
             // C = i gamma_2 gamma_4 => C gamma_5 = - i gamma_1 gamma_3
 	    //auto v2g = v2*Gamma(Gamma::Algebra::SigmaXZ);
             //auto v2g=v2;  
@@ -209,6 +210,13 @@ void A2Autils<FImpl>::NucleonFieldMom(Eigen::Tensor<ComplexD,6> &mat,
                   +           pv1()(s1)(2) * v2()(s2)(0) * gv3()(s2)(1)    
                   -           pv1()(s1)(2) * v2()(s2)(1) * gv3()(s2)(0);    
               }}
+
+              /*if (i+j+k == 0) {
+                Serializable::WriteMember(std::cout, pv1);
+                Serializable::WriteMember(std::cout, v2);
+                Serializable::WriteMember(std::cout, gv3);
+                Serializable::WriteMember(std::cout, vv);
+              }*/
 	    
 	      // After getting the sitewise product do the mom phase loop
               int base = Nmom*i+Nmom*oneBlock*j+Nmom*oneBlock*twoBlock*k+Nmom*oneBlock*twoBlock*threeBlock*r;
diff --git a/Hadrons/Modules/MDistil/Baryon2pt.hpp b/Hadrons/Modules/MDistil/Baryon2pt.hpp
index 02ac892a..d4420ecc 100644
--- a/Hadrons/Modules/MDistil/Baryon2pt.hpp
+++ b/Hadrons/Modules/MDistil/Baryon2pt.hpp
@@ -186,20 +186,16 @@ void TBaryon2pt<FImpl>::execute(void)
       Eigen::array<Eigen::IndexPair<int>, 3> product_dims = { Eigen::IndexPair<int>(0,epsilon[pairs[ipair]][0]),Eigen::IndexPair<int>(1,epsilon[pairs[ipair]][1]) ,Eigen::IndexPair<int>(2,epsilon[pairs[ipair]][2])  };
       for (int imom=0 ; imom < Nmom ; imom++){
         std::cout << imom << std::endl;
-        Eigen::Tensor<Complex,6> B6L = BFieldL.BField.chip(imom,0);
-        Eigen::Tensor<Complex,6> B6R = BFieldR.BField.chip(imom,0);
-        for (int ig=0 ; ig < Ngamma ; ig++){
-          Eigen::Tensor<Complex,5> B5L = B6L.chip(ig,0);
-          Eigen::Tensor<Complex,5> B5R = B6R.chip(ig,0);
-          for (int t=0 ; t < Nt ; t++){
-            Eigen::Tensor<Complex,4> B4L = B5L.chip(t,0);
-            Eigen::Tensor<Complex,4> B4R = B5R.chip(tsrc,0);
-            for (int is=0 ; is < 4 ; is++){
-              Eigen::Tensor<Complex,3> B3L = B4L.chip(is,0);
-              Eigen::Tensor<Complex,3> B3R = B4R.chip(is,0);
-              Eigen::Tensor<Complex,0> C2 = B3L.contract(B3R,product_dims);
-              corr(imom,t) += static_cast<Real>(epsilon_sgn[pairs[ipair]])*C2(0);
-            }
+        Eigen::Tensor<Complex,5> B5L = BFieldL.BField.chip(imom,0);
+        Eigen::Tensor<Complex,5> B5R = BFieldR.BField.chip(imom,0);
+        for (int t=0 ; t < Nt ; t++){
+          Eigen::Tensor<Complex,4> B4L = B5L.chip(t,0);
+          Eigen::Tensor<Complex,4> B4R = B5R.chip(tsrc,0);
+          for (int is=0 ; is < 4 ; is++){
+            Eigen::Tensor<Complex,3> B3L = B4L.chip(is,0);
+            Eigen::Tensor<Complex,3> B3R = B4R.chip(is,0);
+            Eigen::Tensor<Complex,0> C2 = B3L.contract(B3R,product_dims);
+            corr(imom,t) += static_cast<Real>(epsilon_sgn[pairs[ipair]])*C2(0);
           }
         }
       }