diff --git a/Grid/qcd/action/gauge/WilsonGaugeAction.h b/Grid/qcd/action/gauge/WilsonGaugeAction.h
index f535b54f..661b0ca5 100644
--- a/Grid/qcd/action/gauge/WilsonGaugeAction.h
+++ b/Grid/qcd/action/gauge/WilsonGaugeAction.h
@@ -69,15 +69,22 @@ public:
 
     RealD factor = 0.5 * beta / RealD(Nc);
 
-    GaugeLinkField Umu(U.Grid());
+    // GaugeLinkField Umu(U.Grid());
+    std::vector<GaugeLinkField> Uvec(Nd, U.Grid());
+    for (int d = 0; d < Nd; d++) {
+      Uvec[d] = PeekIndex<LorentzIndex>(U, d);
+    }
+
     GaugeLinkField dSdU_mu(U.Grid());
     for (int mu = 0; mu < Nd; mu++) {
 
-      Umu = PeekIndex<LorentzIndex>(U, mu);
-      
+      // Umu = PeekIndex<LorentzIndex>(U, mu);
+
       // Staple in direction mu
-      WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
-      dSdU_mu = Ta(Umu * dSdU_mu) * factor;
+      // WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
+      // dSdU_mu = Ta(Umu * dSdU_mu) * factor;
+      WilsonLoops<Gimpl>::Staple(dSdU_mu, Uvec, mu);
+      dSdU_mu = Ta(Uvec[mu] * dSdU_mu) * factor;
       
       PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
     }
diff --git a/HMC/FTHMC2p1f.cc b/HMC/FTHMC2p1f.cc
index 7d93d168..1e914e87 100644
--- a/HMC/FTHMC2p1f.cc
+++ b/HMC/FTHMC2p1f.cc
@@ -25,13 +25,20 @@ directory
 *************************************************************************************/
 /*  END LEGAL */
 #include <Grid/Grid.h>
+
+#if Nc == 3
 #include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
 #include <Grid/qcd/smearing/JacobianAction.h>
+#endif
 
 using namespace Grid;
 
 int main(int argc, char **argv)
 {
+#if Nc != 3
+#warning FTHMC2p1f will not work for Nc != 3
+  std::cout << "This program will currently only work for Nc == 3." << std::endl;
+#else
   std::cout << std::setprecision(12);
   
   Grid_init(&argc, &argv);
@@ -220,7 +227,6 @@ int main(int argc, char **argv)
   TheHMC.Run(SmearingPolicy); // for smearing
 
   Grid_finalize();
+#endif
 } // main
 
-
-
diff --git a/HMC/FTHMC2p1f_3GeV.cc b/HMC/FTHMC2p1f_3GeV.cc
index a8aa67f8..36d5caa3 100644
--- a/HMC/FTHMC2p1f_3GeV.cc
+++ b/HMC/FTHMC2p1f_3GeV.cc
@@ -24,14 +24,22 @@ See the full license in the file "LICENSE" in the top level distribution
 directory
 *************************************************************************************/
 /*  END LEGAL */
+
 #include <Grid/Grid.h>
+
+#if Nc == 3
 #include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
 #include <Grid/qcd/smearing/JacobianAction.h>
+#endif
 
 using namespace Grid;
 
 int main(int argc, char **argv)
 {
+#if Nc != 3
+#warning FTHMC2p1f_3GeV will not work for Nc != 3
+  std::cout << "This program will currently only work for Nc == 3." << std::endl;
+#else
   std::cout << std::setprecision(12);
   
   Grid_init(&argc, &argv);
@@ -220,6 +228,7 @@ int main(int argc, char **argv)
   TheHMC.Run(SmearingPolicy); // for smearing
 
   Grid_finalize();
+#endif
 } // main
 
 
diff --git a/HMC/HMC2p1f_3GeV.cc b/HMC/HMC2p1f_3GeV.cc
index 4bf088d7..199d4be8 100644
--- a/HMC/HMC2p1f_3GeV.cc
+++ b/HMC/HMC2p1f_3GeV.cc
@@ -25,13 +25,20 @@ directory
 *************************************************************************************/
 /*  END LEGAL */
 #include <Grid/Grid.h>
+
+#if Nc == 3
 #include <Grid/qcd/smearing/GaugeConfigurationMasked.h>
 #include <Grid/qcd/smearing/JacobianAction.h>
+#endif
 
 using namespace Grid;
 
 int main(int argc, char **argv)
 {
+#if Nc != 3
+#warning HMC2p1f_3GeV will not work for Nc != 3
+  std::cout << "This program will currently only work for Nc == 3." << std::endl;
+#else
   std::cout << std::setprecision(12);
   
   Grid_init(&argc, &argv);
@@ -220,6 +227,7 @@ int main(int argc, char **argv)
   TheHMC.Run(SmearingPolicy); // for smearing
 
   Grid_finalize();
+#endif
 } // main