From ee11678b1fe595f5442ff800a655929c4e03a31b Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Tue, 3 Nov 2020 12:41:35 +0000
Subject: [PATCH 01/10] added Xi-to-Sigma rare decays

---
 Grid/qcd/utils/BaryonUtils.h | 169 +++++++++++++++++++++++++++++++++++
 1 file changed, 169 insertions(+)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index 15516b56..e6b52a43 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -1292,4 +1292,173 @@ void BaryonUtils<FImpl>::SigmaToNucleonNonEye(const PropagatorField &qq_ti,
   }  );//end loop over lattice sites
 }
 
+
+/***********************************************************************
+ * The following code is for Xi -> Sigma rare hypeon decays            *
+ **********************************************************************/
+
+/* Dq_loop is a quark line from t_H to t_H
+ * Dd_spec is a quark line from t_i to t_f
+ * Ds_spec is a quark line from t_i to t_f
+ * Dd_tf is a quark line from t_f to t_H
+ * Ds_ti is a quark line from t_i to t_H */
+template <class FImpl>
+template <class mobj, class mobj2, class robj>
+void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
+						 const mobj2 &Dd_spec,
+						 const mobj2 &Ds_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_xi,
+		                 		 const Gamma GammaB_sigma,
+						 robj &result)
+{
+
+  Gamma g5(Gamma::Algebra::Gamma5); 
+
+  auto DdG = Dd_spec * GammaB_sigma;
+  auto GDs = GammaB_xi * Ds_spec;
+  // Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
+  auto DsGDd = Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
+  // DsGDd * GammaB
+  auto DsGDdG = DsGDd * GammaB_sigma;
+  // GammaB * DsGDd
+  auto GDsGDd = GammaB_xi * DsGDd;
+  // GammaB * DsGDd * GammaB
+  auto GDsGDdG = GDsGDd * GammaB_sigma;
+  // \gamma_\mu^L * Dq_loop 
+  auto trGDq = trace(GammaH * Dq_loop);
+
+  for (int ie_s=0; ie_s < 6 ; ie_s++){
+    int a_s = epsilon[ie_s][0]; //a
+    int b_s = epsilon[ie_s][1]; //b
+    int c_s = epsilon[ie_s][2]; //c
+    for (int ie_x=0; ie_x < 6 ; ie_x++){
+      int a_x = epsilon[ie_x][0]; //a'
+      int b_x = epsilon[ie_x][1]; //b'
+      int c_x = epsilon[ie_x][2]; //c'
+      auto ee_GD = epsilon_sgn[ie_s] * epsilon_sgn[ie_x] * trGDq;
+      for (int alpha_x=0; alpha_x<Ns; alpha_x++){
+      for (int beta_s=0; beta_s<Ns; beta_s++){
+        auto GDsGDdG_ab_ba = GDsGDd()(alpha_x,beta_s)(b_x,a_s);
+        auto Ds_ab_ab = Ds_spec()(alpha_x,beta_s)(a_x,b_s);
+        auto DsGDdG_ab_aa = DsGDd()(alpha_x,beta_s)(a_x,a_s);
+        auto GDs_ab_bb = GDs()(alpha_x,beta_s)(b_x,b_s);
+        for (int gamma_x=0; gamma_x<Ns; gamma_x++){
+          auto DdG_cb_ca = DdG()(gamma_x,beta_s)(c_x,a_s);
+          for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+            result()(gamma_x,gamma_s)() -= ee_GD * Dd_spec()(gamma_x,gamma_s)(c_x,c_s) * GDsGDdG_ab_ba * Ds_ab_ab;
+            result()(gamma_x,gamma_s)() += ee_GD * DdG_cb_ca * GDsGDd()(alpha_x,gamma_s)(b_x,c_s) * Ds_ab_ab;
+            result()(gamma_x,gamma_s)() += ee_GD * Dd_spec()(gamma_x,gamma_s)(c_x,c_s) * DsGDdG_ab_aa * GDs_ab_bb;
+            result()(gamma_x,gamma_s)() -= ee_GD * DdG_cb_ca * DsGDd()(alpha_x,gamma_s)(a_x,c_s) * GDs_ab_bb;
+          }
+	}
+      }}
+    }
+  }
+}
+
+
+//Equivalent to "One-trace"
+/* Dq_loop is a quark line from t_H to t_H
+ * Dd_spec is a quark line from t_i to t_f
+ * Ds_spec is a quark line from t_i to t_f
+ * Dd_tf is a quark line from t_f to t_H
+ * Ds_ti is a quark line from t_i to t_H */
+template <class FImpl>
+template <class mobj, class mobj2, class robj>
+void BaryonUtils<FImpl>::XiToSigmaQ2EyeSite(const mobj &Dq_loop,
+						 const mobj2 &Dd_spec,
+						 const mobj2 &Ds_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_xi,
+		                 		 const Gamma GammaB_sigma,
+						 robj &result)
+{
+
+  Gamma g5(Gamma::Algebra::Gamma5); 
+
+  auto DdG = Dd_spec * GammaB_sigma;
+  auto GDs = GammaB_xi * Ds_spec;
+  // Ds * \gamma_\mu^L * Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
+  auto DsGDqGDd = Ds_ti * Gamma_H * Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5;
+  // DsGDd * GammaB
+  auto DsGDqGDdG = DsGDqGDd * GammaB_sigma;
+  // GammaB * DsGDd
+  auto GDsGDqGDd = GammaB_xi * DsGDqGDd;
+  // GammaB * DsGDd * GammaB
+  auto GDsGDqGDdG = GDsGDqGDd * GammaB_sigma;
+
+  for (int ie_s=0; ie_s < 6 ; ie_s++){
+    int a_s = epsilon[ie_s][0]; //a
+    int b_s = epsilon[ie_s][1]; //b
+    int c_s = epsilon[ie_s][2]; //c
+    for (int ie_x=0; ie_x < 6 ; ie_x++){
+      int a_x = epsilon[ie_x][0]; //a'
+      int b_x = epsilon[ie_x][1]; //b'
+      int c_x = epsilon[ie_x][2]; //c'
+      auto ee = epsilon_sgn[ie_s] * epsilon_sgn[ie_x];
+      for (int alpha_x=0; alpha_x<Ns; alpha_x++){
+      for (int beta_s=0; beta_s<Ns; beta_s++){
+        auto GDsGDqGDdG_ab_ba = GDsGDqGDdG()(alpha_x,beta_s)(b_x,a_s);
+        auto Ds_ab_ab = Ds_spec()(alpha_x,beta_s)(a_x,b_s);
+        auto DsGDqGDdG_ab_aa = GDsGDqGDdG()(alpha_x,beta_s)(a_x,a_s);
+        auto GDs_ab_bb = GDs()(alpha_x,beta_s)(b_x,b_s);
+          for (int gamma_x=0; gamma_x<Ns; gamma_x++){
+            auto DdG_cb_ca = DdG()(gamma_x, beta_s)(c_x,a_s);
+          for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+            result()(gamma_x,gamma_s)() -= ee * Dd_spec()(gamma_x,gamma_s)(c_x,c_s) * GDsGDqGDdG_ab_ba * Ds_ab_ab;
+            result()(gamma_x,gamma_s)() += ee * DdG_cb_ca * GDsGDqGDd()(alpha_x,gamma_s)(b_x,c_s) * Ds_ab_ab;
+            result()(gamma_x,gamma_s)() += ee * Dd_spec()(gamma_x,gamma_s)(c_x,c_s) * DsGDqGDdG_ab_aa * GDs_ab_bb;
+            result()(gamma_x,gamma_s)() -= ee * DdG_cb_ca * DsGDqGDd()(alpha_x,gamma_s)(a_x,c_s) * GDs_ab_bb;
+          }}
+      }}
+    }
+  }
+}
+
+template<class FImpl>
+template <class mobj>
+void BaryonUtils<FImpl>::XiToSigmaEye(const PropagatorField &qq_loop,
+						 const mobj &Dd_spec,
+						 const mobj &Ds_spec,
+						 const PropagatorField &qd_tf,
+						 const PropagatorField &qs_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_xi,
+		                 		 const Gamma GammaB_sigma,
+						 const std::string op,
+						 SpinMatrixField &stn_corr)
+{
+
+  assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
+  assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
+
+  GridBase *grid = qs_ti.Grid();
+
+  autoView( vcorr, stn_corr, CpuWrite);
+  autoView( vq_loop , qq_loop, CpuRead);
+  autoView( vd_tf , qd_tf, CpuRead);
+  autoView( vs_ti , qs_ti, CpuRead);
+
+  accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+    auto Dq_loop = vq_loop[ss];
+    auto Dd_tf = vd_tf[ss];
+    auto Ds_ti = vs_ti[ss];
+    sobj result=Zero();
+    if(op == "Q1"){
+      XiToSigmaQ1EyeSite(Dq_loop,Dd_spec,Ds_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
+    } else if(op == "Q2"){
+      XiToSigmaQ2EyeSite(Dq_loop,Dd_spec,Ds_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
+    } else {
+      assert(0 && "Weak Operator not correctly specified");
+    }
+      vcorr[ss] = result; 
+  }  );//end loop over lattice sites
+}
+
+
 NAMESPACE_END(Grid);

From a3de7026c8d091b4d42eba5e5cb8e57419e0963e Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Tue, 3 Nov 2020 12:51:50 +0000
Subject: [PATCH 02/10] bugfix

---
 Grid/qcd/utils/BaryonUtils.h | 31 +++++++++++++++++++++++++++++++
 1 file changed, 31 insertions(+)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index e6b52a43..8e6bf722 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -207,6 +207,26 @@ public:
 				                 const Gamma GammaB_sigma,
 		                 		 const Gamma GammaB_nucl,
 						 robj &result);
+  template <class mobj, class mobj2, class robj>
+  static void XiToSigmaQ1EyeSite(const mobj &Dq_loop,
+						 const mobj2 &Dd_spec,
+						 const mobj2 &Ds_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_sigma,
+		                 		 const Gamma GammaB_nucl,
+						 robj &result);
+  template <class mobj, class mobj2, class robj>
+  static void XiToSigmaQ2EyeSite(const mobj &Dq_loop,
+						 const mobj2 &Dd_spec,
+						 const mobj2 &Ds_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_sigma,
+		                 		 const Gamma GammaB_nucl,
+						 robj &result);
   public:
   template <class mobj>
   static void SigmaToNucleonEye(const PropagatorField &qq_loop,
@@ -229,6 +249,17 @@ public:
 				 const Gamma GammaB_nucl,
 		                 const std::string op,
 				 SpinMatrixField &stn_corr);
+  template <class mobj>
+  static void XiToSigmaEye(const PropagatorField &qq_loop,
+				 const mobj &Dd_spec,
+				 const mobj &Ds_spec,
+				 const PropagatorField &qd_tf,
+				 const PropagatorField &qs_ti,
+				 const Gamma Gamma_H,
+				 const Gamma GammaB_sigma,
+				 const Gamma GammaB_nucl,
+		                 const std::string op,
+				 SpinMatrixField &stn_corr);
 };
 
 template <class FImpl> 

From 67023c334b16938c36af9f41647ba42891008e03 Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Tue, 3 Nov 2020 13:07:37 +0000
Subject: [PATCH 03/10] bugfix

---
 Grid/qcd/utils/BaryonUtils.h | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index 8e6bf722..a281a669 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -1359,7 +1359,7 @@ void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
   // GammaB * DsGDd * GammaB
   auto GDsGDdG = GDsGDd * GammaB_sigma;
   // \gamma_\mu^L * Dq_loop 
-  auto trGDq = trace(GammaH * Dq_loop);
+  auto trGDq = trace(Gamma_H * Dq_loop);
 
   for (int ie_s=0; ie_s < 6 ; ie_s++){
     int a_s = epsilon[ie_s][0]; //a

From 4014dfd5b9a79223023ab2bee3789a7efde45945 Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Tue, 3 Nov 2020 16:13:08 +0000
Subject: [PATCH 04/10] first tested version

---
 Grid/qcd/utils/BaryonUtils.h | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index a281a669..f5bc1480 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -1359,7 +1359,7 @@ void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
   // GammaB * DsGDd * GammaB
   auto GDsGDdG = GDsGDd * GammaB_sigma;
   // \gamma_\mu^L * Dq_loop 
-  auto trGDq = trace(Gamma_H * Dq_loop);
+  auto trGDq = TensorRemove(trace(Gamma_H * Dq_loop)); 
 
   for (int ie_s=0; ie_s < 6 ; ie_s++){
     int a_s = epsilon[ie_s][0]; //a
@@ -1369,7 +1369,7 @@ void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
       int a_x = epsilon[ie_x][0]; //a'
       int b_x = epsilon[ie_x][1]; //b'
       int c_x = epsilon[ie_x][2]; //c'
-      auto ee_GD = epsilon_sgn[ie_s] * epsilon_sgn[ie_x] * trGDq;
+      auto ee_GD = epsilon_sgn[ie_s] * epsilon_sgn[ie_x] * trGDq; 
       for (int alpha_x=0; alpha_x<Ns; alpha_x++){
       for (int beta_s=0; beta_s<Ns; beta_s++){
         auto GDsGDdG_ab_ba = GDsGDd()(alpha_x,beta_s)(b_x,a_s);

From 9bae6b889a94024b60c39bc2fb993eed7456e1b9 Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Tue, 3 Nov 2020 20:03:09 +0000
Subject: [PATCH 05/10] speedup in Sigma-to-nucleon

---
 Grid/qcd/utils/BaryonUtils.h | 193 ++++++++++++++++++++++++++++++++++-
 1 file changed, 189 insertions(+), 4 deletions(-)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index f5bc1480..c0999a4a 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -1027,7 +1027,7 @@ void BaryonUtils<FImpl>::BaryonGamma3pt(
  * Du_spec is a quark line from t_i to t_f
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
-template <class FImpl>
+/*template <class FImpl>
 template <class mobj, class mobj2, class robj>
 void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Du_spec,
@@ -1071,6 +1071,50 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
       }}
     }
   }
+}*/
+template <class FImpl>
+template <class mobj, class mobj2, class robj>
+void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
+						 const mobj2 &Du_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_sigma,
+		                 		 const Gamma GammaB_nucl,
+						 robj &result)
+{
+
+  Gamma g5(Gamma::Algebra::Gamma5); 
+
+  auto DuG = Du_spec * GammaB_nucl;
+  // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
+  auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
+  // Dq_loop * \gamma_\mu^L
+  auto trDqG = TensorRemove(trace(Dq_loop * Gamma_H));
+
+  for (int ie_n=0; ie_n < 6 ; ie_n++){
+    int a_n = epsilon[ie_n][0]; //a
+    int b_n = epsilon[ie_n][1]; //b
+    int c_n = epsilon[ie_n][2]; //c
+    for (int ie_s=0; ie_s < 6 ; ie_s++){
+      int a_s = epsilon[ie_s][0]; //a'
+      int b_s = epsilon[ie_s][1]; //b'
+      int c_s = epsilon[ie_s][2]; //c'
+      for (int alpha_s=0; alpha_s<Ns; alpha_s++){
+      for (int beta_n=0; beta_n<Ns; beta_n++){
+        auto GDsGDd_ab_bb = GDsGDd()(alpha_s,beta_n)(b_s,b_n);
+        auto ee_GDGDDG = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * trDqG;
+	auto DuG_ab_aa = DuG()(alpha_s, beta_n)(a_s,a_n);
+        for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+	  auto DuG_cb_ca = DuG()(gamma_s, beta_n)(c_s,a_n);
+          for (int gamma_n=0; gamma_n<Ns; gamma_n++){
+            result()(gamma_s,gamma_n)() += ee_GDGDDG * DuG_ab_aa * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
+            result()(gamma_s,gamma_n)() -= ee_GDGDDG * DuG_cb_ca * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
+	  }
+	}
+      }}
+    }
+  }
 }
 
 /* Du_ti is a quark line from t_i to t_H
@@ -1078,7 +1122,7 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
  * Du_spec is a quark line from t_i to t_f
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
-template <class FImpl>
+/*template <class FImpl>
 template <class mobj, class mobj2, class robj>
 void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
 						 const mobj &Du_tf,
@@ -1130,6 +1174,55 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
       }}
     }
   }
+}*/
+template <class FImpl>
+template <class mobj, class mobj2, class robj>
+void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
+						 const mobj &Du_tf,
+						 const mobj2 &Du_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_sigma,
+		                 		 const Gamma GammaB_nucl,
+						 robj &result)
+{
+
+  Gamma g5(Gamma::Algebra::Gamma5); 
+
+  auto DuG = Du_spec * GammaB_nucl;
+  // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
+  auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
+  // Du_ti * \gamma_\mu^L * adj(Du_tf)
+  auto DuGHDu = Du_ti * Gamma_H * g5 * adj(Du_tf) * g5;
+  auto DuGHDuG = DuGHDu * GammaB_nucl;
+
+  for (int ie_n=0; ie_n < 6 ; ie_n++){
+    int a_n = epsilon[ie_n][0]; //a
+    int b_n = epsilon[ie_n][1]; //b
+    int c_n = epsilon[ie_n][2]; //c
+    for (int ie_s=0; ie_s < 6 ; ie_s++){
+      int a_s = epsilon[ie_s][0]; //a'
+      int b_s = epsilon[ie_s][1]; //b'
+      int c_s = epsilon[ie_s][2]; //c'
+      for (int alpha_s=0; alpha_s<Ns; alpha_s++){
+      for (int beta_n=0; beta_n<Ns; beta_n++){
+        auto ee_GDGDDG = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd()(alpha_s,beta_n)(b_s,b_n);
+	auto DuGHDuG_ab_aa = DuGHDuG()(alpha_s, beta_n)(a_s,a_n);
+	auto DuG_ab_aa = DuG()(alpha_s,beta_n)(a_s,a_n);
+        for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+	  auto DuGHDuG_cb_ca = DuGHDuG()(gamma_s, beta_n)(c_s,a_n);
+          auto DuG_cb_ca = DuG()(gamma_s,beta_n)(c_s,a_n);
+          for (int gamma_n=0; gamma_n<Ns; gamma_n++){
+            result()(gamma_s,gamma_n)() += ee_GDGDDG * DuG_cb_ca * DuGHDu()(alpha_s,gamma_n)(a_s,c_n);
+            result()(gamma_s,gamma_n)() += ee_GDGDDG * DuGHDuG_cb_ca * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
+            result()(gamma_s,gamma_n)() -= ee_GDGDDG * DuGHDuG_ab_aa * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
+            result()(gamma_s,gamma_n)() -= ee_GDGDDG * DuG_ab_aa * DuGHDu()(gamma_s,gamma_n)(c_s,c_n);
+          }
+        }
+      }}
+    }
+  }
 }
 
 //Equivalent to "One-trace"
@@ -1137,7 +1230,7 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
  * Du_spec is a quark line from t_i to t_f
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
-template <class FImpl>
+/*template <class FImpl>
 template <class mobj, class mobj2, class robj>
 void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Du_spec,
@@ -1181,6 +1274,47 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
       }}}
     }
   }
+}*/
+template <class FImpl>
+template <class mobj, class mobj2, class robj>
+void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
+						 const mobj2 &Du_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_sigma,
+		                 		 const Gamma GammaB_nucl,
+						 robj &result)
+{
+
+  Gamma g5(Gamma::Algebra::Gamma5); 
+
+  auto DuG = Du_spec * GammaB_nucl;
+  // Gamma^B * Ds * \gamma_\mu^L * Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
+  auto GDsGDqGDd = GammaB_sigma * Ds_ti * Gamma_H * Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5;
+
+  for (int ie_n=0; ie_n < 6 ; ie_n++){
+    int a_n = epsilon[ie_n][0]; //a
+    int b_n = epsilon[ie_n][1]; //b
+    int c_n = epsilon[ie_n][2]; //c
+    for (int ie_s=0; ie_s < 6 ; ie_s++){
+      int a_s = epsilon[ie_s][0]; //a'
+      int b_s = epsilon[ie_s][1]; //b'
+      int c_s = epsilon[ie_s][2]; //c'
+      for (int alpha_s=0; alpha_s<Ns; alpha_s++){
+      for (int beta_n=0; beta_n<Ns; beta_n++){
+        auto ee_GDGDGD = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDqGDd()(alpha_s,beta_n)(b_s,b_n);
+	auto DuG_ab_aa = DuG()(alpha_s, beta_n)(a_s,a_n);
+        for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+          auto DuG_cb_ca = DuG()(gamma_s, beta_n)(c_s,a_n);
+          for (int gamma_n=0; gamma_n<Ns; gamma_n++){
+            result()(gamma_s,gamma_n)() -= ee_GDGDGD * DuG_ab_aa * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
+            result()(gamma_s,gamma_n)() += ee_GDGDGD * DuG_cb_ca * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
+          }
+	}
+      }}
+    }
+  }
 }
 
 /* Du_ti is a quark line from t_i to t_H
@@ -1188,7 +1322,7 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
  * Du_spec is a quark line from t_i to t_f
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
-template <class FImpl>
+/*template <class FImpl>
 template <class mobj, class mobj2, class robj>
 void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
 						 const mobj &Du_tf,
@@ -1240,6 +1374,57 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
       }}}
     }
   }
+}*/
+template <class FImpl>
+template <class mobj, class mobj2, class robj>
+void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
+						 const mobj &Du_tf,
+						 const mobj2 &Du_spec,
+						 const mobj &Dd_tf,
+						 const mobj &Ds_ti,
+				                 const Gamma Gamma_H,
+				                 const Gamma GammaB_sigma,
+		                 		 const Gamma GammaB_nucl,
+						 robj &result)
+{
+
+  Gamma g5(Gamma::Algebra::Gamma5); 
+
+  auto DuG = Du_spec * GammaB_nucl;
+  // Gamma^B * Ds * \gamma_\mu^L * adj(Du)
+  auto GDsGDu = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Du_tf) * g5;
+  // GDsGDu * GammaB
+  auto GDsGDuG = GDsGDu * GammaB_nucl;
+  // Du * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
+  auto DuGDd = Du_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
+
+  for (int ie_n=0; ie_n < 6 ; ie_n++){
+    int a_n = epsilon[ie_n][0]; //a
+    int b_n = epsilon[ie_n][1]; //b
+    int c_n = epsilon[ie_n][2]; //c
+    for (int ie_s=0; ie_s < 6 ; ie_s++){
+      int a_s = epsilon[ie_s][0]; //a'
+      int b_s = epsilon[ie_s][1]; //b'
+      int c_s = epsilon[ie_s][2]; //c'
+      auto ee = epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      for (int alpha_s=0; alpha_s<Ns; alpha_s++){
+      for (int beta_n=0; beta_n<Ns; beta_n++){
+        auto DuGDd_ab_ab = DuGDd()(alpha_s,beta_n)(a_s,b_n);
+	auto GDsGDuG_ab_ba = GDsGDuG()(alpha_s, beta_n)(b_s,a_n);
+        auto DuG_ab_aa = DuG()(alpha_s,beta_n)(a_s,a_n);
+        for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+          auto DuGDd_cb_cb = DuGDd()(gamma_s,beta_n)(c_s,b_n);
+	  auto DuG_cb_ca = DuG()(gamma_s,beta_n)(c_s,a_n);
+          for (int gamma_n=0; gamma_n<Ns; gamma_n++){
+            result()(gamma_s,gamma_n)() -= ee * DuGDd_ab_ab * DuG_cb_ca * GDsGDu()(alpha_s,gamma_n)(b_s,c_n);
+            result()(gamma_s,gamma_n)() -= ee * DuGDd_cb_cb * GDsGDuG_ab_ba * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
+            result()(gamma_s,gamma_n)() += ee * DuGDd_ab_ab * GDsGDuG_ab_ba * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
+            result()(gamma_s,gamma_n)() += ee * DuGDd_cb_cb * DuG_ab_aa * GDsGDu()(alpha_s,gamma_n)(b_s,c_n);
+          }
+        }
+      }}
+    }
+  }
 }
 
 template<class FImpl>

From 3594ce877beb22161dcce972f14bd89598972fd6 Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Tue, 3 Nov 2020 20:04:30 +0000
Subject: [PATCH 06/10] speedup in Sigma-to-nucleon

---
 Grid/qcd/utils/BaryonUtils.h | 196 -----------------------------------
 1 file changed, 196 deletions(-)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index c0999a4a..8a4ff6ac 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -1027,51 +1027,6 @@ void BaryonUtils<FImpl>::BaryonGamma3pt(
  * Du_spec is a quark line from t_i to t_f
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
-/*template <class FImpl>
-template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
-						 const mobj2 &Du_spec,
-						 const mobj &Dd_tf,
-						 const mobj &Ds_ti,
-				                 const Gamma Gamma_H,
-				                 const Gamma GammaB_sigma,
-		                 		 const Gamma GammaB_nucl,
-						 robj &result)
-{
-
-  Gamma g5(Gamma::Algebra::Gamma5); 
-
-  auto DuG = Du_spec * GammaB_nucl;
-  // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
-  auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
-  // Dq_loop * \gamma_\mu^L
-  auto DqG = Dq_loop * Gamma_H;
-
-  for (int ie_n=0; ie_n < 6 ; ie_n++){
-    int a_n = epsilon[ie_n][0]; //a
-    int b_n = epsilon[ie_n][1]; //b
-    int c_n = epsilon[ie_n][2]; //c
-    for (int ie_s=0; ie_s < 6 ; ie_s++){
-      int a_s = epsilon[ie_s][0]; //a'
-      int b_s = epsilon[ie_s][1]; //b'
-      int c_s = epsilon[ie_s][2]; //c'
-      for (int alpha_s=0; alpha_s<Ns; alpha_s++){
-      for (int beta_n=0; beta_n<Ns; beta_n++){
-        auto GDsGDd_ab_bb = GDsGDd()(alpha_s,beta_n)(b_s,b_n);
-        for (int tau2=0; tau2<Ns; tau2++){
-        for (int j=0; j<Nc; j++){
-          auto DqG_tt_jj = DqG()(tau2,tau2)(j,j);
-          auto ee_GDGDDG = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * DqG_tt_jj;
-          for (int gamma_s=0; gamma_s<Ns; gamma_s++){
-          for (int gamma_n=0; gamma_n<Ns; gamma_n++){
-            result()(gamma_s,gamma_n)() += ee_GDGDDG * DuG()(alpha_s, beta_n)(a_s,a_n) * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
-            result()(gamma_s,gamma_n)() -= ee_GDGDDG * DuG()(gamma_s, beta_n)(c_s,a_n) * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
-          }}
-	}}
-      }}
-    }
-  }
-}*/
 template <class FImpl>
 template <class mobj, class mobj2, class robj>
 void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
@@ -1122,59 +1077,6 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
  * Du_spec is a quark line from t_i to t_f
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
-/*template <class FImpl>
-template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
-						 const mobj &Du_tf,
-						 const mobj2 &Du_spec,
-						 const mobj &Dd_tf,
-						 const mobj &Ds_ti,
-				                 const Gamma Gamma_H,
-				                 const Gamma GammaB_sigma,
-		                 		 const Gamma GammaB_nucl,
-						 robj &result)
-{
-
-  Gamma g5(Gamma::Algebra::Gamma5); 
-
-  auto DuG = Du_spec * GammaB_nucl;
-  auto adjDu = g5 * adj(Du_tf) * g5;
-  auto adjDuG = adjDu * GammaB_nucl;
-  // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
-  auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
-  // Dq_loop * \gamma_\mu^L
-  auto DuGH = Du_ti * Gamma_H;
-
-  for (int ie_n=0; ie_n < 6 ; ie_n++){
-    int a_n = epsilon[ie_n][0]; //a
-    int b_n = epsilon[ie_n][1]; //b
-    int c_n = epsilon[ie_n][2]; //c
-    for (int ie_s=0; ie_s < 6 ; ie_s++){
-      int a_s = epsilon[ie_s][0]; //a'
-      int b_s = epsilon[ie_s][1]; //b'
-      int c_s = epsilon[ie_s][2]; //c'
-      for (int alpha_s=0; alpha_s<Ns; alpha_s++){
-      for (int beta_n=0; beta_n<Ns; beta_n++){
-        auto GDsGDd_ab_bb = GDsGDd()(alpha_s,beta_n)(b_s,b_n);
-        for (int tau2=0; tau2<Ns; tau2++){
-        for (int j=0; j<Nc; j++){
-          auto DuGH_at_aj = DuGH()(alpha_s,tau2)(a_s,j);
-          auto ee_GDGDDG_a = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * DuGH_at_aj;
-          for (int gamma_s=0; gamma_s<Ns; gamma_s++){
-            auto DuGH_gt_cj = DuGH()(gamma_s,tau2)(c_s,j);
-            auto ee_GDGDDG_c = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsGDd_ab_bb * DuGH_gt_cj;
-            for (int gamma_n=0; gamma_n<Ns; gamma_n++){
-              result()(gamma_s,gamma_n)() += ee_GDGDDG_a * DuG()(gamma_s, beta_n)(c_s,a_n) * adjDu()(tau2,gamma_n)(j,c_n);
-              result()(gamma_s,gamma_n)() += ee_GDGDDG_c * adjDuG()(tau2, beta_n)(j,a_n) * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
-              result()(gamma_s,gamma_n)() -= ee_GDGDDG_a * adjDuG()(tau2, beta_n)(j,a_n) * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
-              result()(gamma_s,gamma_n)() -= ee_GDGDDG_c * DuG()(alpha_s, beta_n)(a_s,a_n) * adjDu()(tau2,gamma_n)(j,c_n);
-            }
-	  }
-	}}
-      }}
-    }
-  }
-}*/
 template <class FImpl>
 template <class mobj, class mobj2, class robj>
 void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
@@ -1230,51 +1132,6 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
  * Du_spec is a quark line from t_i to t_f
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
-/*template <class FImpl>
-template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
-						 const mobj2 &Du_spec,
-						 const mobj &Dd_tf,
-						 const mobj &Ds_ti,
-				                 const Gamma Gamma_H,
-				                 const Gamma GammaB_sigma,
-		                 		 const Gamma GammaB_nucl,
-						 robj &result)
-{
-
-  Gamma g5(Gamma::Algebra::Gamma5); 
-
-  auto DuG = Du_spec * GammaB_nucl;
-  // Gamma^B * Ds * \gamma_\mu^L
-  auto GDsG = GammaB_sigma * Ds_ti * Gamma_H;
-  // Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
-  auto DqGDd = Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5;
-
-  for (int ie_n=0; ie_n < 6 ; ie_n++){
-    int a_n = epsilon[ie_n][0]; //a
-    int b_n = epsilon[ie_n][1]; //b
-    int c_n = epsilon[ie_n][2]; //c
-    for (int ie_s=0; ie_s < 6 ; ie_s++){
-      int a_s = epsilon[ie_s][0]; //a'
-      int b_s = epsilon[ie_s][1]; //b'
-      int c_s = epsilon[ie_s][2]; //c'
-      for (int alpha_s=0; alpha_s<Ns; alpha_s++){
-      for (int tau=0; tau<Ns; tau++){
-      for (int i=0; i<Nc; i++){
-	auto GDsG_at_bi = GDsG()(alpha_s,tau)(b_s,i);
-        for (int beta_n=0; beta_n<Ns; beta_n++){
-          auto DqGDd_tb_ib = DqGDd()(tau,beta_n)(i,b_n);
-	  auto ee_GDGDGD = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DqGDd_tb_ib;
-          for (int gamma_s=0; gamma_s<Ns; gamma_s++){
-          for (int gamma_n=0; gamma_n<Ns; gamma_n++){
-            result()(gamma_s,gamma_n)() -= ee_GDGDGD * DuG()(alpha_s, beta_n)(a_s,a_n) * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
-            result()(gamma_s,gamma_n)() += ee_GDGDGD * DuG()(gamma_s, beta_n)(c_s,a_n) * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
-          }}
-	}
-      }}}
-    }
-  }
-}*/
 template <class FImpl>
 template <class mobj, class mobj2, class robj>
 void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
@@ -1322,59 +1179,6 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
  * Du_spec is a quark line from t_i to t_f
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
-/*template <class FImpl>
-template <class mobj, class mobj2, class robj>
-void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
-						 const mobj &Du_tf,
-						 const mobj2 &Du_spec,
-						 const mobj &Dd_tf,
-						 const mobj &Ds_ti,
-				                 const Gamma Gamma_H,
-				                 const Gamma GammaB_sigma,
-		                 		 const Gamma GammaB_nucl,
-						 robj &result)
-{
-
-  Gamma g5(Gamma::Algebra::Gamma5); 
-
-  auto DuG = Du_spec * GammaB_nucl;
-  auto adjDu = g5 * adj(Du_tf) * g5;
-  auto adjDuG = adjDu * GammaB_nucl;
-  // Gamma^B * Ds * \gamma_\mu^L
-  auto GDsG = GammaB_sigma * Ds_ti * Gamma_H;
-  // Du * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5)
-  auto DuGDd = Du_ti * Gamma_H * g5 * adj(Dd_tf) * g5;
-
-  for (int ie_n=0; ie_n < 6 ; ie_n++){
-    int a_n = epsilon[ie_n][0]; //a
-    int b_n = epsilon[ie_n][1]; //b
-    int c_n = epsilon[ie_n][2]; //c
-    for (int ie_s=0; ie_s < 6 ; ie_s++){
-      int a_s = epsilon[ie_s][0]; //a'
-      int b_s = epsilon[ie_s][1]; //b'
-      int c_s = epsilon[ie_s][2]; //c'
-      for (int alpha_s=0; alpha_s<Ns; alpha_s++){
-      for (int tau=0; tau<Ns; tau++){
-      for (int i=0; i<Nc; i++){
-	auto GDsG_at_bi = GDsG()(alpha_s,tau)(b_s,i);
-        for (int beta_n=0; beta_n<Ns; beta_n++){
-          auto DuGDd_ab_ab = DuGDd()(alpha_s,beta_n)(a_s,b_n);
-	  auto ee_GDGDGD_a = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DuGDd_ab_ab;
-          for (int gamma_s=0; gamma_s<Ns; gamma_s++){
-            auto DuGDd_gb_cb = DuGDd()(gamma_s,beta_n)(c_s,b_n);
-	    auto ee_GDGDGD_c = epsilon_sgn[ie_n] * epsilon_sgn[ie_s] * GDsG_at_bi * DuGDd_gb_cb;
-            for (int gamma_n=0; gamma_n<Ns; gamma_n++){
-              result()(gamma_s,gamma_n)() -= ee_GDGDGD_a * DuG()(gamma_s, beta_n)(c_s,a_n) * adjDu()(tau,gamma_n)(i,c_n);
-              result()(gamma_s,gamma_n)() -= ee_GDGDGD_c * adjDuG()(tau, beta_n)(i,a_n) * Du_spec()(alpha_s,gamma_n)(a_s,c_n);
-              result()(gamma_s,gamma_n)() += ee_GDGDGD_a * adjDuG()(tau, beta_n)(i,a_n) * Du_spec()(gamma_s,gamma_n)(c_s,c_n);
-              result()(gamma_s,gamma_n)() += ee_GDGDGD_c * DuG()(alpha_s, beta_n)(a_s,a_n) * adjDu()(tau,gamma_n)(i,c_n);
-            }
-	  }
-	}
-      }}}
-    }
-  }
-}*/
 template <class FImpl>
 template <class mobj, class mobj2, class robj>
 void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,

From 9b9a53f87066b5c67b607f56cb73c56efa5d4243 Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Tue, 2 Feb 2021 13:06:43 +0000
Subject: [PATCH 07/10] ...

---
 Grid/qcd/utils/BaryonUtils.h | 53 ++++++++++++++++++++++--------------
 1 file changed, 32 insertions(+), 21 deletions(-)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index 8a4ff6ac..94cc07b1 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -1350,15 +1350,19 @@ void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
   // \gamma_\mu^L * Dq_loop 
   auto trGDq = TensorRemove(trace(Gamma_H * Dq_loop)); 
 
+  Real ee;
+
   for (int ie_s=0; ie_s < 6 ; ie_s++){
-    int a_s = epsilon[ie_s][0]; //a
-    int b_s = epsilon[ie_s][1]; //b
-    int c_s = epsilon[ie_s][2]; //c
+      int a_s = (ie_s < 3 ? ie_s       : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a'
+      int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b'
+      int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_s = (ie_s < 3 ? 1 : -1);
     for (int ie_x=0; ie_x < 6 ; ie_x++){
-      int a_x = epsilon[ie_x][0]; //a'
-      int b_x = epsilon[ie_x][1]; //b'
-      int c_x = epsilon[ie_x][2]; //c'
-      auto ee_GD = epsilon_sgn[ie_s] * epsilon_sgn[ie_x] * trGDq; 
+      int a_x = (ie_x < 3 ? ie_x       : (6-ie_x)%3 ); //epsilon[ie_x][0]; //a'
+      int b_x = (ie_x < 3 ? (ie_x+1)%3 : (8-ie_x)%3 ); //epsilon[ie_x][1]; //b'
+      int c_x = (ie_x < 3 ? (ie_x+2)%3 : (7-ie_x)%3 ); //epsilon[ie_x][2]; //c'
+      int eSgn_x = (ie_x < 3 ? 1 : -1);
+      ee = Real(eSgn_s * eSgn_x);
       for (int alpha_x=0; alpha_x<Ns; alpha_x++){
       for (int beta_s=0; beta_s<Ns; beta_s++){
         auto GDsGDdG_ab_ba = GDsGDd()(alpha_x,beta_s)(b_x,a_s);
@@ -1451,7 +1455,7 @@ void BaryonUtils<FImpl>::XiToSigmaEye(const PropagatorField &qq_loop,
 				                 const Gamma GammaB_xi,
 		                 		 const Gamma GammaB_sigma,
 						 const std::string op,
-						 SpinMatrixField &stn_corr)
+						 SpinMatrixField &xts_corr)
 {
 
   assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
@@ -1459,24 +1463,31 @@ void BaryonUtils<FImpl>::XiToSigmaEye(const PropagatorField &qq_loop,
 
   GridBase *grid = qs_ti.Grid();
 
-  autoView( vcorr, stn_corr, CpuWrite);
-  autoView( vq_loop , qq_loop, CpuRead);
-  autoView( vd_tf , qd_tf, CpuRead);
-  autoView( vs_ti , qs_ti, CpuRead);
+  autoView( vcorr   , xts_corr , AcceleratorWrite);
+  autoView( vq_loop , qq_loop  , AcceleratorRead);
+  autoView( vd_tf   , qd_tf    , AcceleratorRead);
+  autoView( vs_ti   , qs_ti    , AcceleratorRead);
+
+  bool doQ1 = (op == "Q1");
+  bool doQ2 = (op == "Q2");
+
+  Vector<mobj> my_Dq_spec{Dd_spec,Ds_spec};
+  mobj * Dq_spec_p = &my_Dq_spec[0];
 
   accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
-    auto Dq_loop = vq_loop[ss];
-    auto Dd_tf = vd_tf[ss];
-    auto Ds_ti = vs_ti[ss];
-    sobj result=Zero();
-    if(op == "Q1"){
-      XiToSigmaQ1EyeSite(Dq_loop,Dd_spec,Ds_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
-    } else if(op == "Q2"){
-      XiToSigmaQ2EyeSite(Dq_loop,Dd_spec,Ds_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
+    auto Dq_loop = vq_loop(ss);
+    auto Dd_tf   = vd_tf(ss);
+    auto Ds_ti   = vs_ti(ss);
+    typedef decltype(coalescedRead(vcorr[0])) spinor;
+    spinor result=Zero();
+    if(doQ1){
+      XiToSigmaQ1EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[1],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
+    } else if(doQ2){
+      XiToSigmaQ2EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
     } else {
       assert(0 && "Weak Operator not correctly specified");
     }
-      vcorr[ss] = result; 
+    coalescedWrite(vcorr[ss],result);
   }  );//end loop over lattice sites
 }
 

From d4b4de8f428c625757c31c6d265c288671a0ef5b Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Thu, 4 Mar 2021 20:01:24 +0000
Subject: [PATCH 08/10] changes

---
 Grid/qcd/utils/BaryonUtils.h | 60 +++++++++++++++++++++---------------
 1 file changed, 35 insertions(+), 25 deletions(-)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index b74a5b20..1a7a4d38 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -252,7 +252,7 @@ public:
 				 const Gamma GammaB_sigma,
 				 const Gamma GammaB_nucl,
 		                 const std::string op,
-				 SpinMatrixField &stn_corr);
+				 SpinMatrixField &xts_corr);
 };
 //This computes a baryon contraction on a lattice site, including the spin-trace of the correlation matrix
 template <class FImpl>
@@ -1435,6 +1435,7 @@ void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
       int c_x = (ie_x < 3 ? (ie_x+2)%3 : (7-ie_x)%3 ); //epsilon[ie_x][2]; //c'
       int eSgn_x = (ie_x < 3 ? 1 : -1);
       ee = Real(eSgn_s * eSgn_x);
+      auto ee_GD = ee * trGDq;
       for (int alpha_x=0; alpha_x<Ns; alpha_x++){
       for (int beta_s=0; beta_s<Ns; beta_s++){
         auto GDsGDdG_ab_ba = GDsGDd()(alpha_x,beta_s)(b_x,a_s);
@@ -1488,15 +1489,19 @@ void BaryonUtils<FImpl>::XiToSigmaQ2EyeSite(const mobj &Dq_loop,
   // GammaB * DsGDd * GammaB
   auto GDsGDqGDdG = GDsGDqGDd * GammaB_sigma;
 
+  Real ee;
+
   for (int ie_s=0; ie_s < 6 ; ie_s++){
-    int a_s = epsilon[ie_s][0]; //a
-    int b_s = epsilon[ie_s][1]; //b
-    int c_s = epsilon[ie_s][2]; //c
+      int a_s = (ie_s < 3 ? ie_s       : (6-ie_s)%3 ); //epsilon[ie_s][0]; //a'
+      int b_s = (ie_s < 3 ? (ie_s+1)%3 : (8-ie_s)%3 ); //epsilon[ie_s][1]; //b'
+      int c_s = (ie_s < 3 ? (ie_s+2)%3 : (7-ie_s)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_s = (ie_s < 3 ? 1 : -1);
     for (int ie_x=0; ie_x < 6 ; ie_x++){
-      int a_x = epsilon[ie_x][0]; //a'
-      int b_x = epsilon[ie_x][1]; //b'
-      int c_x = epsilon[ie_x][2]; //c'
-      auto ee = epsilon_sgn[ie_s] * epsilon_sgn[ie_x];
+      int a_x = (ie_x < 3 ? ie_x       : (6-ie_x)%3 ); //epsilon[ie_x][0]; //a'
+      int b_x = (ie_x < 3 ? (ie_x+1)%3 : (8-ie_x)%3 ); //epsilon[ie_x][1]; //b'
+      int c_x = (ie_x < 3 ? (ie_x+2)%3 : (7-ie_x)%3 ); //epsilon[ie_x][2]; //c'
+      int eSgn_x = (ie_x < 3 ? 1 : -1);
+      ee = Real(eSgn_s * eSgn_x);
       for (int alpha_x=0; alpha_x<Ns; alpha_x++){
       for (int beta_s=0; beta_s<Ns; beta_s++){
         auto GDsGDqGDdG_ab_ba = GDsGDqGDdG()(alpha_x,beta_s)(b_x,a_s);
@@ -1540,27 +1545,32 @@ void BaryonUtils<FImpl>::XiToSigmaEye(const PropagatorField &qq_loop,
   autoView( vd_tf   , qd_tf    , AcceleratorRead);
   autoView( vs_ti   , qs_ti    , AcceleratorRead);
 
-  bool doQ1 = (op == "Q1");
-  bool doQ2 = (op == "Q2");
-
   Vector<mobj> my_Dq_spec{Dd_spec,Ds_spec};
   mobj * Dq_spec_p = &my_Dq_spec[0];
 
-  accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
-    auto Dq_loop = vq_loop(ss);
-    auto Dd_tf   = vd_tf(ss);
-    auto Ds_ti   = vs_ti(ss);
-    typedef decltype(coalescedRead(vcorr[0])) spinor;
-    spinor result=Zero();
-    if(doQ1){
+  if(op == "Q1"){
+    accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+      auto Dq_loop = vq_loop(ss);
+      auto Dd_tf   = vd_tf(ss);
+      auto Ds_ti   = vs_ti(ss);
+      typedef decltype(coalescedRead(vcorr[0])) spinor;
+      spinor result=Zero();
       XiToSigmaQ1EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[1],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
-    } else if(doQ2){
-      XiToSigmaQ2EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[0],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
-    } else {
-      assert(0 && "Weak Operator not correctly specified");
-    }
-    coalescedWrite(vcorr[ss],result);
-  }  );//end loop over lattice sites
+      coalescedWrite(vcorr[ss],result);
+    }  );//end loop over lattice sites
+  } else if(op == "Q2"){
+    accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+      auto Dq_loop = vq_loop(ss);
+      auto Dd_tf   = vd_tf(ss);
+      auto Ds_ti   = vs_ti(ss);
+      typedef decltype(coalescedRead(vcorr[0])) spinor;
+      spinor result=Zero();
+      XiToSigmaQ2EyeSite(Dq_loop,Dq_spec_p[0],Dq_spec_p[1],Dd_tf,Ds_ti,Gamma_H,GammaB_xi,GammaB_sigma,result);
+      coalescedWrite(vcorr[ss],result);
+    }  );//end loop over lattice sites
+  } else {
+    assert(0 && "Weak Operator not correctly specified");
+  }
 }
 
 

From 9b15704290048b936e71c3ddd652a6eec9bf5fff Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Fri, 5 Mar 2021 10:42:32 +0000
Subject: [PATCH 09/10] tested and consitent

---
 Grid/qcd/utils/BaryonUtils.h | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index 1a7a4d38..b69865e8 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -27,7 +27,6 @@
  *************************************************************************************/
 /*  END LEGAL */
 #pragma once
-//#include <Grid/Hadrons/Global.hpp>
 #include <Grid/Eigen/unsupported/CXX11/Tensor>
 
 NAMESPACE_BEGIN(Grid);
@@ -200,7 +199,7 @@ public:
 				                 const Gamma GammaB_sigma,
 		                 		 const Gamma GammaB_nucl,
 						 robj &result);
-  template <class mobj, class mobj2, class robj>
+  template <class mobj, class mobj2, class robj> accelerator_inline
   static void XiToSigmaQ1EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Dd_spec,
 						 const mobj2 &Ds_spec,
@@ -210,7 +209,7 @@ public:
 				                 const Gamma GammaB_sigma,
 		                 		 const Gamma GammaB_nucl,
 						 robj &result);
-  template <class mobj, class mobj2, class robj>
+  template <class mobj, class mobj2, class robj> accelerator_inline
   static void XiToSigmaQ2EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Dd_spec,
 						 const mobj2 &Ds_spec,
@@ -1395,7 +1394,7 @@ void BaryonUtils<FImpl>::SigmaToNucleonNonEye(const PropagatorField &qq_ti,
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
 template <class FImpl>
-template <class mobj, class mobj2, class robj>
+template <class mobj, class mobj2, class robj> accelerator_inline
 void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Dd_spec,
 						 const mobj2 &Ds_spec,
@@ -1464,7 +1463,7 @@ void BaryonUtils<FImpl>::XiToSigmaQ1EyeSite(const mobj &Dq_loop,
  * Dd_tf is a quark line from t_f to t_H
  * Ds_ti is a quark line from t_i to t_H */
 template <class FImpl>
-template <class mobj, class mobj2, class robj>
+template <class mobj, class mobj2, class robj> accelerator_inline
 void BaryonUtils<FImpl>::XiToSigmaQ2EyeSite(const mobj &Dq_loop,
 						 const mobj2 &Dd_spec,
 						 const mobj2 &Ds_spec,

From 7a19432e0b77fff8943ad6701f448ab366dab27f Mon Sep 17 00:00:00 2001
From: Felix Erben <felix.erben@ed.ac.uk>
Date: Fri, 5 Mar 2021 10:57:09 +0000
Subject: [PATCH 10/10] whitespace

---
 Grid/qcd/utils/BaryonUtils.h | 84 ++++++++++++++++++------------------
 1 file changed, 42 insertions(+), 42 deletions(-)

diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index b69865e8..9d9cb508 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -191,34 +191,34 @@ public:
              robj &result);
   template <class mobj, class mobj2, class robj> accelerator_inline
   static void SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
-						 const mobj &Du_tf,
-						 const mobj2 &Du_spec,
-						 const mobj &Dd_tf,
-						 const mobj &Ds_ti,
-				                 const Gamma Gamma_H,
-				                 const Gamma GammaB_sigma,
-		                 		 const Gamma GammaB_nucl,
-						 robj &result);
+	     const mobj &Du_tf,
+	     const mobj2 &Du_spec,
+	     const mobj &Dd_tf,
+	     const mobj &Ds_ti,
+	     const Gamma Gamma_H,
+	     const Gamma GammaB_sigma,
+	     const Gamma GammaB_nucl,
+	     robj &result);
   template <class mobj, class mobj2, class robj> accelerator_inline
   static void XiToSigmaQ1EyeSite(const mobj &Dq_loop,
-						 const mobj2 &Dd_spec,
-						 const mobj2 &Ds_spec,
-						 const mobj &Dd_tf,
-						 const mobj &Ds_ti,
-				                 const Gamma Gamma_H,
-				                 const Gamma GammaB_sigma,
-		                 		 const Gamma GammaB_nucl,
-						 robj &result);
+	     const mobj2 &Dd_spec,
+	     const mobj2 &Ds_spec,
+	     const mobj &Dd_tf,
+	     const mobj &Ds_ti,
+	     const Gamma Gamma_H,
+	     const Gamma GammaB_sigma,
+	     const Gamma GammaB_nucl,
+	     robj &result);
   template <class mobj, class mobj2, class robj> accelerator_inline
   static void XiToSigmaQ2EyeSite(const mobj &Dq_loop,
-						 const mobj2 &Dd_spec,
-						 const mobj2 &Ds_spec,
-						 const mobj &Dd_tf,
-						 const mobj &Ds_ti,
-				                 const Gamma Gamma_H,
-				                 const Gamma GammaB_sigma,
-		                 		 const Gamma GammaB_nucl,
-						 robj &result);
+	     const mobj2 &Dd_spec,
+	     const mobj2 &Ds_spec,
+	     const mobj &Dd_tf,
+	     const mobj &Ds_ti,
+	     const Gamma Gamma_H,
+	     const Gamma GammaB_sigma,
+	     const Gamma GammaB_nucl,
+	     robj &result);
   public:
   template <class mobj>
   static void SigmaToNucleonEye(const PropagatorField &qq_loop,
@@ -232,26 +232,26 @@ public:
          SpinMatrixField &stn_corr);
   template <class mobj>
   static void SigmaToNucleonNonEye(const PropagatorField &qq_ti,
-				 const PropagatorField &qq_tf,
-				 const mobj &Du_spec,
-				 const PropagatorField &qd_tf,
-				 const PropagatorField &qs_ti,
-				 const Gamma Gamma_H,
-				 const Gamma GammaB_sigma,
-				 const Gamma GammaB_nucl,
-		                 const std::string op,
-				 SpinMatrixField &stn_corr);
+	 const PropagatorField &qq_tf,
+	 const mobj &Du_spec,
+	 const PropagatorField &qd_tf,
+	 const PropagatorField &qs_ti,
+	 const Gamma Gamma_H,
+	 const Gamma GammaB_sigma,
+	 const Gamma GammaB_nucl,
+         const std::string op,
+	 SpinMatrixField &stn_corr);
   template <class mobj>
   static void XiToSigmaEye(const PropagatorField &qq_loop,
-				 const mobj &Dd_spec,
-				 const mobj &Ds_spec,
-				 const PropagatorField &qd_tf,
-				 const PropagatorField &qs_ti,
-				 const Gamma Gamma_H,
-				 const Gamma GammaB_sigma,
-				 const Gamma GammaB_nucl,
-		                 const std::string op,
-				 SpinMatrixField &xts_corr);
+	 const mobj &Dd_spec,
+	 const mobj &Ds_spec,
+	 const PropagatorField &qd_tf,
+	 const PropagatorField &qs_ti,
+	 const Gamma Gamma_H,
+	 const Gamma GammaB_sigma,
+	 const Gamma GammaB_nucl,
+         const std::string op,
+	 SpinMatrixField &xts_corr);
 };
 //This computes a baryon contraction on a lattice site, including the spin-trace of the correlation matrix
 template <class FImpl>