diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h
index 25c71e3a..35358d05 100644
--- a/Grid/qcd/utils/BaryonUtils.h
+++ b/Grid/qcd/utils/BaryonUtils.h
@@ -40,14 +40,7 @@ public:
   typedef typename FImpl::FermionField FermionField;
   typedef typename FImpl::PropagatorField PropagatorField;
 
-  typedef typename FImpl::SitePropagator pobj;
-  typedef typename ComplexField::vector_object vobj;
-
   typedef Lattice<iSpinMatrix<typename FImpl::Simd>> SpinMatrixField;
-  //typedef typename SpinMatrixField::vector_object sobj;
-
-  //static const int epsilon[6][3] ;
-  //static const Real epsilon_sgn[6];
 
   private: 
   template <class mobj, class robj> accelerator_inline
@@ -122,7 +115,7 @@ public:
   static void BaryonGamma3ptGroup1Site(
            const mobj &Dq1_ti,
            const mobj2 &Dq2_spec,
-          // const mobj2 &Dq3_spec,
+           const mobj2 &Dq3_spec,
            const mobj &Dq4_tf,
                    const Gamma GammaJ,
                    const Gamma GammaBi,
@@ -134,7 +127,7 @@ public:
   static void BaryonGamma3ptGroup2Site(
            const mobj2 &Dq1_spec,
            const mobj &Dq2_ti,
-           //const mobj2 &Dq3_spec,
+           const mobj2 &Dq3_spec,
            const mobj &Dq4_tf,
                    const Gamma GammaJ,
                    const Gamma GammaBi,
@@ -145,7 +138,7 @@ public:
   template <class mobj, class mobj2, class robj> accelerator_inline
   static void BaryonGamma3ptGroup3Site(
            const mobj2 &Dq1_spec,
-           //const mobj2 &Dq2_spec,
+           const mobj2 &Dq2_spec,
            const mobj &Dq3_ti,
            const mobj &Dq4_tf,
                    const Gamma GammaJ,
@@ -230,13 +223,7 @@ public:
                      const std::string op,
          SpinMatrixField &stn_corr);
 };
-/*
-template <class FImpl> 
-const int BaryonUtils<FImpl>::epsilon[6][3] = {{0,1,2},{1,2,0},{2,0,1},{0,2,1},{2,1,0},{1,0,2}};
-template <class FImpl> 
-const Real BaryonUtils<FImpl>::epsilon_sgn[6] = {1.,1.,1.,-1.,-1.,-1.};
-*/
-//This is the old version
+//This computes a baryon contraction on a lattice site, including the spin-trace of the correlation matrix
 template <class FImpl>
 template <class mobj, class robj> accelerator_inline
 void BaryonUtils<FImpl>::BaryonSite(const mobj &D1,
@@ -251,20 +238,20 @@ void BaryonUtils<FImpl>::BaryonSite(const mobj &D1,
                 robj &result)
 {
 
-    Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4)
+  Gamma g4(Gamma::Algebra::GammaT); //needed for parity P_\pm = 0.5*(1 \pm \gamma_4)
     
-    auto D1_GAi =  D1 * GammaA_i;
-    auto D1_GAi_g4 = D1_GAi * g4;
-    auto D1_GAi_P = 0.5*(D1_GAi + (Real)parity * D1_GAi_g4);
-    auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P;
-    auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P;
+  auto D1_GAi =  D1 * GammaA_i;
+  auto D1_GAi_g4 = D1_GAi * g4;
+  auto D1_GAi_P = 0.5*(D1_GAi + (Real)parity * D1_GAi_g4);
+  auto GAf_D1_GAi_P = GammaA_f * D1_GAi_P;
+  auto GBf_D1_GAi_P = GammaB_f * D1_GAi_P;
 
-    auto D2_GBi = D2 * GammaB_i;
-    auto GBf_D2_GBi = GammaB_f * D2_GBi;
-    auto GAf_D2_GBi = GammaA_f * D2_GBi;
+  auto D2_GBi = D2 * GammaB_i;
+  auto GBf_D2_GBi = GammaB_f * D2_GBi;
+  auto GAf_D2_GBi = GammaA_f * D2_GBi;
 
-    auto GBf_D3 = GammaB_f * D3;
-    auto GAf_D3 = GammaA_f * D3;
+  auto GBf_D3 = GammaB_f * D3;
+  auto GAf_D3 = GammaA_f * D3;
 
   Real ee;
 
@@ -273,86 +260,87 @@ void BaryonUtils<FImpl>::BaryonSite(const mobj &D1,
     int b_f    = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b
     int c_f    = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_f = (ie_f < 3 ? 1 : -1);
-  for (int ie_i=0; ie_i < 6 ; ie_i++){
-    int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
-    int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
-    int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
-    int eSgn_i = (ie_i < 3 ? 1 : -1);
+    for (int ie_i=0; ie_i < 6 ; ie_i++){
+      int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
+      int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
+      int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_i = (ie_i < 3 ? 1 : -1);
 
-    ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
-        //This is the \delta_{456}^{123} part
-        if (wick_contraction[0]){
-            for (int rho=0; rho<Ns; rho++){
-                auto GAf_D1_GAi_P_rr_cc = GAf_D1_GAi_P()(rho,rho)(c_f,c_i);
-                for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                    result()()() += ee  * GAf_D1_GAi_P_rr_cc
+      ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      //This is the \delta_{456}^{123} part
+      if (wick_contraction[0]){
+        for (int rho=0; rho<Ns; rho++){
+          auto GAf_D1_GAi_P_rr_cc = GAf_D1_GAi_P()(rho,rho)(c_f,c_i);
+          for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            result()()() += ee  * GAf_D1_GAi_P_rr_cc
                                         * D2_GBi    ()(alpha_f,beta_i)(a_f,a_i)
                                         * GBf_D3    ()(alpha_f,beta_i)(b_f,b_i);
-                }}
-            }
-        }   
-        //This is the \delta_{456}^{231} part
-        if (wick_contraction[1]){
-            for (int rho=0; rho<Ns; rho++){
-            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                auto D1_GAi_P_ar_ac = D1_GAi_P()(alpha_f,rho)(a_f,c_i);
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                    result()()() += ee  * D1_GAi_P_ar_ac
+          }}
+        }
+      }   
+      //This is the \delta_{456}^{231} part
+      if (wick_contraction[1]){
+        for (int rho=0; rho<Ns; rho++){
+        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          auto D1_GAi_P_ar_ac = D1_GAi_P()(alpha_f,rho)(a_f,c_i);
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            result()()() += ee  * D1_GAi_P_ar_ac
                                         * GBf_D2_GBi    ()(alpha_f,beta_i)(b_f,a_i)
                                         * GAf_D3        ()(rho,beta_i)(c_f,b_i);
-                }
-            }}
-        }   
-        //This is the \delta_{456}^{312} part
-        if (wick_contraction[2]){
-            for (int rho=0; rho<Ns; rho++){
-                for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                auto GBf_D1_GAi_P_ar_bc = GBf_D1_GAi_P()(alpha_f,rho)(b_f,c_i);
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                    result()()() += ee  * GBf_D1_GAi_P_ar_bc
+          }
+        }}
+      }   
+      //This is the \delta_{456}^{312} part
+      if (wick_contraction[2]){
+        for (int rho=0; rho<Ns; rho++){
+        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          auto GBf_D1_GAi_P_ar_bc = GBf_D1_GAi_P()(alpha_f,rho)(b_f,c_i);
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            result()()() += ee  * GBf_D1_GAi_P_ar_bc
                                         * GAf_D2_GBi    ()(rho,beta_i)(c_f,a_i)
                                         * D3            ()(alpha_f,beta_i)(a_f,b_i);
-                }
-            }}
-        }   
-        //This is the \delta_{456}^{132} part
-        if (wick_contraction[3]){
-            for (int rho=0; rho<Ns; rho++){
-                auto GAf_D1_GAi_P_rr_cc = GAf_D1_GAi_P()(rho,rho)(c_f,c_i);
-                for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                    result()()() -= ee  * GAf_D1_GAi_P_rr_cc
+          }
+        }}
+      }   
+      //This is the \delta_{456}^{132} part
+      if (wick_contraction[3]){
+        for (int rho=0; rho<Ns; rho++){
+          auto GAf_D1_GAi_P_rr_cc = GAf_D1_GAi_P()(rho,rho)(c_f,c_i);
+          for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            result()()() -= ee  * GAf_D1_GAi_P_rr_cc
                                         * GBf_D2_GBi    ()(alpha_f,beta_i)(b_f,a_i)
                                         * D3            ()(alpha_f,beta_i)(a_f,b_i);
-              }
-        }}
-        }   
-        //This is the \delta_{456}^{321} part
-        if (wick_contraction[4]){
-            for (int rho=0; rho<Ns; rho++){
-            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                auto GBf_D1_GAi_P_ar_bc = GBf_D1_GAi_P()(alpha_f,rho)(b_f,c_i);
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                    result()()() -= ee  * GBf_D1_GAi_P_ar_bc
+          }}
+        }
+      }   
+      //This is the \delta_{456}^{321} part
+      if (wick_contraction[4]){
+        for (int rho=0; rho<Ns; rho++){
+        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          auto GBf_D1_GAi_P_ar_bc = GBf_D1_GAi_P()(alpha_f,rho)(b_f,c_i);
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            result()()() -= ee  * GBf_D1_GAi_P_ar_bc
                                         * D2_GBi    ()(alpha_f,beta_i)(a_f,a_i)
                                         * GAf_D3    ()(rho,beta_i)(c_f,b_i);
-                }
-            }}
-        }   
-        //This is the \delta_{456}^{213} part
-        if (wick_contraction[5]){
-            for (int rho=0; rho<Ns; rho++){
-            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                auto D1_GAi_P_ar_ac = D1_GAi_P()(alpha_f,rho)(a_f,c_i);
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                    result()()() -= ee  * D1_GAi_P_ar_ac
+          }
+        }}
+      }   
+      //This is the \delta_{456}^{213} part
+      if (wick_contraction[5]){
+        for (int rho=0; rho<Ns; rho++){
+        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          auto D1_GAi_P_ar_ac = D1_GAi_P()(alpha_f,rho)(a_f,c_i);
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            result()()() -= ee  * D1_GAi_P_ar_ac
                                         * GAf_D2_GBi    ()(rho,beta_i)(c_f,a_i)
                                         * GBf_D3        ()(alpha_f,beta_i)(b_f,b_i);
-                }
-            }}
-        }
-    }}
+          }
+        }}
+      }
+    }
+  }
 }
 
 //New version without parity projection or trace
@@ -369,17 +357,16 @@ void BaryonUtils<FImpl>::BaryonSiteMatrix(const mobj &D1,
                 robj &result)
 {
 
-    auto D1_GAi =  D1 * GammaA_i;
-    auto GAf_D1_GAi = GammaA_f * D1_GAi;
-    auto GBf_D1_GAi = GammaB_f * D1_GAi;
+  auto D1_GAi =  D1 * GammaA_i;
+  auto GAf_D1_GAi = GammaA_f * D1_GAi;
+  auto GBf_D1_GAi = GammaB_f * D1_GAi;
 
-    auto D2_GBi = D2 * GammaB_i;
-    auto GBf_D2_GBi = GammaB_f * D2_GBi;
-    auto GAf_D2_GBi = GammaA_f * D2_GBi;
-
-    auto GBf_D3 = GammaB_f * D3;
-    auto GAf_D3 = GammaA_f * D3;
+  auto D2_GBi = D2 * GammaB_i;
+  auto GBf_D2_GBi = GammaB_f * D2_GBi;
+  auto GAf_D2_GBi = GammaA_f * D2_GBi;
 
+  auto GBf_D3 = GammaB_f * D3;
+  auto GAf_D3 = GammaA_f * D3;
 
   Real ee;
 
@@ -388,96 +375,101 @@ void BaryonUtils<FImpl>::BaryonSiteMatrix(const mobj &D1,
     int b_f    = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b
     int c_f    = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_f = (ie_f < 3 ? 1 : -1);
-  for (int ie_i=0; ie_i < 6 ; ie_i++){
-    int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
-    int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
-    int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
-    int eSgn_i = (ie_i < 3 ? 1 : -1);
+    for (int ie_i=0; ie_i < 6 ; ie_i++){
+      int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
+      int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
+      int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_i = (ie_i < 3 ? 1 : -1);
 
-    ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
-        //This is the \delta_{456}^{123} part
-        if (wick_contraction[0]){
-            for (int rho_i=0; rho_i<Ns; rho_i++){
-            for (int rho_f=0; rho_f<Ns; rho_f++){
-                auto GAf_D1_GAi_rr_cc = GAf_D1_GAi()(rho_f,rho_i)(c_f,c_i);
-                for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                    result()(rho_f,rho_i)() += ee  * GAf_D1_GAi_rr_cc
+      ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      //This is the \delta_{456}^{123} part
+      if (wick_contraction[0]){
+        for (int rho_i=0; rho_i<Ns; rho_i++){
+        for (int rho_f=0; rho_f<Ns; rho_f++){
+          auto GAf_D1_GAi_rr_cc = GAf_D1_GAi()(rho_f,rho_i)(c_f,c_i);
+          for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            result()(rho_f,rho_i)() += ee  * GAf_D1_GAi_rr_cc
                                         * D2_GBi    ()(alpha_f,beta_i)(a_f,a_i)
                                         * GBf_D3    ()(alpha_f,beta_i)(b_f,b_i);
-                }}
-            }}
-        }   
-        //This is the \delta_{456}^{231} part
-        if (wick_contraction[1]){
-            for (int rho_i=0; rho_i<Ns; rho_i++){
-            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                auto D1_GAi_ar_ac = D1_GAi()(alpha_f,rho_i)(a_f,c_i);
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                  auto GBf_D2_GBi_ab_ba = GBf_D2_GBi ()(alpha_f,beta_i)(b_f,a_i);
-                for (int rho_f=0; rho_f<Ns; rho_f++){
-                    result()(rho_f,rho_i)() += ee  * D1_GAi_ar_ac
+          }}
+        }}
+      }   
+      //This is the \delta_{456}^{231} part
+      if (wick_contraction[1]){
+        for (int rho_i=0; rho_i<Ns; rho_i++){
+        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          auto D1_GAi_ar_ac = D1_GAi()(alpha_f,rho_i)(a_f,c_i);
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            auto GBf_D2_GBi_ab_ba = GBf_D2_GBi ()(alpha_f,beta_i)(b_f,a_i);
+            for (int rho_f=0; rho_f<Ns; rho_f++){
+              result()(rho_f,rho_i)() += ee  * D1_GAi_ar_ac
                                         * GBf_D2_GBi_ab_ba
                                         * GAf_D3        ()(rho_f,beta_i)(c_f,b_i);
-                }}
-            }}
-        }   
-        //This is the \delta_{456}^{312} part
-        if (wick_contraction[2]){
-            for (int rho_i=0; rho_i<Ns; rho_i++){
-            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                auto GBf_D1_GAi_ar_bc = GBf_D1_GAi()(alpha_f,rho_i)(b_f,c_i);
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                  auto D3_ab_ab = D3 ()(alpha_f,beta_i)(a_f,b_i);
-                for (int rho_f=0; rho_f<Ns; rho_f++){
-                    result()(rho_f,rho_i)() += ee  * GBf_D1_GAi_ar_bc
+            }
+          }
+        }}
+      }   
+      //This is the \delta_{456}^{312} part
+      if (wick_contraction[2]){
+        for (int rho_i=0; rho_i<Ns; rho_i++){
+        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          auto GBf_D1_GAi_ar_bc = GBf_D1_GAi()(alpha_f,rho_i)(b_f,c_i);
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            auto D3_ab_ab = D3 ()(alpha_f,beta_i)(a_f,b_i);
+            for (int rho_f=0; rho_f<Ns; rho_f++){
+              result()(rho_f,rho_i)() += ee  * GBf_D1_GAi_ar_bc
                                         * GAf_D2_GBi    ()(rho_f,beta_i)(c_f,a_i)
                                         * D3_ab_ab;
-                }}
-            }}
-        }   
-        //This is the \delta_{456}^{132} part
-        if (wick_contraction[3]){
-            for (int rho_i=0; rho_i<Ns; rho_i++){
-            for (int rho_f=0; rho_f<Ns; rho_f++){
-                auto GAf_D1_GAi_rr_cc = GAf_D1_GAi()(rho_f,rho_i)(c_f,c_i);
-                for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                    result()(rho_f,rho_i)() -= ee  * GAf_D1_GAi_rr_cc
+            }
+          }
+        }}
+      }   
+      //This is the \delta_{456}^{132} part
+      if (wick_contraction[3]){
+        for (int rho_i=0; rho_i<Ns; rho_i++){
+        for (int rho_f=0; rho_f<Ns; rho_f++){
+          auto GAf_D1_GAi_rr_cc = GAf_D1_GAi()(rho_f,rho_i)(c_f,c_i);
+          for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            result()(rho_f,rho_i)() -= ee  * GAf_D1_GAi_rr_cc
                                         * GBf_D2_GBi    ()(alpha_f,beta_i)(b_f,a_i)
                                         * D3            ()(alpha_f,beta_i)(a_f,b_i);
-                }}
-            }}
-        }   
-        //This is the \delta_{456}^{321} part
-        if (wick_contraction[4]){
-            for (int rho_i=0; rho_i<Ns; rho_i++){
-            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                auto GBf_D1_GAi_ar_bc = GBf_D1_GAi()(alpha_f,rho_i)(b_f,c_i);
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                  auto D2_GBi_ab_aa = D2_GBi()(alpha_f,beta_i)(a_f,a_i);
-                for (int rho_f=0; rho_f<Ns; rho_f++){
-                    result()(rho_f,rho_i)() -= ee  * GBf_D1_GAi_ar_bc
+          }}
+        }}
+      }   
+      //This is the \delta_{456}^{321} part
+      if (wick_contraction[4]){
+        for (int rho_i=0; rho_i<Ns; rho_i++){
+        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          auto GBf_D1_GAi_ar_bc = GBf_D1_GAi()(alpha_f,rho_i)(b_f,c_i);
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            auto D2_GBi_ab_aa = D2_GBi()(alpha_f,beta_i)(a_f,a_i);
+            for (int rho_f=0; rho_f<Ns; rho_f++){
+              result()(rho_f,rho_i)() -= ee  * GBf_D1_GAi_ar_bc
                                         * D2_GBi_ab_aa
                                         * GAf_D3    ()(rho_f,beta_i)(c_f,b_i);
-                }}
-            }}
-        }   
-        //This is the \delta_{456}^{213} part
-        if (wick_contraction[5]){
-            for (int rho_i=0; rho_i<Ns; rho_i++){
-            for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-                auto D1_GAi_ar_ac = D1_GAi()(alpha_f,rho_i)(a_f,c_i);
-                for (int beta_i=0; beta_i<Ns; beta_i++){
-                  auto GBf_D3_ab_bb = GBf_D3()(alpha_f,beta_i)(b_f,b_i);
-                for (int rho_f=0; rho_f<Ns; rho_f++){
-                    result()(rho_f,rho_i)() -= ee  * D1_GAi_ar_ac
+            }
+          }
+        }}
+      }   
+      //This is the \delta_{456}^{213} part
+      if (wick_contraction[5]){
+        for (int rho_i=0; rho_i<Ns; rho_i++){
+        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+          auto D1_GAi_ar_ac = D1_GAi()(alpha_f,rho_i)(a_f,c_i);
+          for (int beta_i=0; beta_i<Ns; beta_i++){
+            auto GBf_D3_ab_bb = GBf_D3()(alpha_f,beta_i)(b_f,b_i);
+            for (int rho_f=0; rho_f<Ns; rho_f++){
+              result()(rho_f,rho_i)() -= ee  * D1_GAi_ar_ac
                                         * GAf_D2_GBi    ()(rho_f,beta_i)(c_f,a_i)
                                         * GBf_D3_ab_bb;
-                }}
-            }}
-        }
-    }}
+            }
+          }
+        }}
+      }
+    }
+  }
 }
 
 /* Computes which wick contractions should be performed for a    *
@@ -530,8 +522,7 @@ void BaryonUtils<FImpl>::ContractBaryons(const PropagatorField &q1_left,
   for (int ie=0; ie < 6 ; ie++){
     if(ie==0 or ie==3){
        bytes += grid->oSites() * (4.*sizeof(int) + 4752.*sizeof(vComplex)) * wick_contractions[ie];
-    }
-    else{
+    } else{
        bytes += grid->oSites() * (64.*sizeof(int) + 5184.*sizeof(vComplex)) * wick_contractions[ie];
     }
   }
@@ -653,7 +644,7 @@ template <class mobj, class mobj2, class robj> accelerator_inline
 void BaryonUtils<FImpl>::BaryonGamma3ptGroup1Site(
                         const mobj &Dq1_ti,
                         const mobj2 &Dq2_spec,
-           //             const mobj2 &Dq3_spec,
+                        const mobj2 &Dq3_spec,
                         const mobj &Dq4_tf,
                                 const Gamma GammaJ,
                                 const Gamma GammaBi,
@@ -661,18 +652,14 @@ void BaryonUtils<FImpl>::BaryonGamma3ptGroup1Site(
                         int wick_contraction,
                         robj &result)
 {
-    Gamma g5(Gamma::Algebra::Gamma5); 
-
-//    auto adjD4_g_D1     = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq1_ti;
-    auto adjD4          = g5 * adj(Dq4_tf) * g5 ;
-    auto adjD4_g_D1     = adjD4 * GammaJ * Dq1_ti;
-    auto Gf_adjD4_g_D1  = GammaBf * adjD4_g_D1;
-    auto D2_Gi          = Dq2_spec * GammaBi;
-    auto Gf_D2_Gi       = GammaBf * D2_Gi;
-
-//    auto Gf_D3          = GammaBf * Dq3_spec;  // including a second mobj2 parameter leads to compilation error
-    auto Gf_D3          = GammaBf * Dq2_spec;  //WRONG!!!!!
+  Gamma g5(Gamma::Algebra::Gamma5); 
 
+  auto adjD4          = g5 * adj(Dq4_tf) * g5 ;
+  auto adjD4_g_D1     = adjD4 * GammaJ * Dq1_ti;
+  auto Gf_adjD4_g_D1  = GammaBf * adjD4_g_D1;
+  auto D2_Gi          = Dq2_spec * GammaBi;
+  auto Gf_D2_Gi       = GammaBf * D2_Gi;
+  auto Gf_D3          = GammaBf * Dq3_spec;  
 
   Real ee;
 
@@ -681,65 +668,65 @@ void BaryonUtils<FImpl>::BaryonGamma3ptGroup1Site(
     int b_f    = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b
     int c_f    = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_f = (ie_f < 3 ? 1 : -1);
-  for (int ie_i=0; ie_i < 6 ; ie_i++){
-    int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
-    int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
-    int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
-    int eSgn_i = (ie_i < 3 ? 1 : -1);
+    for (int ie_i=0; ie_i < 6 ; ie_i++){
+      int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
+      int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
+      int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_i = (ie_i < 3 ? 1 : -1);
 
-    ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      ee = Real(eSgn_f * eSgn_i);
 
-        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-        for (int beta_i=0; beta_i<Ns; beta_i++){
-            auto D2_Gi_ab_aa        = D2_Gi     ()(alpha_f,beta_i)(a_f,a_i);
-            auto Gf_D3_ab_bb        = Gf_D3     ()(alpha_f,beta_i)(b_f,b_i);
-            auto Gf_D2_Gi_ab_ba     = Gf_D2_Gi  ()(alpha_f,beta_i)(b_f,a_i);
-            //auto Dq3_spec_ab_ab     = Dq3_spec  ()(alpha_f,beta_i)(a_f,b_i);
-            auto Dq3_spec_ab_ab     = Dq2_spec  ()(alpha_f,beta_i)(a_f,b_i); // WRONG!!!!
+      for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+      for (int beta_i=0; beta_i<Ns; beta_i++){
+        auto D2_Gi_ab_aa        = D2_Gi     ()(alpha_f,beta_i)(a_f,a_i);
+        auto Gf_D3_ab_bb        = Gf_D3     ()(alpha_f,beta_i)(b_f,b_i);
+        auto Gf_D2_Gi_ab_ba     = Gf_D2_Gi  ()(alpha_f,beta_i)(b_f,a_i);
+        auto Dq3_spec_ab_ab     = Dq3_spec  ()(alpha_f,beta_i)(a_f,b_i);
 
-            for (int gamma_i=0; gamma_i<Ns; gamma_i++){
-                auto ee_adjD4_g_D1_ag_ac        = ee * adjD4_g_D1   ()(alpha_f,gamma_i)(a_f,c_i);
-                auto ee_Gf_adjD4_g_D1_ag_bc     = ee * Gf_adjD4_g_D1()(alpha_f,gamma_i)(b_f,c_i);
-            for (int gamma_f=0; gamma_f<Ns; gamma_f++){
-                auto ee_adjD4_g_D1_gg_cc        = ee * adjD4_g_D1   ()(gamma_f,gamma_i)(c_f,c_i);
-                //auto Dq3_spec_gb_cb             = Dq3_spec          ()(gamma_f,beta_i)(c_f,b_i);
-                auto Dq3_spec_gb_cb             = Dq2_spec          ()(gamma_f,beta_i)(c_f,b_i); //WRONG!!!!
-                auto D2_Gi_gb_ca                = D2_Gi             ()(gamma_f,beta_i)(c_f,a_i);
+        for (int gamma_i=0; gamma_i<Ns; gamma_i++){
+          auto ee_adjD4_g_D1_ag_ac        = ee * adjD4_g_D1   ()(alpha_f,gamma_i)(a_f,c_i);
+          auto ee_Gf_adjD4_g_D1_ag_bc     = ee * Gf_adjD4_g_D1()(alpha_f,gamma_i)(b_f,c_i);
+          for (int gamma_f=0; gamma_f<Ns; gamma_f++){
+            auto ee_adjD4_g_D1_gg_cc        = ee * adjD4_g_D1   ()(gamma_f,gamma_i)(c_f,c_i);
+            auto Dq3_spec_gb_cb             = Dq3_spec          ()(gamma_f,beta_i)(c_f,b_i);
+            auto D2_Gi_gb_ca                = D2_Gi             ()(gamma_f,beta_i)(c_f,a_i);
 
 
-                if(wick_contraction == 1) { // Do contraction I1
-                    result()(gamma_f,gamma_i)() -= ee_adjD4_g_D1_gg_cc
+            if(wick_contraction == 1) { // Do contraction I1
+              result()(gamma_f,gamma_i)() -= ee_adjD4_g_D1_gg_cc
                                                         * D2_Gi_ab_aa
                                                         * Gf_D3_ab_bb;
-                }
-                if(wick_contraction == 2) { // Do contraction I2
-                    result()(gamma_f,gamma_i)() -= ee_adjD4_g_D1_ag_ac
+            }
+            if(wick_contraction == 2) { // Do contraction I2
+              result()(gamma_f,gamma_i)() -= ee_adjD4_g_D1_ag_ac
                                                         * Gf_D2_Gi_ab_ba
                                                         * Dq3_spec_gb_cb;
-                }
-                if(wick_contraction == 3) { // Do contraction I3
-                    result()(gamma_f,gamma_i)() -= ee_Gf_adjD4_g_D1_ag_bc
+            }
+            if(wick_contraction == 3) { // Do contraction I3
+              result()(gamma_f,gamma_i)() -= ee_Gf_adjD4_g_D1_ag_bc
                                                         * D2_Gi_gb_ca
                                                         * Dq3_spec_ab_ab;
-                }
-                if(wick_contraction == 4) { // Do contraction I4
-                    result()(gamma_f,gamma_i)() += ee_adjD4_g_D1_gg_cc
+            }
+            if(wick_contraction == 4) { // Do contraction I4
+              result()(gamma_f,gamma_i)() += ee_adjD4_g_D1_gg_cc
                                                         * Gf_D2_Gi_ab_ba
                                                         * Dq3_spec_ab_ab;
-                }
-                if(wick_contraction == 5) { // Do contraction I5
-                    result()(gamma_f,gamma_i)() += ee_Gf_adjD4_g_D1_ag_bc
+            }
+            if(wick_contraction == 5) { // Do contraction I5
+              result()(gamma_f,gamma_i)() += ee_Gf_adjD4_g_D1_ag_bc
                                                         * D2_Gi_ab_aa
                                                         * Dq3_spec_gb_cb;
-                }
-                if(wick_contraction == 6) { // Do contraction I6
-                    result()(gamma_f,gamma_i)() += ee_adjD4_g_D1_ag_ac
+            }
+            if(wick_contraction == 6) { // Do contraction I6
+              result()(gamma_f,gamma_i)() += ee_adjD4_g_D1_ag_ac
                                                         * D2_Gi_gb_ca
                                                         * Gf_D3_ab_bb;
-                }
-            }}
-        }}
-    }}
+            }
+          }
+        }
+      }}
+    }
+  }
 }
 
 /* Dq1_spec is a quark line from t_i to t_f
@@ -751,7 +738,7 @@ template <class mobj, class mobj2, class robj> accelerator_inline
 void BaryonUtils<FImpl>::BaryonGamma3ptGroup2Site(
                         const mobj2 &Dq1_spec,
                         const mobj &Dq2_ti,
-                    //    const mobj2 &Dq3_spec,
+                        const mobj2 &Dq3_spec,
                         const mobj &Dq4_tf,
                                 const Gamma GammaJ,
                                 const Gamma GammaBi,
@@ -759,14 +746,12 @@ void BaryonUtils<FImpl>::BaryonGamma3ptGroup2Site(
                         int wick_contraction,
                         robj &result)
 {
-    Gamma g5(Gamma::Algebra::Gamma5); 
-
-    auto adjD4_g_D2_Gi      = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi;
-    auto Gf_adjD4_g_D2_Gi   = GammaBf * adjD4_g_D2_Gi;
-    auto Gf_D1              = GammaBf * Dq1_spec;
-    //auto Gf_D3              = GammaBf * Dq3_spec;
-    auto Gf_D3              = GammaBf * Dq1_spec; // WRONG!!!!!
+  Gamma g5(Gamma::Algebra::Gamma5); 
 
+  auto adjD4_g_D2_Gi      = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq2_ti * GammaBi;
+  auto Gf_adjD4_g_D2_Gi   = GammaBf * adjD4_g_D2_Gi;
+  auto Gf_D1              = GammaBf * Dq1_spec;
+  auto Gf_D3              = GammaBf * Dq3_spec;
 
   Real ee;
 
@@ -775,64 +760,64 @@ void BaryonUtils<FImpl>::BaryonGamma3ptGroup2Site(
     int b_f    = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b
     int c_f    = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_f = (ie_f < 3 ? 1 : -1);
-  for (int ie_i=0; ie_i < 6 ; ie_i++){
-    int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
-    int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
-    int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
-    int eSgn_i = (ie_i < 3 ? 1 : -1);
+    for (int ie_i=0; ie_i < 6 ; ie_i++){
+      int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
+      int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
+      int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_i = (ie_i < 3 ? 1 : -1);
 
-    ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
 
-        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-        for (int beta_i=0; beta_i<Ns; beta_i++){
-            auto adjD4_g_D2_Gi_ab_aa        = adjD4_g_D2_Gi     ()(alpha_f,beta_i)(a_f,a_i);
-            auto Gf_D3_ab_bb                = Gf_D3             ()(alpha_f,beta_i)(b_f,b_i);
-            auto Gf_adjD4_g_D2_Gi_ab_ba     = Gf_adjD4_g_D2_Gi  ()(alpha_f,beta_i)(b_f,a_i);
-            //auto Dq3_spec_ab_ab             = Dq3_spec          ()(alpha_f,beta_i)(a_f,b_i);
-            auto Dq3_spec_ab_ab             = Dq1_spec          ()(alpha_f,beta_i)(a_f,b_i); //WRONG!!!
+      for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+      for (int beta_i=0; beta_i<Ns; beta_i++){
+        auto adjD4_g_D2_Gi_ab_aa        = adjD4_g_D2_Gi     ()(alpha_f,beta_i)(a_f,a_i);
+        auto Gf_D3_ab_bb                = Gf_D3             ()(alpha_f,beta_i)(b_f,b_i);
+        auto Gf_adjD4_g_D2_Gi_ab_ba     = Gf_adjD4_g_D2_Gi  ()(alpha_f,beta_i)(b_f,a_i);
+        auto Dq3_spec_ab_ab             = Dq3_spec          ()(alpha_f,beta_i)(a_f,b_i);
 
-            for (int gamma_i=0; gamma_i<Ns; gamma_i++){ 
-                auto ee_Dq1_spec_ag_ac      = ee * Dq1_spec     ()(alpha_f,gamma_i)(a_f,c_i);
-                auto ee_Gf_D1_ag_bc         = ee * Gf_D1        ()(alpha_f,gamma_i)(b_f,c_i);
-            for (int gamma_f=0; gamma_f<Ns; gamma_f++){
-                auto ee_Dq1_spec_gg_cc      = ee * Dq1_spec     ()(gamma_f,gamma_i)(c_f,c_i);
-                //auto Dq3_spec_gb_cb         = Dq3_spec          ()(gamma_f,beta_i)(c_f,b_i);
-                auto Dq3_spec_gb_cb         = Dq1_spec          ()(gamma_f,beta_i)(c_f,b_i); //WRONG!!!
-                auto adjD4_g_D2_Gi_gb_ca    = adjD4_g_D2_Gi     ()(gamma_f,beta_i)(c_f,a_i);
+        for (int gamma_i=0; gamma_i<Ns; gamma_i++){ 
+          auto ee_Dq1_spec_ag_ac      = ee * Dq1_spec     ()(alpha_f,gamma_i)(a_f,c_i);
+          auto ee_Gf_D1_ag_bc         = ee * Gf_D1        ()(alpha_f,gamma_i)(b_f,c_i);
+          for (int gamma_f=0; gamma_f<Ns; gamma_f++){
+            auto ee_Dq1_spec_gg_cc      = ee * Dq1_spec     ()(gamma_f,gamma_i)(c_f,c_i);
+            auto Dq3_spec_gb_cb         = Dq3_spec          ()(gamma_f,beta_i)(c_f,b_i);
+            auto adjD4_g_D2_Gi_gb_ca    = adjD4_g_D2_Gi     ()(gamma_f,beta_i)(c_f,a_i);
 
-                if(wick_contraction == 1) { // Do contraction II1
-                    result()(gamma_f,gamma_i)() -= ee_Dq1_spec_gg_cc
+            if(wick_contraction == 1) { // Do contraction II1
+              result()(gamma_f,gamma_i)() -= ee_Dq1_spec_gg_cc
                                                         * adjD4_g_D2_Gi_ab_aa
                                                         * Gf_D3_ab_bb;
-                }
-                if(wick_contraction == 2) { // Do contraction II2
-                    result()(gamma_f,gamma_i)() -= ee_Dq1_spec_ag_ac
+            }
+            if(wick_contraction == 2) { // Do contraction II2
+              result()(gamma_f,gamma_i)() -= ee_Dq1_spec_ag_ac
                                                         * Gf_adjD4_g_D2_Gi_ab_ba
                                                         * Dq3_spec_gb_cb;
-                }
-                if(wick_contraction == 3) { // Do contraction II3
-                    result()(gamma_f,gamma_i)() -= ee_Gf_D1_ag_bc
+            }
+            if(wick_contraction == 3) { // Do contraction II3
+              result()(gamma_f,gamma_i)() -= ee_Gf_D1_ag_bc
                                                         * adjD4_g_D2_Gi_gb_ca
                                                         * Dq3_spec_ab_ab;
-                }
-                if(wick_contraction == 4) { // Do contraction II4
-                    result()(gamma_f,gamma_i)() += ee_Dq1_spec_gg_cc
+            }
+            if(wick_contraction == 4) { // Do contraction II4
+              result()(gamma_f,gamma_i)() += ee_Dq1_spec_gg_cc
                                                         * Gf_adjD4_g_D2_Gi_ab_ba
                                                         * Dq3_spec_ab_ab;
-                }
-                if(wick_contraction == 5) { // Do contraction II5
-                    result()(gamma_f,gamma_i)() += ee_Gf_D1_ag_bc
+            }
+            if(wick_contraction == 5) { // Do contraction II5
+              result()(gamma_f,gamma_i)() += ee_Gf_D1_ag_bc
                                                         * adjD4_g_D2_Gi_ab_aa
                                                         * Dq3_spec_gb_cb;
-                }
-                if(wick_contraction == 6) { // Do contraction II6
-                    result()(gamma_f,gamma_i)() += ee_Dq1_spec_ag_ac
+            }
+            if(wick_contraction == 6) { // Do contraction II6
+              result()(gamma_f,gamma_i)() += ee_Dq1_spec_ag_ac
                                                         * adjD4_g_D2_Gi_gb_ca
                                                         * Gf_D3_ab_bb;
-                }
-            }}
-        }}
-    }}
+            }
+          }
+        }
+      }}
+    }
+  }
 }
 
 /* Dq1_spec is a quark line from t_i to t_f
@@ -843,7 +828,7 @@ template<class FImpl>
 template <class mobj, class mobj2, class robj> accelerator_inline
 void BaryonUtils<FImpl>::BaryonGamma3ptGroup3Site(
                         const mobj2 &Dq1_spec,
-                       // const mobj2 &Dq2_spec,
+                        const mobj2 &Dq2_spec,
                         const mobj &Dq3_ti,
                         const mobj &Dq4_tf,
                                 const Gamma GammaJ,
@@ -852,15 +837,13 @@ void BaryonUtils<FImpl>::BaryonGamma3ptGroup3Site(
                         int wick_contraction,
                         robj &result)
 {
-    Gamma g5(Gamma::Algebra::Gamma5);
-
-    auto adjD4_g_D3     = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti;
-    auto Gf_adjD4_g_D3  = GammaBf * adjD4_g_D3;
-    auto Gf_D1          = GammaBf * Dq1_spec;
-    //auto D2_Gi          = Dq2_spec * GammaBi;
-    auto D2_Gi          = Dq1_spec * GammaBi; //WRONG!!!!!!!!!!!!!!!!!
-    auto Gf_D2_Gi       = GammaBf * D2_Gi;
+  Gamma g5(Gamma::Algebra::Gamma5);
 
+  auto adjD4_g_D3     = g5 * adj(Dq4_tf) * g5 * GammaJ * Dq3_ti;
+  auto Gf_adjD4_g_D3  = GammaBf * adjD4_g_D3;
+  auto Gf_D1          = GammaBf * Dq1_spec;
+  auto D2_Gi          = Dq2_spec * GammaBi;
+  auto Gf_D2_Gi       = GammaBf * D2_Gi;
 
   Real ee;
 
@@ -869,62 +852,64 @@ void BaryonUtils<FImpl>::BaryonGamma3ptGroup3Site(
     int b_f    = (ie_f < 3 ? (ie_f+1)%3 : (8-ie_f)%3 ); //epsilon[ie_n][1]; //b
     int c_f    = (ie_f < 3 ? (ie_f+2)%3 : (7-ie_f)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_f = (ie_f < 3 ? 1 : -1);
-  for (int ie_i=0; ie_i < 6 ; ie_i++){
-    int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
-    int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
-    int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
-    int eSgn_i = (ie_i < 3 ? 1 : -1);
+    for (int ie_i=0; ie_i < 6 ; ie_i++){
+      int a_i = (ie_i < 3 ? ie_i       : (6-ie_i)%3 ); //epsilon[ie_s][0]; //a'
+      int b_i = (ie_i < 3 ? (ie_i+1)%3 : (8-ie_i)%3 ); //epsilon[ie_s][1]; //b'
+      int c_i = (ie_i < 3 ? (ie_i+2)%3 : (7-ie_i)%3 ); //epsilon[ie_s][2]; //c'
+      int eSgn_i = (ie_i < 3 ? 1 : -1);
 
-    ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      ee = Real(eSgn_f * eSgn_i); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
 
-        for (int alpha_f=0; alpha_f<Ns; alpha_f++){
-        for (int beta_i=0; beta_i<Ns; beta_i++){
-            auto D2_Gi_ab_aa            = D2_Gi         ()(alpha_f,beta_i)(a_f,a_i);
-            auto Gf_adjD4_g_D3_ab_bb    = Gf_adjD4_g_D3 ()(alpha_f,beta_i)(b_f,b_i);
-            auto Gf_D2_Gi_ab_ba         = Gf_D2_Gi      ()(alpha_f,beta_i)(b_f,a_i);
-            auto adjD4_g_D3_ab_ab       = adjD4_g_D3    ()(alpha_f,beta_i)(a_f,b_i);
+      for (int alpha_f=0; alpha_f<Ns; alpha_f++){
+      for (int beta_i=0; beta_i<Ns; beta_i++){
+        auto D2_Gi_ab_aa            = D2_Gi         ()(alpha_f,beta_i)(a_f,a_i);
+        auto Gf_adjD4_g_D3_ab_bb    = Gf_adjD4_g_D3 ()(alpha_f,beta_i)(b_f,b_i);
+        auto Gf_D2_Gi_ab_ba         = Gf_D2_Gi      ()(alpha_f,beta_i)(b_f,a_i);
+        auto adjD4_g_D3_ab_ab       = adjD4_g_D3    ()(alpha_f,beta_i)(a_f,b_i);
 
-            for (int gamma_i=0; gamma_i<Ns; gamma_i++) {
-                auto ee_Dq1_spec_ag_ac  = ee * Dq1_spec ()(alpha_f,gamma_i)(a_f,c_i);
-                auto ee_Gf_D1_ag_bc     = ee * Gf_D1    ()(alpha_f,gamma_i)(b_f,c_i);
-            for (int gamma_f=0; gamma_f<Ns; gamma_f++) {
-                auto ee_Dq1_spec_gg_cc  = ee * Dq1_spec ()(gamma_f,gamma_i)(c_f,c_i);
-                auto adjD4_g_D3_gb_cb   = adjD4_g_D3    ()(gamma_f,beta_i)(c_f,b_i);
-                auto D2_Gi_gb_ca        = D2_Gi         ()(gamma_f,beta_i)(c_f,a_i);
+        for (int gamma_i=0; gamma_i<Ns; gamma_i++) {
+          auto ee_Dq1_spec_ag_ac  = ee * Dq1_spec ()(alpha_f,gamma_i)(a_f,c_i);
+          auto ee_Gf_D1_ag_bc     = ee * Gf_D1    ()(alpha_f,gamma_i)(b_f,c_i);
+          for (int gamma_f=0; gamma_f<Ns; gamma_f++) {
+            auto ee_Dq1_spec_gg_cc  = ee * Dq1_spec ()(gamma_f,gamma_i)(c_f,c_i);
+            auto adjD4_g_D3_gb_cb   = adjD4_g_D3    ()(gamma_f,beta_i)(c_f,b_i);
+            auto D2_Gi_gb_ca        = D2_Gi         ()(gamma_f,beta_i)(c_f,a_i);
 
-                if(wick_contraction == 1) { // Do contraction III1
-                    result()(gamma_f,gamma_i)() -= ee_Dq1_spec_gg_cc
+            if(wick_contraction == 1) { // Do contraction III1
+              result()(gamma_f,gamma_i)() -= ee_Dq1_spec_gg_cc
                                                         * D2_Gi_ab_aa
                                                         * Gf_adjD4_g_D3_ab_bb;
-                }
-                if(wick_contraction == 2) { // Do contraction III2
-                    result()(gamma_f,gamma_i)() -= ee_Dq1_spec_ag_ac
+            }
+            if(wick_contraction == 2) { // Do contraction III2
+              result()(gamma_f,gamma_i)() -= ee_Dq1_spec_ag_ac
                                                         * Gf_D2_Gi_ab_ba
                                                         * adjD4_g_D3_gb_cb;
-                }
-                if(wick_contraction == 3) { // Do contraction III3
-                    result()(gamma_f,gamma_i)() -= ee_Gf_D1_ag_bc
+            }
+            if(wick_contraction == 3) { // Do contraction III3
+              result()(gamma_f,gamma_i)() -= ee_Gf_D1_ag_bc
                                                         * D2_Gi_gb_ca
                                                         * adjD4_g_D3_ab_ab;
-                }
-                if(wick_contraction == 4) { // Do contraction III4
-                    result()(gamma_f,gamma_i)() += ee_Dq1_spec_gg_cc
+            }
+            if(wick_contraction == 4) { // Do contraction III4
+              result()(gamma_f,gamma_i)() += ee_Dq1_spec_gg_cc
                                                         * Gf_D2_Gi_ab_ba
                                                         * adjD4_g_D3_ab_ab;
-                }
-                if(wick_contraction == 5) { // Do contraction III5
-                    result()(gamma_f,gamma_i)() += ee_Gf_D1_ag_bc
+            }
+            if(wick_contraction == 5) { // Do contraction III5
+              result()(gamma_f,gamma_i)() += ee_Gf_D1_ag_bc
                                                         * D2_Gi_ab_aa
                                                         * adjD4_g_D3_gb_cb;
-                }
-                if(wick_contraction == 6) { // Do contraction III6
-                    result()(gamma_f,gamma_i)() += ee_Dq1_spec_ag_ac
+            }
+            if(wick_contraction == 6) { // Do contraction III6
+              result()(gamma_f,gamma_i)() += ee_Dq1_spec_ag_ac
                                                         * D2_Gi_gb_ca
                                                         * Gf_adjD4_g_D3_ab_bb;
-                }
-            }}
-        }}
-    }}
+            }
+          }
+        }
+      }}
+    }
+  }
 }
 
 /* The group indicates which inital state quarks the current is  * 
@@ -949,81 +934,44 @@ void BaryonUtils<FImpl>::BaryonGamma3pt(
   assert(Ns==4 && "Baryon code only implemented for N_spin = 4");
   assert(Nc==3 && "Baryon code only implemented for N_colour = 3");
 
-    GridBase *grid = q_tf.Grid();
+  GridBase *grid = q_tf.Grid();
 
-    // autoView( vcorr, stn_corr, CpuWrite);
-    // autoView( vq_ti , q_ti, CpuRead);
-    // autoView( vq_tf , q_tf, CpuRead);
+  autoView( vcorr , stn_corr , AcceleratorWrite);
+  autoView( vq_ti , q_ti     , AcceleratorRead);
+  autoView( vq_tf , q_tf     , AcceleratorRead);
 
-    // if (group == 1) {
-    //     accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
-    //         auto Dq_ti = vq_ti[ss];
-    //         auto Dq_tf = vq_tf[ss];
-    //         sobj result=Zero();
-    //         BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
-    //         vcorr[ss] += result; 
-    //     });//end loop over lattice sites
-    // } else if (group == 2) {
-    //     accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
-    //         auto Dq_ti = vq_ti[ss];
-    //         auto Dq_tf = vq_tf[ss];
-    //         sobj result=Zero();
-    //         BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
-    //         vcorr[ss] += result; 
-    //     });//end loop over lattice sites
-    // } else if (group == 3) {
-    //     accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
-    //         auto Dq_ti = vq_ti[ss];
-    //         auto Dq_tf = vq_tf[ss];
-    //         sobj result=Zero();
-    //         BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
+  Vector<mobj> my_Dq_spec{Dq_spec1,Dq_spec2};
+  mobj * Dq_spec_p = &my_Dq_spec[0];
 
-    //         vcorr[ss] += result; 
-    //     });//end loop over lattice sites
-    // }
+  if (group == 1) {
+    accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+      auto Dq_ti = vq_ti(ss);
+      auto Dq_tf = vq_tf(ss);
+      typedef decltype(coalescedRead(vcorr[0])) spinor;
+      spinor result=Zero();
+      BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec_p[0],Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
+      coalescedWrite(vcorr[ss],result);
+    });//end loop over lattice sites
 
-    autoView( vcorr , stn_corr , AcceleratorWrite);
-    autoView( vq_ti , q_ti     , AcceleratorRead);
-    autoView( vq_tf , q_tf     , AcceleratorRead);
-
-    if (group == 1) {
-        accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
-            auto Dq_ti = vq_ti(ss);
-            auto Dq_tf = vq_tf(ss);
-            //sobj result=Zero();
-            typedef decltype(coalescedRead(vcorr[0])) spinor;
-            spinor result=Zero();
-            //BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
-            BaryonGamma3ptGroup1Site(Dq_ti,Dq_spec1,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG
-            // vcorr[ss] += result; 
-            coalescedWrite(vcorr[ss],result);
-        });//end loop over lattice sites
-
-    } else if (group == 2) {
-        accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
-            auto Dq_ti = vq_ti(ss);
-            auto Dq_tf = vq_tf(ss);
-            //sobj result=Zero();
-            typedef decltype(coalescedRead(vcorr[0])) spinor;
-            spinor result=Zero();
-           // BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_spec2,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
-            BaryonGamma3ptGroup2Site(Dq_spec1,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG
-            // vcorr[ss] += result; 
-            coalescedWrite(vcorr[ss],result);
-        });//end loop over lattice sites
-    } else if (group == 3) {
-        accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
-            auto Dq_ti = vq_ti(ss);
-            auto Dq_tf = vq_tf(ss);
-            //sobj result=Zero();
-            typedef decltype(coalescedRead(vcorr[0])) spinor;
-            spinor result=Zero();
-            //BaryonGamma3ptGroup3Site(Dq_spec1,Dq_spec2,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result);
-            BaryonGamma3ptGroup3Site(Dq_spec1,Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); //WRONG
-            // vcorr[ss] += result; 
-            coalescedWrite(vcorr[ss],result);
-        });//end loop over lattice sites
-    }
+  } else if (group == 2) {
+    accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+      auto Dq_ti = vq_ti(ss);
+      auto Dq_tf = vq_tf(ss);
+      typedef decltype(coalescedRead(vcorr[0])) spinor;
+      spinor result=Zero();
+      BaryonGamma3ptGroup2Site(Dq_spec_p[0],Dq_ti,Dq_spec_p[1],Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); 
+      coalescedWrite(vcorr[ss],result);
+    });//end loop over lattice sites
+  } else if (group == 3) {
+    accelerator_for(ss, grid->oSites(), grid->Nsimd(), {
+      auto Dq_ti = vq_ti(ss);
+      auto Dq_tf = vq_tf(ss);
+      typedef decltype(coalescedRead(vcorr[0])) spinor;
+      spinor result=Zero();
+      BaryonGamma3ptGroup3Site(Dq_spec_p[0],Dq_spec_p[1],Dq_ti,Dq_tf,GammaJ,GammaBi,GammaBf,wick_contraction,result); 
+      coalescedWrite(vcorr[ss],result);
+    });//end loop over lattice sites
+  }
 
 }
 
@@ -1052,7 +1000,6 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
 
   Gamma g5(Gamma::Algebra::Gamma5); 
 
-  //auto Gn_adjDd_GH_Ds   = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti;
   auto adjDd_GH_Ds      = g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti;
   auto Gn_adjDd_GH_Ds   = GammaB_nucl * adjDd_GH_Ds;
   auto Du_Gs            = Du_spec * GammaB_sigma;
@@ -1066,33 +1013,33 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1EyeSite(const mobj &Dq_loop,
     int b_n    = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b
     int c_n    = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_n = (ie_n < 3 ? 1 : -1);
-  for (int ie_s=0; ie_s < 6 ; ie_s++){
-    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_s=0; ie_s < 6 ; ie_s++){
+      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);
 
-    ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      ee = Real(eSgn_n * eSgn_s); 
+      for (int alpha_n=0; alpha_n<Ns; alpha_n++){
+      for (int beta_s=0;  beta_s<Ns;  beta_s++){
 
-    for (int alpha_n=0; alpha_n<Ns; alpha_n++){
-    for (int beta_s=0;  beta_s<Ns;  beta_s++){
+        auto Gn_adjDd_GH_Ds_ab_bb = Gn_adjDd_GH_Ds ()(alpha_n, beta_s)(b_n,b_s);
 
-      auto Gn_adjDd_GH_Ds_ab_bb = Gn_adjDd_GH_Ds ()(alpha_n, beta_s)(b_n,b_s);
-
-      for (int gamma_s=0; gamma_s<Ns; gamma_s++){
-      for (int gamma_n=0; gamma_n<Ns; gamma_n++){
-        result()(gamma_n,gamma_s)() += ee   * Gn_adjDd_GH_Ds_ab_bb
+        for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+        for (int gamma_n=0; gamma_n<Ns; gamma_n++){
+          result()(gamma_n,gamma_s)() += ee   * Gn_adjDd_GH_Ds_ab_bb
                                             * Du_spec         ()(gamma_n,gamma_s)(c_n,c_s)
                                             * Du_Gs           ()(alpha_n, beta_s)(a_n,a_s) 
                                             * Tr_Dq_GH;
 
-        result()(gamma_n,gamma_s)() -= ee   * Gn_adjDd_GH_Ds_ab_bb
+          result()(gamma_n,gamma_s)() -= ee   * Gn_adjDd_GH_Ds_ab_bb
                                             * Du_spec         ()(alpha_n,gamma_s)(a_n,c_s)
                                             * Du_Gs           ()(gamma_n, beta_s)(c_n,a_s) 
                                             * Tr_Dq_GH;
+        }}
       }}
-    }}
-  }}
+    }
+  }
 }
 
 /* Du_ti is a quark line from t_i to t_H
@@ -1116,7 +1063,6 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
   Gamma g5(Gamma::Algebra::Gamma5); 
 
   auto Du_Gs          = Du_spec * GammaB_sigma;
-  //auto Gn_adjDd_GH_Ds = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti;
   auto adjDd_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Ds_ti;
   auto Gn_adjDd_GH_Ds = GammaB_nucl * adjDd_GH_Ds;
   auto adjDu_GH_Du    = g5 * adj(Du_tf) * g5 * Gamma_H * Du_ti;
@@ -1129,40 +1075,41 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ1NonEyeSite(const mobj &Du_ti,
     int b_n    = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b
     int c_n    = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_n = (ie_n < 3 ? 1 : -1);
-  for (int ie_s=0; ie_s < 6 ; ie_s++){
-    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_s=0; ie_s < 6 ; ie_s++){
+      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);
 
-    ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
 
-    for (int alpha_n=0; alpha_n<Ns; alpha_n++){
-    for (int beta_s=0;  beta_s<Ns;  beta_s++){
+      for (int alpha_n=0; alpha_n<Ns; alpha_n++){
+      for (int beta_s=0;  beta_s<Ns;  beta_s++){
 
-      auto Gn_adjDd_GH_Ds_ab_bb = Gn_adjDd_GH_Ds  ()(alpha_n, beta_s)(b_n,b_s);
+        auto Gn_adjDd_GH_Ds_ab_bb = Gn_adjDd_GH_Ds  ()(alpha_n, beta_s)(b_n,b_s);
 
-      for (int gamma_s=0; gamma_s<Ns; gamma_s++){
-      for (int gamma_n=0; gamma_n<Ns; gamma_n++){
+        for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+        for (int gamma_n=0; gamma_n<Ns; gamma_n++){
 
-        result()(gamma_n,gamma_s)() += ee * Gn_adjDd_GH_Ds_ab_bb
+          result()(gamma_n,gamma_s)() += ee * Gn_adjDd_GH_Ds_ab_bb
                                           * adjDu_GH_Du     ()(alpha_n,gamma_s)(a_n,c_s)
                                           * Du_Gs           ()(gamma_n, beta_s)(c_n,a_s);
 
-        result()(gamma_n,gamma_s)() += ee * Gn_adjDd_GH_Ds_ab_bb
+          result()(gamma_n,gamma_s)() += ee * Gn_adjDd_GH_Ds_ab_bb
                                           * adjDu_GH_Du_Gs  ()(gamma_n, beta_s)(c_n,a_s)
                                           * Du_spec         ()(alpha_n,gamma_s)(a_n,c_s);
 
-        result()(gamma_n,gamma_s)() -= ee * Gn_adjDd_GH_Ds_ab_bb
+          result()(gamma_n,gamma_s)() -= ee * Gn_adjDd_GH_Ds_ab_bb
                                           * adjDu_GH_Du_Gs  ()(alpha_n, beta_s)(a_n,a_s)
                                           * Du_spec         ()(gamma_n,gamma_s)(c_n,c_s);
 
-        result()(gamma_n,gamma_s)() -= ee * Gn_adjDd_GH_Ds_ab_bb
+          result()(gamma_n,gamma_s)() -= ee * Gn_adjDd_GH_Ds_ab_bb
                                           * adjDu_GH_Du     ()(gamma_n,gamma_s)(c_n,c_s)
                                           * Du_Gs           ()(alpha_n, beta_s)(a_n,a_s);
+        }}
       }}
-    }}
-  }}
+    }
+  }
 }
 
 //Equivalent to "One-trace"
@@ -1184,7 +1131,6 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
 
   Gamma g5(Gamma::Algebra::Gamma5); 
 
-  //auto Gn_adjDd_GH_Duloop_GH_Ds = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Dq_loop * Gamma_H * Ds_ti;
   auto adjDd_GH_Duloop_GH_Ds = g5 * adj(Dd_tf) * g5 * Gamma_H * Dq_loop * Gamma_H * Ds_ti;
   auto Gn_adjDd_GH_Duloop_GH_Ds = GammaB_nucl * adjDd_GH_Duloop_GH_Ds;
   auto Du_Gs = Du_spec * GammaB_sigma;
@@ -1196,32 +1142,33 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2EyeSite(const mobj &Dq_loop,
     int b_n    = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b
     int c_n    = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_n = (ie_n < 3 ? 1 : -1);
-  for (int ie_s=0; ie_s < 6 ; ie_s++){
-    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_s=0; ie_s < 6 ; ie_s++){
+      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);
 
-    ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
 
-    for (int alpha_n=0; alpha_n<Ns; alpha_n++){
-    for (int beta_s=0; beta_s<Ns; beta_s++){
+      for (int alpha_n=0; alpha_n<Ns; alpha_n++){
+      for (int beta_s=0; beta_s<Ns; beta_s++){
 
-      auto Gn_adjDd_GH_Duloop_GH_Ds_ab_bb = Gn_adjDd_GH_Duloop_GH_Ds ()(alpha_n,beta_s)(b_n,b_s);
+        auto Gn_adjDd_GH_Duloop_GH_Ds_ab_bb = Gn_adjDd_GH_Duloop_GH_Ds ()(alpha_n,beta_s)(b_n,b_s);
 
-      for (int gamma_s=0; gamma_s<Ns; gamma_s++){
-      for (int gamma_n=0; gamma_n<Ns; gamma_n++){
+        for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+        for (int gamma_n=0; gamma_n<Ns; gamma_n++){
 
-        result()(gamma_n,gamma_s)() -= ee   * Du_spec         ()(gamma_n,gamma_s)(c_n,c_s)
+          result()(gamma_n,gamma_s)() -= ee   * Du_spec         ()(gamma_n,gamma_s)(c_n,c_s)
                                             * Du_Gs           ()(alpha_n,beta_s)(a_n,a_s) 
                                             * Gn_adjDd_GH_Duloop_GH_Ds_ab_bb;
 
-        result()(gamma_n,gamma_s)() += ee   * Du_Gs         ()(alpha_n,gamma_s)(a_n,c_s)
+          result()(gamma_n,gamma_s)() += ee   * Du_Gs         ()(alpha_n,gamma_s)(a_n,c_s)
                                             * Du_spec       ()(gamma_n,beta_s)(c_n,a_s) 
                                             * Gn_adjDd_GH_Duloop_GH_Ds_ab_bb;
+        }}
       }}
-    }}
-  }}
+    }
+  }
 }
 
 /* Du_ti is a quark line from t_i to t_H
@@ -1246,7 +1193,6 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
 
   auto Du_Gs              = Du_spec * GammaB_sigma;
   auto adjDu_GH_Ds        = g5 * adj(Du_tf) * g5 * Gamma_H * Ds_ti;
-  //auto Gn_adjDd_GH_Du     = GammaB_nucl * g5 * adj(Dd_tf) * g5 * Gamma_H * Du_ti;
   auto adjDd_GH_Du        = g5 * adj(Dd_tf) * g5 * Gamma_H * Du_ti;
   auto Gn_adjDd_GH_Du     = GammaB_nucl * adjDd_GH_Du; // for some reason I needed to split this into two lines to avoid the compilation error 'error: identifier "Grid::Gamma::mul" is undefined in device code'
 
@@ -1259,43 +1205,45 @@ void BaryonUtils<FImpl>::SigmaToNucleonQ2NonEyeSite(const mobj &Du_ti,
     int b_n    = (ie_n < 3 ? (ie_n+1)%3 : (8-ie_n)%3 ); //epsilon[ie_n][1]; //b
     int c_n    = (ie_n < 3 ? (ie_n+2)%3 : (7-ie_n)%3 ); //epsilon[ie_n][2]; //c
     int eSgn_n = (ie_n < 3 ? 1 : -1);
-  for (int ie_s=0; ie_s < 6 ; ie_s++){
-    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_s=0; ie_s < 6 ; ie_s++){
+      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);
 
-    ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
+      ee = Real(eSgn_n * eSgn_s); //epsilon_sgn[ie_n] * epsilon_sgn[ie_s];
 
-    for (int alpha_n=0; alpha_n<Ns; alpha_n++){
-    for (int beta_s=0;  beta_s<Ns;   beta_s++){
+      for (int alpha_n=0; alpha_n<Ns; alpha_n++){
+      for (int beta_s=0;  beta_s<Ns;   beta_s++){
 
-      auto adjDu_GH_Ds_ab_ab = adjDu_GH_Ds()(alpha_n, beta_s)(a_n,b_s);
-      auto Gn_adjDd_GH_Du_Gs_ab_ba = Gn_adjDd_GH_Du_Gs()(alpha_n, beta_s)(b_n,a_s);
+        auto adjDu_GH_Ds_ab_ab = adjDu_GH_Ds()(alpha_n, beta_s)(a_n,b_s);
+        auto Gn_adjDd_GH_Du_Gs_ab_ba = Gn_adjDd_GH_Du_Gs()(alpha_n, beta_s)(b_n,a_s);
 
-      for (int gamma_s=0; gamma_s<Ns; gamma_s++){
-        auto Gn_adjDd_GH_Du_ag_bc = Gn_adjDd_GH_Du()(alpha_n,gamma_s)(b_n,c_s);
-      for (int gamma_n=0; gamma_n<Ns; gamma_n++){
-        auto adjDu_GH_Ds_gb_cb = adjDu_GH_Ds()(gamma_n, beta_s)(c_n,b_s);
+        for (int gamma_s=0; gamma_s<Ns; gamma_s++){
+          auto Gn_adjDd_GH_Du_ag_bc = Gn_adjDd_GH_Du()(alpha_n,gamma_s)(b_n,c_s);
+          for (int gamma_n=0; gamma_n<Ns; gamma_n++){
+            auto adjDu_GH_Ds_gb_cb = adjDu_GH_Ds()(gamma_n, beta_s)(c_n,b_s);
 
-        result()(gamma_n,gamma_s)() += ee * adjDu_GH_Ds_ab_ab
+            result()(gamma_n,gamma_s)() += ee * adjDu_GH_Ds_ab_ab
                                           * Gn_adjDd_GH_Du_Gs_ab_ba
                                           * Du_spec()(gamma_n,gamma_s)(c_n,c_s);
 
-        result()(gamma_n,gamma_s)() -= ee * adjDu_GH_Ds_gb_cb
+            result()(gamma_n,gamma_s)() -= ee * adjDu_GH_Ds_gb_cb
                                           * Gn_adjDd_GH_Du_Gs_ab_ba
                                           * Du_spec()(alpha_n,gamma_s)(a_n,c_s);
 
-        result()(gamma_n,gamma_s)() += ee * adjDu_GH_Ds_gb_cb
+            result()(gamma_n,gamma_s)() += ee * adjDu_GH_Ds_gb_cb
                                           * Gn_adjDd_GH_Du_ag_bc
                                           * Du_Gs()(alpha_n, beta_s)(a_n,a_s);
 
-        result()(gamma_n,gamma_s)() -= ee * adjDu_GH_Ds_ab_ab
+            result()(gamma_n,gamma_s)() -= ee * adjDu_GH_Ds_ab_ab
                                           * Gn_adjDd_GH_Du_ag_bc
                                           * Du_Gs()(gamma_n, beta_s)(c_n,a_s);
+          }
+        }
       }}
-    }}
-  }}
+    }
+  }
 }
 
 template<class FImpl>