diff --git a/lib/qcd/action/gauge/PlaqPlusRectangleAction.h b/lib/qcd/action/gauge/PlaqPlusRectangleAction.h
index 82e21aed..0f39b8f3 100644
--- a/lib/qcd/action/gauge/PlaqPlusRectangleAction.h
+++ b/lib/qcd/action/gauge/PlaqPlusRectangleAction.h
@@ -34,29 +34,43 @@ namespace Grid{
 	return action;
       };
 
-      virtual void deriv(const GaugeField &U,GaugeField & dSdU) {
+      virtual void deriv(const GaugeField &Umu,GaugeField & dSdU) {
 	//extend Ta to include Lorentz indexes
 	RealD factor_p = c_plaq/RealD(Nc)*0.5;
 	RealD factor_r =   c_rect/RealD(Nc)*0.5;
 
-	GaugeLinkField Umu(U._grid);
-	GaugeLinkField dSdU_mu(U._grid);
-	GaugeLinkField staple(U._grid);
+	GridBase *grid = Umu._grid;
+
+	std::vector<GaugeLinkField> U (Nd,grid);
+	std::vector<GaugeLinkField> U2(Nd,grid);
+
+	for(int mu=0;mu<Nd;mu++){
+	  U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+	  WilsonLoops<GaugeField>::RectStapleDouble(U2[mu],U[mu],mu);
+	}
+
+	GaugeLinkField dSdU_mu(grid);
+	GaugeLinkField staple(grid);
+
 	for (int mu=0; mu < Nd; mu++){
 
-	  Umu = PeekIndex<LorentzIndex>(U,mu);
-
 	  // Staple in direction mu
-	  WilsonLoops<GaugeField>::Staple(staple,U,mu);
-	  dSdU_mu = Ta(Umu*staple)*factor_p;
-	  
-	  WilsonLoops<GaugeField>::RectStaple(staple,U,mu);
-	  dSdU_mu = dSdU_mu + Ta(Umu*staple)*factor_r;
+
+	  WilsonLoops<GaugeField>::Staple(staple,Umu,mu);
+
+	  dSdU_mu = Ta(U[mu]*staple)*factor_p;
+
+	  //	  WilsonLoops<GaugeField>::RectStaple(staple,Umu,mu);
+
+	  WilsonLoops<GaugeField>::RectStapleOptimised(staple,U2,U,mu);
+
+	  dSdU_mu = dSdU_mu + Ta(U[mu]*staple)*factor_r;
 	  
 	  PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
-
 	}
+
       };
+
     };
 
     // Convenience for common physically defined cases.
diff --git a/lib/qcd/utils/CovariantCshift.h b/lib/qcd/utils/CovariantCshift.h
index e090f497..8bc0986f 100644
--- a/lib/qcd/utils/CovariantCshift.h
+++ b/lib/qcd/utils/CovariantCshift.h
@@ -5,6 +5,7 @@ namespace QCD {
 ////////////////////////////////////////////////////////////////////////
 // Low performance implementation of CovariantCshift API
 ////////////////////////////////////////////////////////////////////////
+// Make these members of an Impl class for BC's.
 template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link, 
 									    int mu,
 									    const Lattice<covariant> &field)
@@ -19,5 +20,6 @@ template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const
   tmp = adj(Link)*field;
   return Cshift(tmp,mu,-1);// moves towards positive mu
 }
+
 }}
 #endif
diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h
index cc6eab0b..36fbe51c 100644
--- a/lib/qcd/utils/WilsonLoops.h
+++ b/lib/qcd/utils/WilsonLoops.h
@@ -206,8 +206,96 @@ public:
   //////////////////////////////////////////////////
   // the sum over all staples on each site
   //////////////////////////////////////////////////
+  static void RectStapleDouble(GaugeMat &U2,const GaugeMat & U,int mu){
+    U2 = U * Cshift(U,mu,1);
+  }
+  static void RectStapleOptimised(GaugeMat &Stap,std::vector<GaugeMat> &U2,std::vector<GaugeMat> &U,int mu){
+
+    Stap = zero;
+
+    GridBase *grid = U[0]._grid;
+
+    GaugeMat Staple2x1 (grid);
+    GaugeMat tmp (grid);
+
+    for(int nu=0;nu<Nd;nu++){
+      if ( nu!=mu) {
+
+	// Up staple    ___ ___ 
+	//             |       |
+	tmp = Cshift(adj(U[nu]),nu,-1); 
+	tmp = adj(U2[mu])*tmp;
+	tmp = Cshift(tmp,mu,-2);
+
+	Staple2x1 = CovShiftForward (U[nu],nu,tmp);
+
+
+	// Down staple
+	//             |___ ___|
+	//
+	tmp = adj(U2[mu])*U[nu];
+	Staple2x1+= CovShiftBackward(U[nu],nu,Cshift(tmp,mu,-2));
+
+
+	//              ___ ___
+	//             |    ___|
+	//             |___ ___|
+	//
+
+	Stap+= Cshift(CovShiftForward (U[mu],mu,Staple2x1),mu,1);
+
+	//              ___ ___
+	//             |___    |
+	//             |___ ___|
+	//
+
+	//	tmp= Staple2x1* Cshift(U[mu],mu,-2);
+	//	Stap+= Cshift(tmp,mu,1) ;
+	Stap+= Cshift(Staple2x1,mu,1)*Cshift(U[mu],mu,-1); ;
+
+	//       --    
+	//      |  |              
+	//          
+	//      |  | 
+	
+	tmp = Cshift(adj(U2[nu]),nu,-2);
+	tmp = CovShiftBackward(U[mu],mu,tmp);
+	tmp = U2[nu]*Cshift(tmp,nu,2);
+	Stap+= Cshift(tmp, mu, 1);
+
+	//      |  |              
+	//          
+	//      |  | 
+	//       -- 
+	
+	tmp = CovShiftBackward(U[mu],mu,U2[nu]);
+	tmp = adj(U2[nu])*tmp;
+	tmp = Cshift(tmp,nu,-2);
+	Stap+=Cshift(tmp, mu, 1);
+    }}
+
+
+  }
   static void RectStaple(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
 
+#if 1
+    RectStapleUnoptimised(Stap,Umu,mu);
+#else
+    GridBase *grid = Umu._grid;
+
+    std::vector<GaugeMat> U (Nd,grid);
+    std::vector<GaugeMat> U2(Nd,grid);
+
+    for(int mu=0;mu<Nd;mu++){
+      U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
+      RectStapleDouble(U2[mu],U[mu],mu);
+    }
+    
+    RectStapleOptimised(Stap,U2,U,mu);
+#endif
+  }
+
+  static void RectStapleUnoptimised(GaugeMat &Stap,const GaugeLorentz &Umu,int mu){
     GridBase *grid = Umu._grid;
 
     std::vector<GaugeMat> U(4,grid);