diff --git a/lib/qcd/action/gauge/GaugeActionBase.h b/lib/qcd/action/gauge/GaugeActionBase.h
new file mode 100644
index 00000000..55141965
--- /dev/null
+++ b/lib/qcd/action/gauge/GaugeActionBase.h
@@ -0,0 +1,20 @@
+#ifndef QCD_GAUGE_ACTION_BASE
+#define QCD_GAUGE_ACTION_BASE
+namespace Grid {
+namespace QCD{
+
+template<class GaugeField>
+class GaugeActionBase { // derive this from TermInAction?
+
+ public:
+  virtual RealD S(const GaugeField &U)                        = 0;  // evaluate the action
+  virtual void  deriv(const GaugeField &U,GaugeField & dSdU ) = 0;  // evaluate the action derivative
+  virtual void  staple(const GaugeField &U,GaugeField & dSdU ) = 0;  // evaluate the action derivative
+  virtual void refresh(const GaugeField & ) {}; 
+  // Boundary conditions?
+  // Heatbath?
+  virtual ~GaugeActionBase() {};
+};
+
+}}
+#endif
diff --git a/lib/qcd/action/gauge/WilsonGaugeAction.h b/lib/qcd/action/gauge/WilsonGaugeAction.h
new file mode 100644
index 00000000..8ecb1f81
--- /dev/null
+++ b/lib/qcd/action/gauge/WilsonGaugeAction.h
@@ -0,0 +1,24 @@
+#ifndef QCD_WILSON_GAUGE_ACTION_H
+#define QCD_WILSON_GAUGE_ACTION_H
+
+////////////////////////////////////////////////////////////////////////
+// Wilson Gauge Action .. should I template the Nc etc..
+////////////////////////////////////////////////////////////////////////
+template<class GaugeField,class MatrixField>
+class WilsonGaugeAction : public GaugeActionBase<GaugeField> {
+ public:
+  virtual RealD S(const GaugeField &U) {
+    return WilsonLoops<MatrixField,GaugeField>::sumPlaquette(U);
+  };
+  virtual RealD deriv(GaugeField &U,GaugeField & dSdU ) {
+    WilsonLoops<MatrixField,GaugeField>::Staple(dSdU,U,mu);
+  };
+  virtual void  staple(const MatrixField &stap,GaugeField & U,int mu ) {
+    WilsonLoops<MatrixField,GaugeField>::Staple(stap,U,mu);
+  };
+};
+
+
+#endif
+
+#endif
diff --git a/lib/qcd/utils/CovariantCshift.h b/lib/qcd/utils/CovariantCshift.h
new file mode 100644
index 00000000..0c621069
--- /dev/null
+++ b/lib/qcd/utils/CovariantCshift.h
@@ -0,0 +1,22 @@
+#ifndef QCD_UTILS_COVARIANT_CSHIFT_H
+#define QCD_UTILS_COVARIANT_CSHIFT_H
+namespace Grid {
+namespace QCD {
+////////////////////////////////////////////////////////////////////////
+// Low performance implementation of CovariantCshift API
+////////////////////////////////////////////////////////////////////////
+template<class covariant,class gauge> Lattice<covariant> CovShiftForward(const Lattice<gauge> &Link, 
+									    int mu,
+									    const Lattice<covariant> &field)
+{
+  return Link*Cshift(field,mu,1);// moves towards negative mu
+}
+template<class covariant,class gauge> Lattice<covariant> CovShiftBackward(const Lattice<gauge> &Link, 
+									     int mu,
+									     const Lattice<covariant> &field)
+{
+  Lattice<covariant> tmp(field._grid);
+  tmp = adj(Link)*field;
+  return Cshift(tmp,mu,-1);// moves towards positive mu
+}
+}}
diff --git a/lib/qcd/utils/WilsonLoops.h b/lib/qcd/utils/WilsonLoops.h
new file mode 100644
index 00000000..f2d6d4de
--- /dev/null
+++ b/lib/qcd/utils/WilsonLoops.h
@@ -0,0 +1,130 @@
+#ifndef QCD_UTILS_WILSON_LOOPS_H
+#define QCD_UTILS_WILSON_LOOPS_H
+namespace Grid {
+namespace QCD {
+
+// Common wilson loop observables
+template<class GaugeMat,class GaugeLorentz>
+class WilsonLoops {
+public:
+
+  //////////////////////////////////////////////////
+  // directed plaquette oriented in mu,nu plane
+  //////////////////////////////////////////////////
+  static void dirPlaquette(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu)
+  {
+    plaq=CovShiftForward(U[mu],mu,U[nu])*adj(CovShiftForward(U[nu],nu,U[mu]));
+  }
+  //////////////////////////////////////////////////
+  // trace of directed plaquette oriented in mu,nu plane
+  //////////////////////////////////////////////////
+  static void traceDirPlaquette(LatticeComplex &plaq, const std::vector<GaugeMat> &U, const int mu, const int nu)
+  {
+    GaugeMat sp(U[0]._grid);
+    dirPlaquette(sp,U,mu,nu);
+    plaq=trace(sp);
+  }
+  //////////////////////////////////////////////////
+  // sum over all planes of plaquette
+  //////////////////////////////////////////////////
+  static void sitePlaquette(LatticeComplex &Plaq,const std::vector<GaugeMat> &U)
+  {
+    LatticeComplex sitePlaq(U[0]._grid);
+    for(int mu=1;mu<Nd;mu++){
+      for(int nu=0;nu<mu;nu++){
+	traceDirPlaquette(sitePlaq,U,mu,nu);
+	Plaq = Plaq + sitePlaq;
+      }
+    }
+  }
+  //////////////////////////////////////////////////
+  // sum over all x,y,z,t and over all planes of plaquette
+  //////////////////////////////////////////////////
+  static RealD sumPlaquette(const GaugeLorentz &Umu){
+    std::vector<GaugeMat> U(4,Umu._grid);
+    for(int mu=0;mu<Nd;mu++){
+      U[mu] = peekIndex<LorentzIndex>(Umu,mu);
+    }
+
+    LatticeComplex Plaq(Umu._grid);
+    
+    sitePlaquette(Plaq,U);
+    
+    TComplex Tp = sum(Plaq);
+    Complex p  = TensorRemove(Tp);
+    return p.real();
+  }
+  //////////////////////////////////////////////////
+  // average over all x,y,z,t and over all planes of plaquette
+  //////////////////////////////////////////////////
+  static RealD avgPlaquette(const GaugeLorentz &Umu){
+
+    RealD sumplaq = sumPlaquette(Umu);
+    
+    double vol = Umu._grid->gSites();
+    
+    double faces = (1.0*Nd*(Nd-1))/2.0;
+    
+    return sumplaq/vol/faces/Nc; // Nd , Nc dependent... FIXME
+  }
+  //////////////////////////////////////////////////
+  // the sum over all staples on each site
+  //////////////////////////////////////////////////
+  static void Staple(GaugeMat &staple,GaugeLorentz &Umu,int mu){
+
+    std::vector<GaugeMat> U(4,Umu._grid);
+    for(int d=0;d<Nd;d++){
+      U[d] = peekIndex<LorentzIndex>(Umu,d);
+    }
+
+    staple = zero;
+    GaugeMat tmp(U[0]._grid);
+
+    for(int nu=0;nu<Nd;nu++){
+
+      if(nu != mu) {
+
+      // mu
+      // ^
+      // |__  nu
+
+      //    __                                         __
+      //      |         |                          
+      //    __|  =    __|                         *
+      //
+
+	staple   += CovShiftForward(U[nu],nu,U[mu])*Cshift(adj(U[nu]),mu,+1);
+
+      //
+      //  __        __                         
+      // |         |                          
+      // |__     = |                              *       __
+      //
+      //
+	tmp    = CovShiftForward (U[mu],mu,U[nu]);
+	staple+= CovShiftBackward(U[nu],nu,tmp);
+      
+      }
+    }
+  }
+
+  //////////////////////////////////////////////////////
+  // Similar to above for rectangle is required
+  //////////////////////////////////////////////////////
+  /*
+void siteRectangle(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){
+RealD avgRectangle(const std::vector<GaugeMat> &U){}
+RealD avgRectangle(const std::vector<GaugeMat> &U, const int mu, const int nu){}
+void traceRectangle(LatticeComplex &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){}
+void siteRectangle(GaugeMat &plaq,const std::vector<GaugeMat> &U, const int mu, const int nu){}
+  */
+
+};
+}}
+
+
+
+
+
+
+#endif