diff --git a/lib/qcd/action/ActionBase.h b/lib/qcd/action/ActionBase.h
index 8b91fb79..7c9f6f1c 100644
--- a/lib/qcd/action/ActionBase.h
+++ b/lib/qcd/action/ActionBase.h
@@ -35,6 +35,7 @@ template<class GaugeField>
 class Action { 
 
  public:
+  bool is_smeared = false;
   // Boundary conditions? // Heatbath?
   virtual void  refresh(const GaugeField &U, GridParallelRNG& pRNG) = 0;// refresh pseudofermions
   virtual RealD S    (const GaugeField &U)                        = 0;  // evaluate the action
diff --git a/lib/qcd/hmc/HMC.h b/lib/qcd/hmc/HMC.h
index 0b87b03e..b471eb3a 100644
--- a/lib/qcd/hmc/HMC.h
+++ b/lib/qcd/hmc/HMC.h
@@ -160,7 +160,7 @@ namespace Grid{
       /////////////////////////////////////////
       // Constructor
       /////////////////////////////////////////
-      HybridMonteCarlo(HMCparameters Pms,  IntegratorType &_Int, GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, GaugeField &_U ) :
+      HybridMonteCarlo(HMCparameters Pms,  IntegratorType &_Int, GridSerialRNG &_sRNG, GridParallelRNG &_pRNG, GaugeField &_U) :
         Params(Pms), 
 	TheIntegrator(_Int), 
 	sRNG(_sRNG),
diff --git a/lib/qcd/hmc/HmcRunner.h b/lib/qcd/hmc/HmcRunner.h
index ce92b1d2..6dac0d53 100644
--- a/lib/qcd/hmc/HmcRunner.h
+++ b/lib/qcd/hmc/HmcRunner.h
@@ -81,11 +81,22 @@ public:
       NumTraj = ivec[0];
     }
 
-    // Create integrator
-    typedef MinimumNorm2<GaugeField>  IntegratorType;// change here to change the algorithm
-    IntegratorParameters MDpar(20);
-    IntegratorType MDynamics(UGrid,MDpar, TheAction);
 
+    GridSerialRNG    sRNG;
+    GridParallelRNG  pRNG(UGrid);
+    LatticeGaugeField  U(UGrid); // change this to an extended field (smearing class)
+
+    std::vector<int> SerSeed({1,2,3,4,5});
+    std::vector<int> ParSeed({6,7,8,9,10});
+
+    
+    // Create integrator, including the smearing policy
+    SmearedConfiguration<Gimpl> SmearingPolicy; // simplest empty smearer, construct here more complex smearers
+    typedef MinimumNorm2<GaugeField, SmearedConfiguration<Gimpl> >  IntegratorType;// change here to change the algorithm
+    IntegratorParameters MDpar(20);
+    IntegratorType MDynamics(UGrid, MDpar, TheAction, SmearingPolicy);
+
+    
     // Checkpoint strategy
     NerscHmcCheckpointer<Gimpl> Checkpoint(std::string("ckpoint_lat"),std::string("ckpoint_rng"),1);
     PlaquetteLogger<Gimpl>      PlaqLog(std::string("plaq"));
@@ -94,12 +105,6 @@ public:
     HMCpar.StartTrajectory = StartTraj;
     HMCpar.Trajectories    = NumTraj;
     
-    GridSerialRNG    sRNG;
-    GridParallelRNG  pRNG(UGrid);
-    LatticeGaugeField  U(UGrid); // change this to an extended field (smearing class)
-
-    std::vector<int> SerSeed({1,2,3,4,5});
-    std::vector<int> ParSeed({6,7,8,9,10});
 
     if ( StartType == HotStart ) {
       // Hot start
@@ -129,7 +134,11 @@ public:
       Checkpoint.CheckpointRestore(StartTraj, U, sRNG, pRNG);
     }
 
-    HybridMonteCarlo<GaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U); // pass the extended field
+    // Attach the gauge field to the smearing Policy and create the fill the smeared set
+    // notice that the unit configuration is singular in this procedure
+    SmearingPolicy.set_GaugeField(U);
+    
+    HybridMonteCarlo<GaugeField,IntegratorType>  HMC(HMCpar, MDynamics,sRNG,pRNG,U); 
     HMC.AddObservable(&Checkpoint);
     HMC.AddObservable(&PlaqLog);
     
diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h
index e6f6c1cd..5a70a1ec 100644
--- a/lib/qcd/hmc/integrators/Integrator.h
+++ b/lib/qcd/hmc/integrators/Integrator.h
@@ -68,7 +68,7 @@ namespace Grid{
     };
 
     /*! @brief Class for Molecular Dynamics management */   
-    template<class GaugeField>
+    template<class GaugeField, class SmearingPolicy>
     class Integrator {
 
     protected:
@@ -85,6 +85,8 @@ namespace Grid{
 
       GaugeField P;
 
+      SmearingPolicy &Smearer;
+      
       // Should match any legal (SU(n)) gauge field
       // Need to use this template to match Ncol to pass to SU<N> class
       template<int Ncol,class vec> void generate_momenta(Lattice< iVector< iScalar< iMatrix<vec,Ncol> >, Nd> > & P,GridParallelRNG& pRNG){
@@ -113,7 +115,9 @@ namespace Grid{
       void update_P(GaugeField &Mom,GaugeField&U, int level,double ep){
 	for(int a=0; a<as[level].actions.size(); ++a){
 	  GaugeField force(U._grid);
-	  as[level].actions.at(a)->deriv(U,force);
+	  as[level].actions.at(a)->deriv(U,force); // deriv should not include Ta
+	  if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force);
+	  force = Ta(force);
 	  Mom = Mom - force*ep;
 	}
       }
@@ -137,6 +141,8 @@ namespace Grid{
 	  ProjectOnGroup(Umu);
 	  PokeIndex<LorentzIndex>(U, Umu, mu);
 	}
+	// Update the smeared fields, can be implemented as observer
+	Smearer.set_GaugeField(U);
       }
       
       virtual void step (GaugeField& U,int level, int first,int last)=0;
@@ -145,14 +151,17 @@ namespace Grid{
 
       Integrator(GridBase* grid, 
 		 IntegratorParameters Par,
-		 ActionSet<GaugeField> & Aset):
+		 ActionSet<GaugeField> & Aset,
+		 SmearingPolicy &Sm):
           Params(Par),
     	  as(Aset),
 	  P(grid),
-	  levels(Aset.size())
+	  levels(Aset.size()),
+	  Smearer(Sm)
       {
 	t_P.resize(levels,0.0);
 	t_U=0.0;
+	// initialization of smearer delegated outside of Integrator
       };
       
       virtual ~Integrator(){}
@@ -163,7 +172,10 @@ namespace Grid{
 	generate_momenta(P,pRNG);
 	for(int level=0; level< as.size(); ++level){
 	  for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
-	    as[level].actions.at(actionID)->refresh(U, pRNG);
+	    // get gauge field from the SmearingPolicy and
+	    // based on the boolean is_smeared in actionID
+	    GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+	    as[level].actions.at(actionID)->refresh(Us, pRNG);
 	  }
 	}
       }
@@ -186,7 +198,10 @@ namespace Grid{
 	// Actions
 	for(int level=0; level<as.size(); ++level){
 	  for(int actionID=0; actionID<as[level].actions.size(); ++actionID){
-	    Hterm = as[level].actions.at(actionID)->S(U);
+	    // get gauge field from the SmearingPolicy and
+	    // based on the boolean is_smeared in actionID
+	    GaugeField& Us = Smearer.get_U(as[level].actions.at(actionID)->is_smeared);
+	    Hterm = as[level].actions.at(actionID)->S(Us);
 	    std::cout<<GridLogMessage << "Level "<<level<<" term "<<actionID<<" H = "<<Hterm<<std::endl;
 	    H += Hterm;
 	  }
diff --git a/lib/qcd/hmc/integrators/Integrator_algorithm.h b/lib/qcd/hmc/integrators/Integrator_algorithm.h
index eb1b30ad..66390742 100644
--- a/lib/qcd/hmc/integrators/Integrator_algorithm.h
+++ b/lib/qcd/hmc/integrators/Integrator_algorithm.h
@@ -91,14 +91,17 @@ namespace Grid{
     *  P 1/2                            P 1/2
     */    
 
-    template<class GaugeField> class LeapFrog : public Integrator<GaugeField> {
+    template<class GaugeField, class SmearingPolicy> class LeapFrog :
+      public Integrator<GaugeField, SmearingPolicy> {
     public:
 
-      typedef LeapFrog<GaugeField> Algorithm;
+      typedef LeapFrog<GaugeField, SmearingPolicy> Algorithm;
 
       LeapFrog(GridBase* grid, 
 	       IntegratorParameters Par,
-	       ActionSet<GaugeField> & Aset): Integrator<GaugeField>(grid,Par,Aset) {};
+	       ActionSet<GaugeField> & Aset,
+	       SmearingPolicy & Sm):
+	Integrator<GaugeField, SmearingPolicy>(grid,Par,Aset,Sm) {};
 
 
       void step (GaugeField& U, int level,int _first, int _last){
@@ -135,7 +138,8 @@ namespace Grid{
       }
     };
 
-    template<class GaugeField> class MinimumNorm2 : public Integrator<GaugeField> {
+    template<class GaugeField, class SmearingPolicy> class MinimumNorm2 :
+      public Integrator<GaugeField, SmearingPolicy> {
     private:
       const RealD lambda = 0.1931833275037836;
 
@@ -143,7 +147,9 @@ namespace Grid{
 
       MinimumNorm2(GridBase* grid, 
 		   IntegratorParameters Par,
-		   ActionSet<GaugeField> & Aset): Integrator<GaugeField>(grid,Par,Aset) {};
+		   ActionSet<GaugeField> & Aset,
+		   SmearingPolicy& Sm):
+	Integrator<GaugeField, SmearingPolicy>(grid,Par,Aset,Sm) {};
 
       void step (GaugeField& U, int level, int _first,int _last){
 
@@ -191,7 +197,8 @@ namespace Grid{
     };
 
 
-    template<class GaugeField> class ForceGradient : public Integrator<GaugeField> {
+    template<class GaugeField, class SmearingPolicy> class ForceGradient :
+      public Integrator<GaugeField, SmearingPolicy> {
     private:
       const RealD lambda = 1.0/6.0;;
       const RealD chi    = 1.0/72.0;
@@ -202,7 +209,9 @@ namespace Grid{
       // Looks like dH scales as dt^4. tested wilson/wilson 2 level.
     ForceGradient(GridBase* grid, 
 		  IntegratorParameters Par,
-		  ActionSet<GaugeField> & Aset): Integrator<GaugeField>(grid,Par,Aset) {};
+		  ActionSet<GaugeField> & Aset,
+		  SmearingPolicy &Sm):
+      Integrator<GaugeField, SmearingPolicy>(grid,Par,Aset, Sm) {};
 
 
       void FG_update_P(GaugeField&U, int level,double fg_dt,double ep){
diff --git a/lib/qcd/smearing/APEsmearing.h b/lib/qcd/smearing/APEsmearing.h
index 2580d296..0636fc2c 100644
--- a/lib/qcd/smearing/APEsmearing.h
+++ b/lib/qcd/smearing/APEsmearing.h
@@ -35,21 +35,22 @@ namespace Grid {
       Smear_APE():rho(set_rho(1.0)){}
       ~Smear_APE(){}
 
-
       void smear(GaugeField& u_smr, const GaugeField& U)const{
 	GridBase *grid = U._grid;
 	double d_rho;
 	GaugeLinkField Cup(grid), tmp_stpl(grid);
 	WilsonLoops<Gimpl> WL;
-	u_smr = zero; // probably unecessary
+	u_smr = zero; 
 
 	for(int mu=0; mu<Nd; ++mu){
 	  Cup = zero;
 	  for(int nu=0; nu<Nd; ++nu){
 	    d_rho = rho[mu + Nd * nu];
+	    // get the staple in direction mu, nu
 	    WL.Staple(tmp_stpl, U, mu, nu);  //nb staple conventions of IroIro and Grid differ by a dag
 	    Cup += tmp_stpl*d_rho;
 	  }
+	  // save the Cup link-field on the u_smr gauge-field
 	  pokeLorentz(u_smr, adj(Cup), mu); // u_smr[mu] = Cup^dag
 	}
       }
@@ -66,8 +67,10 @@ namespace Grid {
 	int vol = U._grid->gSites();
 	
 	WilsonLoops<Gimpl> WL;
-	GaugeLinkField staple(grid), u_tmp(grid), iLambda_mu(grid), iLambda_nu(grid);
-	GaugeLinkField U_mu(grid), U_nu(grid), sh_field(grid), temp_Sigma(grid);
+	GaugeLinkField staple(grid), u_tmp(grid);
+	GaugeLinkField iLambda_mu(grid), iLambda_nu(grid);
+	GaugeLinkField U_mu(grid), U_nu(grid);
+	GaugeLinkField sh_field(grid), temp_Sigma(grid);
 	Real rho_munu, rho_numu;
     
 	for(int mu = 0; mu < Nd; ++mu){
diff --git a/lib/qcd/smearing/BaseSmearing.h b/lib/qcd/smearing/BaseSmearing.h
index 377bcbc9..40e82faa 100644
--- a/lib/qcd/smearing/BaseSmearing.h
+++ b/lib/qcd/smearing/BaseSmearing.h
@@ -7,7 +7,7 @@
 template <class Gimpl> 
 class Smear{
 public:
-  INHERIT_GIMPL_TYPES(Gimpl)
+  INHERIT_GIMPL_TYPES(Gimpl) // inherits the types for the gauge fields
 
   virtual ~Smear(){}
   virtual void smear     (GaugeField&,const GaugeField&)const = 0;
diff --git a/lib/qcd/smearing/GaugeConfiguration.h b/lib/qcd/smearing/GaugeConfiguration.h
index 42baced8..efb5c34d 100644
--- a/lib/qcd/smearing/GaugeConfiguration.h
+++ b/lib/qcd/smearing/GaugeConfiguration.h
@@ -22,7 +22,7 @@ namespace Grid {
       It stores a list of smeared configurations.
     */
     template <class Gimpl>
-    class SmearedConfiguration {
+      class SmearedConfiguration {
     public:
       INHERIT_GIMPL_TYPES(Gimpl) ;
       
@@ -33,17 +33,23 @@ namespace Grid {
       
       // Member functions
       //====================================================================
-      void fill_smearedSet(){
-	GaugeField previous_u;
+      void fill_smearedSet(GaugeField& U){
+	ThinLinks = &U; //attach the smearing routine to the field U
+	//check that the pointer is not null
+	if (ThinLinks==NULL) 
+	  std::cout << GridLogError << "[SmearedConfiguration] Error in ThinLinks pointer\n";
 	
-	std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n";
-	
-	previous_u = *ThinLinks;
-	for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){
-	  StoutSmearing.smear(SmearedSet[smearLvl],previous_u);
-	  previous_u = SmearedSet[smearLvl];
+	if (smearingLevels > 0){
+	  GaugeField previous_u(U._grid);
+	  std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n";
+	  
+	  previous_u = *ThinLinks;
+	  for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){
+	    StoutSmearing.smear(SmearedSet[smearLvl],previous_u);
+	    previous_u = SmearedSet[smearLvl];
+	  }
+
 	}
-	
       }
       //====================================================================
       GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, 
@@ -86,9 +92,9 @@ namespace Grid {
 	
 	LatticeReal u(grid), w(grid);
 	LatticeComplex f0(grid), f1(grid), f2(grid);
-	LatticeReal xi0(grid), xi1(grid), fden(grid), tmp(grid);
+	LatticeReal xi0(grid), xi1(grid), tmp(grid);
 	LatticeReal u2(grid), w2(grid), cosw(grid);
-	LatticeComplex emiu(grid), e2iu(grid), qt(grid);
+	LatticeComplex emiu(grid), e2iu(grid), qt(grid), fden(grid);
 	LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid);
 	LatticeComplex r22(grid), tr1(grid), tr2(grid);
 	LatticeComplex b10(grid), b11(grid), b12(grid), b20(grid), b21(grid), b22(grid);
@@ -173,66 +179,70 @@ namespace Grid {
 				    links configuration */
       
       /*! @brief Standard constructor */
-      SmearedConfiguration(GridCartesian * UGrid,
-			   unsigned int Nsmear, 
-			   Smear_Stout<Gimpl>& Stout):
-	smearingLevels(Nsmear),
+    SmearedConfiguration(GridCartesian * UGrid,
+			 unsigned int Nsmear, 
+			 Smear_Stout<Gimpl>& Stout):
+      smearingLevels(Nsmear),
 	StoutSmearing(Stout),
-	ThinLinks(new GaugeField(UGrid)){
+	ThinLinks(NULL){
 	for (unsigned int i=0; i< smearingLevels; ++i)
 	  SmearedSet.push_back(*(new GaugeField(UGrid)));
       }
       
       /*! For just thin links */
-      SmearedConfiguration(GridCartesian * UGrid):
-	smearingLevels(0),
+    SmearedConfiguration():
+      smearingLevels(0),
 	StoutSmearing(),
-	SmearedSet(0),
-	ThinLinks(new GaugeField(UGrid)){}
+	SmearedSet(),
+	ThinLinks(NULL){}
+      
+      
+      // attach the smeared routines to the thin links U and fill the smeared set
+      void set_GaugeField(GaugeField& U){ fill_smearedSet(U);}
       
-      void set_GaugeField(){ fill_smearedSet(); }
       void smeared_force(GaugeField& SigmaTilde) const{
 	GaugeField force = SigmaTilde;//actually = U*SigmaTilde, check this for Grid
 	GaugeLinkField tmp_mu(SigmaTilde._grid);
 	
 	for (int mu = 0; mu < Nd; mu++){
 	  tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu);
-	  pokeLorentz(force, tmp_mu, mu);
+	    pokeLorentz(force, tmp_mu, mu);
 	}
 	for(int ismr = smearingLevels - 1; ismr > 0; --ismr)
 	  force = AnalyticSmearedForce(force,get_smeared_conf(ismr-1));
-
+	
 	force = AnalyticSmearedForce(force,*ThinLinks);
-
+	
 	for (int mu = 0; mu < Nd; mu++){
 	  tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu);
 	  pokeLorentz(SigmaTilde, tmp_mu, mu);
 	}
       }
-
-
-      GaugeField* get_SmearedU() const{ 
-	return const_cast<GaugeField*>(&(SmearedSet[smearingLevels-1]));
+      
+      
+      GaugeField& get_SmearedU(){ 
+	return SmearedSet[smearingLevels-1];
       }
-
-      GaugeField* get_U(bool smeared=false) const { 
+      
+      GaugeField& get_U(bool smeared=false) { 
 	// get the config, thin links by default
 	if (smeared){
 	  if (smearingLevels) return get_SmearedU();
-	  else                return ThinLinks;
+	  else                return *ThinLinks;
 	}
-	else return ThinLinks;
+	else return *ThinLinks;
       }
       
     };
     
     
   }
-
+  
 }
 
 
 
 
 
+
 #endif
diff --git a/lib/qcd/smearing/StoutSmearing.h b/lib/qcd/smearing/StoutSmearing.h
index 3538417a..dc4b29be 100644
--- a/lib/qcd/smearing/StoutSmearing.h
+++ b/lib/qcd/smearing/StoutSmearing.h
@@ -20,10 +20,14 @@ namespace Grid {
     public:
       INHERIT_GIMPL_TYPES(Gimpl)
       
-      Smear_Stout(Smear < Gimpl >* base):SmearBase(base){}
+      Smear_Stout(Smear < Gimpl >* base):SmearBase(base){
+	static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3");
+      }
       
       /*! Default constructor */
-      Smear_Stout():SmearBase(new Smear_APE < Gimpl > ()){}
+      Smear_Stout():SmearBase(new Smear_APE < Gimpl > ()){
+	static_assert(Nc==3, "Stout smearing currently implemented only for Nc==3");
+      }
       
       ~Smear_Stout(){}
       
diff --git a/tests/Test_main.cc b/tests/Test_main.cc
index 467c8a5e..9236fae1 100644
--- a/tests/Test_main.cc
+++ b/tests/Test_main.cc
@@ -608,7 +608,7 @@ int main (int argc, char ** argv)
 
 
 
- // Testing Smearing routine compilation
+ // Testing Smearing routine compilation, separate in a different file
  GridCartesian           Fine(latt_size,simd_layout,mpi_layout);
  Smear_APE< PeriodicGimplR > APEsmearing; // periodic gauge implemetation
  Smear_Stout< PeriodicGimplR > StoutSmearing(&APEsmearing);