diff --git a/Hadrons/Environment.hpp b/Hadrons/Environment.hpp
index 8498cf81..9841d665 100644
--- a/Hadrons/Environment.hpp
+++ b/Hadrons/Environment.hpp
@@ -94,20 +94,20 @@ public:
     void                    createGrid(const unsigned int Ls);
     template <typename VType = vComplex>
     void                    createCoarseGrid(const std::vector<int> &blockSize,
-                                             const unsigned int Ls = 1);
+                                             const unsigned int Ls);
     template <typename VType = vComplex>
-    GridCartesian *         getGrid(void) const;
+    GridCartesian *         getGrid(void);
     template <typename VType = vComplex>
-    GridRedBlackCartesian * getRbGrid(void) const;
+    GridRedBlackCartesian * getRbGrid(void);
     template <typename VType = vComplex>
-    GridCartesian *         getCoarseGrid(const std::vector<int> &blockSize) const;
+    GridCartesian *         getCoarseGrid(const std::vector<int> &blockSize);
     template <typename VType = vComplex>
-    GridCartesian *         getGrid(const unsigned int Ls) const;
+    GridCartesian *         getGrid(const unsigned int Ls);
     template <typename VType = vComplex>
-    GridRedBlackCartesian * getRbGrid(const unsigned int Ls) const;
+    GridRedBlackCartesian * getRbGrid(const unsigned int Ls);
     template <typename VType = vComplex>
     GridCartesian *         getCoarseGrid(const std::vector<int> &blockSize,
-                                          const unsigned int Ls) const;
+                                          const unsigned int Ls);
     std::vector<int>        getDim(void) const;
     int                     getDim(const unsigned int mu) const;
     unsigned int            getNd(void) const;
@@ -220,6 +220,14 @@ void Holder<T>::reset(T *pt)
  *                     Environment template implementation                    *
  ******************************************************************************/
 // grids ///////////////////////////////////////////////////////////////////////
+#define HADRONS_DUMP_GRID(...)\
+LOG(Debug) << "New grid " << (__VA_ARGS__) << std::endl;\
+LOG(Debug) << " - cb  : " << (__VA_ARGS__)->_isCheckerBoarded << std::endl;\
+LOG(Debug) << " - fdim: " << (__VA_ARGS__)->_fdimensions << std::endl;\
+LOG(Debug) << " - gdim: " << (__VA_ARGS__)->_gdimensions << std::endl;\
+LOG(Debug) << " - ldim: " << (__VA_ARGS__)->_ldimensions << std::endl;\
+LOG(Debug) << " - rdim: " << (__VA_ARGS__)->_rdimensions << std::endl;
+
 template <typename VType>
 void Environment::createGrid(const unsigned int Ls)
 {
@@ -231,15 +239,19 @@ void Environment::createGrid(const unsigned int Ls)
             SpaceTimeGrid::makeFourDimGrid(getDim(), 
                                         GridDefaultSimd(getNd(), VType::Nsimd()),
                                         GridDefaultMpi()));
+        HADRONS_DUMP_GRID(grid4d_[{hash, 1}].get());
         gridRb4d_[{hash, 1}].reset(
             SpaceTimeGrid::makeFourDimRedBlackGrid(grid4d_[{hash, 1}].get()));
+        HADRONS_DUMP_GRID(gridRb4d_[{hash, 1}].get());
     }
     if (grid5d_.find({hash, Ls}) == grid5d_.end())
     {
         auto g = grid4d_[{hash, 1}].get();
         
         grid5d_[{hash, Ls}].reset(SpaceTimeGrid::makeFiveDimGrid(Ls, g));
+        HADRONS_DUMP_GRID(grid5d_[{hash, Ls}].get());
         gridRb5d_[{hash, Ls}].reset(SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, g));
+        HADRONS_DUMP_GRID(gridRb5d_[{hash, Ls}].get());
     }
 }
 
@@ -289,18 +301,24 @@ void Environment::createCoarseGrid(const std::vector<int> &blockSize,
         gridCoarse4d_[hkey4d].reset(
             SpaceTimeGrid::makeFourDimGrid(coarseDim, 
                 GridDefaultSimd(nd, VType::Nsimd()), GridDefaultMpi()));
+        HADRONS_DUMP_GRID(gridCoarse4d_[hkey4d].get());
     }
     if (gridCoarse5d_.find(hkey5d) == gridCoarse5d_.end())
     {
         gridCoarse5d_[hkey5d].reset(
             SpaceTimeGrid::makeFiveDimGrid(cLs, gridCoarse4d_[hkey4d].get()));
+        HADRONS_DUMP_GRID(gridCoarse5d_[hkey5d].get());
     }
 }
 
+#undef HADRONS_DUMP_GRID
+
 template <typename VType>
-GridCartesian * Environment::getGrid(void) const
+GridCartesian * Environment::getGrid(void)
 {
-    auto it = grid4d_.find({typeHash<VType>(), 1});
+    FineGridKey key = {typeHash<VType>(), 1};
+
+    auto it = grid4d_.find(key);
 
     if (it != grid4d_.end())
     {
@@ -308,15 +326,17 @@ GridCartesian * Environment::getGrid(void) const
     }
     else
     {
-        HADRONS_ERROR(Definition, "no 4D grid for SIMD type '" 
-                                  + typeName<VType>() + "'");
+        createGrid<VType>(1);
+
+        return grid4d_.at(key).get();
     }
 }
 
 template <typename VType>
-GridRedBlackCartesian * Environment::getRbGrid(void) const
+GridRedBlackCartesian * Environment::getRbGrid(void)
 {
-    auto it = gridRb4d_.find({typeHash<VType>(), 1});
+    FineGridKey key = {typeHash<VType>(), 1};
+    auto        it  = gridRb4d_.find(key);
 
     if (it != gridRb4d_.end())
     {
@@ -324,34 +344,39 @@ GridRedBlackCartesian * Environment::getRbGrid(void) const
     }
     else
     {
-        HADRONS_ERROR(Definition, "no 4D red-black grid for SIMD type '" 
-                                  + typeName<VType>() + "'");
+        createGrid<VType>(1);
+
+        return gridRb4d_.at(key).get();
     }
 }
 
 template <typename VType>
-GridCartesian * Environment::getCoarseGrid(const std::vector<int> &blockSize) const
+GridCartesian * Environment::getCoarseGrid(const std::vector<int> &blockSize)
 {
     std::vector<int> s = blockSize;
 
     s.resize(getNd());
-    auto it = gridCoarse4d_.find({typeHash<VType>(), s});
+
+    CoarseGridKey key = {typeHash<VType>(), s};
+    auto          it  = gridCoarse4d_.find(key);
+
     if (it != gridCoarse4d_.end())
     {
         return it->second.get();
     }
     else
     {
-        HADRONS_ERROR(Definition, "no 4D coarse grid for SIMD type '" 
-                                  + typeName<VType>() + "' and block size "
-                                  + vecToStr(blockSize));
+        createCoarseGrid<VType>(blockSize, 1);
+        
+        return gridCoarse4d_.at(key).get();
     }
 }
 
 template <typename VType>
-GridCartesian * Environment::getGrid(const unsigned int Ls) const
+GridCartesian * Environment::getGrid(const unsigned int Ls)
 {
-    auto it = grid5d_.find({typeHash<VType>(), Ls});
+    FineGridKey key = {typeHash<VType>(), Ls};
+    auto        it  = grid5d_.find(key);
 
     if (it != grid5d_.end())
     {
@@ -359,16 +384,17 @@ GridCartesian * Environment::getGrid(const unsigned int Ls) const
     }
     else
     {
-        HADRONS_ERROR(Definition, "no 5D grid for SIMD type '" 
-                                  + typeName<VType>() + "' and Ls = "
-                                  + std::to_string(Ls));
+        createGrid<VType>(Ls);
+
+        return grid5d_.at(key).get();
     }
 }
 
 template <typename VType>
-GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
+GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls)
 {
-    auto it = gridRb5d_.find({typeHash<VType>(), Ls});
+    FineGridKey key = {typeHash<VType>(), Ls};
+    auto        it  = gridRb5d_.find(key);
 
     if (it != gridRb5d_.end())
     {
@@ -376,30 +402,32 @@ GridRedBlackCartesian * Environment::getRbGrid(const unsigned int Ls) const
     }
     else
     {
-        HADRONS_ERROR(Definition, "no 5D red-black grid for SIMD type '" 
-                                  + typeName<VType>() + "' and Ls = "
-                                  + std::to_string(Ls));
+        createGrid<VType>(Ls);
+
+        return gridRb5d_.at(key).get();
     }
 }
 
 template <typename VType>
 GridCartesian * Environment::getCoarseGrid(const std::vector<int> &blockSize,
-                                           const unsigned int Ls) const
+                                           const unsigned int Ls)
 {
     std::vector<int> s = blockSize;
 
     s.push_back(Ls);
-    auto it = gridCoarse5d_.find({typeHash<VType>(), s});
+
+    CoarseGridKey key = {typeHash<VType>(), s};
+
+    auto it = gridCoarse5d_.find(key);
     if (it != gridCoarse5d_.end())
     {
         return it->second.get();
     }
     else
     {
-        HADRONS_ERROR(Definition, "no 5D coarse grid for SIMD type '" 
-                                  + typeName<VType>() + "', block size "
-                                  + vecToStr(blockSize)
-                                  + " and Ls = " + std::to_string(Ls));
+        createCoarseGrid<VType>(blockSize, Ls);
+
+        return gridCoarse5d_.at(key).get();
     }
 }
 
diff --git a/Hadrons/Global.hpp b/Hadrons/Global.hpp
index 222990f6..e4003e73 100644
--- a/Hadrons/Global.hpp
+++ b/Hadrons/Global.hpp
@@ -106,7 +106,9 @@ typedef std::vector<SitePropagator##suffix>          SlicedPropagator##suffix;
 BASIC_TYPE_ALIASES(FImpl, suffix);\
 typedef FermionOperator<FImpl>            FMat##suffix;\
 typedef typename FImpl::FermionField      FermionField##suffix;\
-typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;
+typedef typename FImpl::GaugeField        GaugeField##suffix;\
+typedef typename FImpl::DoubledGaugeField DoubledGaugeField##suffix;\
+typedef typename FImpl::ComplexField      ComplexField##suffix;
 
 #define GAUGE_TYPE_ALIASES(GImpl, suffix)\
 typedef typename GImpl::GaugeField GaugeField##suffix;
diff --git a/Hadrons/Module.hpp b/Hadrons/Module.hpp
index e9c1c449..083040ce 100644
--- a/Hadrons/Module.hpp
+++ b/Hadrons/Module.hpp
@@ -65,14 +65,27 @@ extern template class base;\
 MODULE_REGISTER(mod, ARG(base), ns);
 
 #define ARG(...) __VA_ARGS__
-#define MACRO_REDIRECT(arg1, arg2, arg3, macro, ...) macro
+#define HADRONS_MACRO_REDIRECT_12(arg1, arg2, macro, ...) macro
+#define HADRONS_MACRO_REDIRECT_23(arg1, arg2, arg3, macro, ...) macro
 
-#define envCreateGrid(latticeType, Ls)\
-env().template createGrid<typename latticeType::vector_type>(Ls)
+#define envGetGrid4(latticeType)\
+env().template getGrid<typename latticeType::vector_type>()
 
-#define envGetGrid(latticeType, Ls)\
+#define envGetGrid5(latticeType, Ls)\
 env().template getGrid<typename latticeType::vector_type>(Ls)
 
+#define envGetGrid(...)\
+HADRONS_MACRO_REDIRECT_12(__VA_ARGS__, envGetGrid5, envGetGrid4)(__VA_ARGS__)
+
+#define envGetRbGrid4(latticeType)\
+env().template getRbGrid<typename latticeType::vector_type>()
+
+#define envGetRbGrid5(latticeType, Ls)\
+env().template getRbGrid<typename latticeType::vector_type>(Ls)
+
+#define envGetRbGrid(...)\
+HADRONS_MACRO_REDIRECT_12(__VA_ARGS__, envGetRbGrid5, envGetRbGrid4)(__VA_ARGS__)
+
 #define envGet(type, name)\
 *env().template getObject<type>(name)
 
@@ -92,38 +105,38 @@ env().template createObject<type>(name, Environment::Storage::object, Ls, __VA_A
 env().template createDerivedObject<base, type>(name, Environment::Storage::object, Ls, __VA_ARGS__)
 
 #define envCreateLat4(type, name)\
-envCreate(type, name, 1, env().template getGrid<typename type::vector_type>())
+envCreate(type, name, 1, envGetGrid(type))
 
 #define envCreateLat5(type, name, Ls)\
-envCreate(type, name, Ls, env().template getGrid<typename type::vector_type>(Ls))
+envCreate(type, name, Ls, envGetGrid(type, Ls))
 
 #define envCreateLat(...)\
-MACRO_REDIRECT(__VA_ARGS__, envCreateLat5, envCreateLat4)(__VA_ARGS__)
+HADRONS_MACRO_REDIRECT_23(__VA_ARGS__, envCreateLat5, envCreateLat4)(__VA_ARGS__)
 
 #define envCache(type, name, Ls, ...)\
 env().template createObject<type>(name, Environment::Storage::cache, Ls, __VA_ARGS__)
 
 #define envCacheLat4(type, name)\
-envCache(type, name, 1, env().template getGrid<typename type::vector_type>())
+envCache(type, name, 1, envGetGrid(type))
 
 #define envCacheLat5(type, name, Ls)\
-envCache(type, name, Ls, env().template getGrid<typename type::vector_type>(Ls))
+envCache(type, name, Ls, envGetGrid(type, Ls))
 
 #define envCacheLat(...)\
-MACRO_REDIRECT(__VA_ARGS__, envCacheLat5, envCacheLat4)(__VA_ARGS__)
+HADRONS_MACRO_REDIRECT_23(__VA_ARGS__, envCacheLat5, envCacheLat4)(__VA_ARGS__)
 
 #define envTmp(type, name, Ls, ...)\
 env().template createObject<type>(getName() + "_tmp_" + name,         \
                                   Environment::Storage::temporary, Ls, __VA_ARGS__)
 
 #define envTmpLat4(type, name)\
-envTmp(type, name, 1, env().template getGrid<typename type::vector_type>())
+envTmp(type, name, 1, envGetGrid(type))
 
 #define envTmpLat5(type, name, Ls)\
-envTmp(type, name, Ls, env().template getGrid<typename type::vector_type>(Ls))
+envTmp(type, name, Ls, envGetGrid(type, Ls))
 
 #define envTmpLat(...)\
-MACRO_REDIRECT(__VA_ARGS__, envTmpLat5, envTmpLat4)(__VA_ARGS__)
+HADRONS_MACRO_REDIRECT_23(__VA_ARGS__, envTmpLat5, envTmpLat4)(__VA_ARGS__)
 
 #define saveResult(ioStem, name, result)\
 if (env().getGrid()->IsBoss() and !ioStem.empty())\
diff --git a/Hadrons/Modules/MAction/DWF.cc b/Hadrons/Modules/MAction/DWF.cc
index 35ff60ac..6cda84ac 100644
--- a/Hadrons/Modules/MAction/DWF.cc
+++ b/Hadrons/Modules/MAction/DWF.cc
@@ -32,4 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TDWF<FIMPL>;
-
+template class Grid::Hadrons::MAction::TDWF<FIMPLF>;
diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp
index 714bccc3..67cfeb1b 100644
--- a/Hadrons/Modules/MAction/DWF.hpp
+++ b/Hadrons/Modules/MAction/DWF.hpp
@@ -73,6 +73,7 @@ protected:
 };
 
 MODULE_REGISTER_TMP(DWF, TDWF<FIMPL>, MAction);
+MODULE_REGISTER_TMP(DWFF, TDWF<FIMPLF>, MAction);
 
 /******************************************************************************
  *                        DWF template implementation                         *
@@ -111,12 +112,11 @@ void TDWF<FImpl>::setup(void)
     LOG(Message) << "Fermion boundary conditions: " << par().boundary
                  << std::endl;
                  
-    env().createGrid(par().Ls);
-    auto &U    = envGet(LatticeGaugeField, par().gauge);
-    auto &g4   = *env().getGrid();
-    auto &grb4 = *env().getRbGrid();
-    auto &g5   = *env().getGrid(par().Ls);
-    auto &grb5 = *env().getRbGrid(par().Ls);
+    auto &U    = envGet(GaugeField, par().gauge);
+    auto &g4   = *envGetGrid(FermionField);
+    auto &grb4 = *envGetRbGrid(FermionField);
+    auto &g5   = *envGetGrid(FermionField, par().Ls);
+    auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
     typename DomainWallFermion<FImpl>::ImplParams implParams(boundary);
     envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
diff --git a/Hadrons/Modules/MAction/MobiusDWF.cc b/Hadrons/Modules/MAction/MobiusDWF.cc
index ba45462f..41683920 100644
--- a/Hadrons/Modules/MAction/MobiusDWF.cc
+++ b/Hadrons/Modules/MAction/MobiusDWF.cc
@@ -32,3 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TMobiusDWF<FIMPL>;
+template class Grid::Hadrons::MAction::TMobiusDWF<FIMPLF>;
diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp
index 4db2ebb6..22964c9a 100644
--- a/Hadrons/Modules/MAction/MobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/MobiusDWF.hpp
@@ -72,6 +72,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(MobiusDWF, TMobiusDWF<FIMPL>, MAction);
+MODULE_REGISTER_TMP(MobiusDWFF, TMobiusDWF<FIMPLF>, MAction);
 
 /******************************************************************************
  *                      TMobiusDWF implementation                             *
@@ -111,12 +112,11 @@ void TMobiusDWF<FImpl>::setup(void)
     LOG(Message) << "Fermion boundary conditions: " << par().boundary
                  << std::endl;
 
-    env().createGrid(par().Ls);
-    auto &U    = envGet(LatticeGaugeField, par().gauge);
-    auto &g4   = *env().getGrid();
-    auto &grb4 = *env().getRbGrid();
-    auto &g5   = *env().getGrid(par().Ls);
-    auto &grb5 = *env().getRbGrid(par().Ls);
+    auto &U    = envGet(GaugeField, par().gauge);
+    auto &g4   = *envGetGrid(FermionField);
+    auto &grb4 = *envGetRbGrid(FermionField);
+    auto &g5   = *envGetGrid(FermionField, par().Ls);
+    auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
     typename MobiusFermion<FImpl>::ImplParams implParams(boundary);
     envCreateDerived(FMat, MobiusFermion<FImpl>, getName(), par().Ls, U, g5,
diff --git a/Hadrons/Modules/MAction/ScaledDWF.cc b/Hadrons/Modules/MAction/ScaledDWF.cc
index b3f7acdb..70eafac5 100644
--- a/Hadrons/Modules/MAction/ScaledDWF.cc
+++ b/Hadrons/Modules/MAction/ScaledDWF.cc
@@ -32,3 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TScaledDWF<FIMPL>;
+template class Grid::Hadrons::MAction::TScaledDWF<FIMPLF>;
diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp
index 850a8662..c890f0e9 100644
--- a/Hadrons/Modules/MAction/ScaledDWF.hpp
+++ b/Hadrons/Modules/MAction/ScaledDWF.hpp
@@ -71,6 +71,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(ScaledDWF, TScaledDWF<FIMPL>, MAction);
+MODULE_REGISTER_TMP(ScaledDWFF, TScaledDWF<FIMPLF>, MAction);
 
 /******************************************************************************
  *                      TScaledDWF implementation                             *
@@ -110,12 +111,11 @@ void TScaledDWF<FImpl>::setup(void)
     LOG(Message) << "Fermion boundary conditions: " << par().boundary
                  << std::endl;
 
-    env().createGrid(par().Ls);
-    auto &U    = envGet(LatticeGaugeField, par().gauge);
-    auto &g4   = *env().getGrid();
-    auto &grb4 = *env().getRbGrid();
-    auto &g5   = *env().getGrid(par().Ls);
-    auto &grb5 = *env().getRbGrid(par().Ls);
+    auto &U    = envGet(GaugeField, par().gauge);
+    auto &g4   = *envGetGrid(FermionField);
+    auto &grb4 = *envGetRbGrid(FermionField);
+    auto &g5   = *envGetGrid(FermionField, par().Ls);
+    auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
     typename MobiusFermion<FImpl>::ImplParams implParams(boundary);
     envCreateDerived(FMat, ScaledShamirFermion<FImpl>, getName(), par().Ls, U, g5,
diff --git a/Hadrons/Modules/MAction/Wilson.cc b/Hadrons/Modules/MAction/Wilson.cc
index faceb76e..4ce3e7ef 100644
--- a/Hadrons/Modules/MAction/Wilson.cc
+++ b/Hadrons/Modules/MAction/Wilson.cc
@@ -32,4 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TWilson<FIMPL>;
-
+template class Grid::Hadrons::MAction::TWilson<FIMPLF>;
diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp
index cccd648f..a9327c2f 100644
--- a/Hadrons/Modules/MAction/Wilson.hpp
+++ b/Hadrons/Modules/MAction/Wilson.hpp
@@ -71,6 +71,7 @@ protected:
 };
 
 MODULE_REGISTER_TMP(Wilson, TWilson<FIMPL>, MAction);
+MODULE_REGISTER_TMP(WilsonF, TWilson<FIMPLF>, MAction);
 
 /******************************************************************************
  *                     TWilson template implementation                        *
@@ -107,9 +108,9 @@ void TWilson<FImpl>::setup(void)
     LOG(Message) << "Fermion boundary conditions: " << par().boundary
                  << std::endl;
                  
-    auto &U      = envGet(LatticeGaugeField, par().gauge);
-    auto &grid   = *env().getGrid();
-    auto &gridRb = *env().getRbGrid();
+    auto &U      = envGet(GaugeField, par().gauge);
+    auto &grid   = *envGetGrid(FermionField);
+    auto &gridRb = *envGetRbGrid(FermionField);
     std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
     typename WilsonFermion<FImpl>::ImplParams implParams(boundary);
     envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
diff --git a/Hadrons/Modules/MAction/WilsonClover.cc b/Hadrons/Modules/MAction/WilsonClover.cc
index 601388aa..2c5c0e66 100644
--- a/Hadrons/Modules/MAction/WilsonClover.cc
+++ b/Hadrons/Modules/MAction/WilsonClover.cc
@@ -32,4 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TWilsonClover<FIMPL>;
-
+template class Grid::Hadrons::MAction::TWilsonClover<FIMPLF>;
diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp
index 48f4f44e..349abe84 100644
--- a/Hadrons/Modules/MAction/WilsonClover.hpp
+++ b/Hadrons/Modules/MAction/WilsonClover.hpp
@@ -75,9 +75,10 @@ public:
 };
 
 MODULE_REGISTER_TMP(WilsonClover, TWilsonClover<FIMPL>, MAction);
+MODULE_REGISTER_TMP(WilsonCloverF, TWilsonClover<FIMPLF>, MAction);
 
 /******************************************************************************
- *                     TWilsonClover template implementation                        *
+ *                    TWilsonClover template implementation                   *
  ******************************************************************************/
 // constructor /////////////////////////////////////////////////////////////////
 template <typename FImpl>
@@ -113,16 +114,14 @@ void TWilsonClover<FImpl>::setup(void)
     LOG(Message) << "Clover term csw_r: " << par().csw_r
                  << " csw_t: " << par().csw_t
                  << std::endl;
-    auto &U      = envGet(LatticeGaugeField, par().gauge);
-    auto &grid   = *env().getGrid();
-    auto &gridRb = *env().getRbGrid();
+    auto &U      = envGet(GaugeField, par().gauge);
+    auto &grid   = *envGetGrid(FermionField);
+    auto &gridRb = *envGetRbGrid(FermionField);
     std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
     typename WilsonCloverFermion<FImpl>::ImplParams implParams(boundary);
-    envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid, gridRb, par().mass,
-						  par().csw_r,
-						  par().csw_t,
-					      par().clover_anisotropy,
-						  implParams); 
+    envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid,
+                     gridRb, par().mass, par().csw_r, par().csw_t, 
+                     par().clover_anisotropy, implParams); 
 }
 
 // execution ///////////////////////////////////////////////////////////////////
diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.cc b/Hadrons/Modules/MAction/ZMobiusDWF.cc
index 0eb06130..ef8e4799 100644
--- a/Hadrons/Modules/MAction/ZMobiusDWF.cc
+++ b/Hadrons/Modules/MAction/ZMobiusDWF.cc
@@ -32,4 +32,4 @@ using namespace Hadrons;
 using namespace MAction;
 
 template class Grid::Hadrons::MAction::TZMobiusDWF<ZFIMPL>;
-
+template class Grid::Hadrons::MAction::TZMobiusDWF<ZFIMPLF>;
diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
index 6a02a1ec..f7959127 100644
--- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
@@ -73,6 +73,7 @@ public:
 };
 
 MODULE_REGISTER_TMP(ZMobiusDWF, TZMobiusDWF<ZFIMPL>, MAction);
+MODULE_REGISTER_TMP(ZMobiusDWFF, TZMobiusDWF<ZFIMPLF>, MAction);
 
 /******************************************************************************
  *                     TZMobiusDWF implementation                             *
@@ -118,11 +119,11 @@ void TZMobiusDWF<FImpl>::setup(void)
                  << std::endl;
 
     env().createGrid(par().Ls);
-    auto &U    = envGet(LatticeGaugeField, par().gauge);
-    auto &g4   = *env().getGrid();
-    auto &grb4 = *env().getRbGrid();
-    auto &g5   = *env().getGrid(par().Ls);
-    auto &grb5 = *env().getRbGrid(par().Ls);
+    auto &U    = envGet(GaugeField, par().gauge);
+    auto &g4   = *envGetGrid(FermionField);
+    auto &grb4 = *envGetRbGrid(FermionField);
+    auto &g5   = *envGetGrid(FermionField, par().Ls);
+    auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     auto omega = par().omega;
     std::vector<Complex> boundary = strToVec<Complex>(par().boundary);
     typename ZMobiusFermion<FImpl>::ImplParams implParams(boundary);
diff --git a/Hadrons/Modules/MContraction/A2AMesonField.hpp b/Hadrons/Modules/MContraction/A2AMesonField.hpp
index 1a922b49..9e07ec6c 100644
--- a/Hadrons/Modules/MContraction/A2AMesonField.hpp
+++ b/Hadrons/Modules/MContraction/A2AMesonField.hpp
@@ -187,9 +187,9 @@ void TA2AMesonField<FImpl>::setup(void)
         mom_.push_back(p);
     }
     
-    envCache(std::vector<LatticeComplex>, momphName_, 1, 
-             par().mom.size(), env().getGrid());
-    envTmpLat(LatticeComplex, "coor");
+    envCache(std::vector<ComplexField>, momphName_, 1, 
+             par().mom.size(), envGetGrid(ComplexField));
+    envTmpLat(ComplexField, "coor");
     // preallocate memory for meson field block
     auto tgp = env().getDim().back()*gamma_.size()*mom_.size();
 
@@ -231,7 +231,7 @@ void TA2AMesonField<FImpl>::execute(void)
     ///////////////////////////////////////////////
     // Momentum setup
     ///////////////////////////////////////////////
-    auto &ph = envGet(std::vector<LatticeComplex>, momphName_);
+    auto &ph = envGet(std::vector<ComplexField>, momphName_);
 
     if (!hasPhase_)
     {
@@ -241,7 +241,7 @@ void TA2AMesonField<FImpl>::execute(void)
             Complex           i(0.0,1.0);
             std::vector<Real> p;
 
-            envGetTmp(LatticeComplex, coor);
+            envGetTmp(ComplexField, coor);
             ph[j] = zero;
             for(unsigned int mu = 0; mu < mom_[j].size(); mu++)
             {
diff --git a/Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp b/Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp
index 7decb142..da634d6c 100644
--- a/Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp
+++ b/Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp
@@ -43,6 +43,8 @@ BEGIN_MODULE_NAMESPACE(MNoise)
 template <typename FImpl>
 class TTimeDilutedSpinColorDiagonal: public Module<NoPar>
 {
+public:
+    FERM_TYPE_ALIASES(FImpl,);
 public:
     // constructor
     TTimeDilutedSpinColorDiagonal(const std::string name);
@@ -92,7 +94,7 @@ void TTimeDilutedSpinColorDiagonal<FImpl>::setup(void)
 {
     envCreateDerived(DilutedNoise<FImpl>, 
                      TimeDilutedSpinColorDiagonalNoise<FImpl>,
-                     getName(), 1, env().getGrid());
+                     getName(), 1, envGetGrid(FermionField));
 }
 
 // execution ///////////////////////////////////////////////////////////////////
diff --git a/Hadrons/Modules/MScalarSUN/StochFreeField.hpp b/Hadrons/Modules/MScalarSUN/StochFreeField.hpp
index 65b6bf21..141eeafd 100644
--- a/Hadrons/Modules/MScalarSUN/StochFreeField.hpp
+++ b/Hadrons/Modules/MScalarSUN/StochFreeField.hpp
@@ -132,7 +132,7 @@ void TStochFreeField<SImpl>::execute(void)
     auto    &w                = envGet(ComplexField, "_" + getName() + "_weight");
     auto    &rng              = rng4d();
     double  trphi2;
-    FFT     fft(env().getGrid());
+    FFT     fft(envGetGrid(Field));
     Integer vol;
 
     vol = 1;
@@ -169,11 +169,6 @@ void TStochFreeField<SImpl>::execute(void)
     phi = 0.5*(phi - adj(phi));
     trphi2 = -TensorRemove(sum(trace(phi*phi))).real()/vol;
     LOG(Message) << "tr(phi^2)= " << trphi2 << std::endl;
-
-    // ComplexField phi2(env().getGrid());
-
-    // phi2=trace(phi*phi);
-    // std::cout << phi2 << std::endl;
 }
 
 END_MODULE_NAMESPACE
diff --git a/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp b/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp
index 2b31fbe9..a73da8de 100644
--- a/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp
+++ b/Hadrons/Modules/MScalarSUN/TimeMomProbe.hpp
@@ -146,7 +146,7 @@ void TTimeMomProbe<SImpl>::execute(void)
     std::set<std::vector<int>>                   timeMomSet;
     std::vector<std::vector<std::vector<int>>>   timeMom;
     std::vector<std::vector<int>>                transferMom;
-    FFT                                          fft(env().getGrid());
+    FFT                                          fft(envGetGrid(Field));
     std::vector<int>                             dMask(nd, 1);
     std::vector<TimeMomProbeResult>              result;
     std::map<std::string, std::vector<SlicedOp>> slicedOp;
diff --git a/Hadrons/Modules/MScalarSUN/TrKinetic.hpp b/Hadrons/Modules/MScalarSUN/TrKinetic.hpp
index 0fd931d6..8c82075f 100644
--- a/Hadrons/Modules/MScalarSUN/TrKinetic.hpp
+++ b/Hadrons/Modules/MScalarSUN/TrKinetic.hpp
@@ -126,7 +126,7 @@ void TTrKinetic<SImpl>::setup(void)
         envCreateLat(ComplexField, varName(getName(), mu, nu));
     }
     envCreateLat(ComplexField, varName(getName(), "sum"));
-    envTmp(std::vector<Field>, "der", 1, env().getNd(), env().getGrid());
+    envTmp(std::vector<Field>, "der", 1, env().getNd(), envGetGrid(Field));
 }
 
 // execution ///////////////////////////////////////////////////////////////////
diff --git a/Hadrons/Modules/MScalarSUN/TwoPoint.hpp b/Hadrons/Modules/MScalarSUN/TwoPoint.hpp
index 5cb7e867..df0d2cba 100644
--- a/Hadrons/Modules/MScalarSUN/TwoPoint.hpp
+++ b/Hadrons/Modules/MScalarSUN/TwoPoint.hpp
@@ -168,7 +168,7 @@ void TTwoPoint<SImpl>::execute(void)
     std::set<std::string>                        ops;
     std::vector<TwoPointResult>                  result;
     std::map<std::string, std::vector<SlicedOp>> slicedOp;
-    FFT                                          fft(env().getGrid());
+    FFT                                          fft(envGetGrid(Field));
     TComplex                                     buf;
 
     envGetTmp(ComplexField, ftBuf);
diff --git a/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp b/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp
index e13ca8b4..6a18f6cb 100644
--- a/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp
+++ b/Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp
@@ -136,7 +136,7 @@ void TTwoPointNPR<SImpl>::execute(void)
     const unsigned int             nd   = env().getNd();
     const unsigned int             nl   = env().getDim(0);
     const Real                     invV = 1./env().getVolume();
-    FFT                            fft(env().getGrid());
+    FFT                            fft(envGetGrid(Field));
     std::vector<TwoPointNPRResult> result;
     TwoPointNPRResult              twoPtp1, twoPtp2, twoPtDisc;
     auto                           &phi    = envGet(Field, par().field);
diff --git a/Hadrons/Modules/MSink/Point.hpp b/Hadrons/Modules/MSink/Point.hpp
index f7970346..02189296 100644
--- a/Hadrons/Modules/MSink/Point.hpp
+++ b/Hadrons/Modules/MSink/Point.hpp
@@ -132,7 +132,7 @@ void TPoint<FImpl>::execute(void)
         for(unsigned int mu = 0; mu < p.size(); mu++)
         {
             LatticeCoordinate(coor, mu);
-            ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
+            ph = ph + (p[mu]/env().getDim(mu))*coor;
         }
         ph = exp((Real)(2*M_PI)*i*ph);
         hasPhase_ = true;
diff --git a/Hadrons/Modules/MSolver/A2AVectors.hpp b/Hadrons/Modules/MSolver/A2AVectors.hpp
index 1de00142..608e05b9 100644
--- a/Hadrons/Modules/MSolver/A2AVectors.hpp
+++ b/Hadrons/Modules/MSolver/A2AVectors.hpp
@@ -135,9 +135,9 @@ void TA2AVectors<FImpl, Pack>::setup(void)
         Nl_ = epack.evec.size();
     }
     envCreate(std::vector<FermionField>, getName() + "_v", 1, 
-              Nl_ + noise.size(), FermionField(env().getGrid()));
+              Nl_ + noise.size(), envGetGrid(FermionField));
     envCreate(std::vector<FermionField>, getName() + "_w", 1, 
-              Nl_ + noise.size(), FermionField(env().getGrid()));
+              Nl_ + noise.size(), envGetGrid(FermionField));
     if (Ls > 1)
     {
         envTmpLat(FermionField, "f5", Ls);
diff --git a/Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp b/Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp
index 5be1e3ce..67e1be19 100644
--- a/Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp
+++ b/Hadrons/Modules/MSolver/MixedPrecisionRBPrecCG.hpp
@@ -170,7 +170,7 @@ void TMixedPrecisionRBPrecCG<FImplInner, FImplOuter, nBasis>
             MixedPrecisionConjugateGradient<FermionFieldOuter, FermionFieldInner> 
                 mpcg(par().residual, par().maxInnerIteration, 
                      par().maxOuterIteration, 
-                     env().template getGrid<VTypeInner>(Ls),
+                     env().template getRbGrid<VTypeInner>(Ls),
                      simat, somat);
             OperatorFunctionWrapper<FermionFieldOuter> wmpcg(mpcg);
             HADRONS_DEFAULT_SCHUR_SOLVE<FermionFieldOuter> schurSolver(wmpcg);
diff --git a/Hadrons/Modules/MSource/SeqConserved.hpp b/Hadrons/Modules/MSource/SeqConserved.hpp
index 48a811d2..d34e427e 100644
--- a/Hadrons/Modules/MSource/SeqConserved.hpp
+++ b/Hadrons/Modules/MSource/SeqConserved.hpp
@@ -187,7 +187,7 @@ void TSeqConserved<FImpl>::execute(void)
         for(unsigned int mu = 0; mu < env().getNd(); mu++)
         {
             LatticeCoordinate(coor, mu);
-            mom_phase = mom_phase + (mom[mu]/env().getGrid()->_fdimensions[mu])*coor;
+            mom_phase = mom_phase + (mom[mu]/env().getDim(mu))*coor;
         }
         mom_phase = exp((Real)(2*M_PI)*i*mom_phase);
         SeqhasPhase_ = true;
diff --git a/Hadrons/Modules/MSource/SeqGamma.hpp b/Hadrons/Modules/MSource/SeqGamma.hpp
index 59fdb7ae..b9fb9b1b 100644
--- a/Hadrons/Modules/MSource/SeqGamma.hpp
+++ b/Hadrons/Modules/MSource/SeqGamma.hpp
@@ -125,7 +125,7 @@ template <typename FImpl>
 void TSeqGamma<FImpl>::setup(void)
 {
     envCreateLat(PropagatorField, getName());
-    envCache(Lattice<iScalar<vInteger>>, tName_, 1, env().getGrid());
+    envCache(Lattice<iScalar<vInteger>>, tName_, 1, envGetGrid(LatticeComplex));
     envCacheLat(LatticeComplex, momphName_);
     envTmpLat(LatticeComplex, "coor");
 }
@@ -162,7 +162,7 @@ void TSeqGamma<FImpl>::execute(void)
         for(unsigned int mu = 0; mu < env().getNd(); mu++)
         {
             LatticeCoordinate(coor, mu);
-            ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
+            ph = ph + (p[mu]/env().getDim(mu))*coor;
         }
         ph = exp((Real)(2*M_PI)*i*ph);
         LatticeCoordinate(t, Tp);
diff --git a/Hadrons/Modules/MSource/Wall.hpp b/Hadrons/Modules/MSource/Wall.hpp
index c4b99448..3dfbc199 100644
--- a/Hadrons/Modules/MSource/Wall.hpp
+++ b/Hadrons/Modules/MSource/Wall.hpp
@@ -143,7 +143,7 @@ void TWall<FImpl>::execute(void)
         for(unsigned int mu = 0; mu < env().getNd(); mu++)
         {
             LatticeCoordinate(coor, mu);
-            ph = ph + (p[mu]/env().getGrid()->_fdimensions[mu])*coor;
+            ph = ph + (p[mu]/env().getDim(mu))*coor;
         }
         ph = exp((Real)(2*M_PI)*i*ph);
         LatticeCoordinate(t, Tp);
diff --git a/Hadrons/Modules/MSource/Z2.hpp b/Hadrons/Modules/MSource/Z2.hpp
index ec57fbd4..e16114e0 100644
--- a/Hadrons/Modules/MSource/Z2.hpp
+++ b/Hadrons/Modules/MSource/Z2.hpp
@@ -120,7 +120,7 @@ template <typename FImpl>
 void TZ2<FImpl>::setup(void)
 {
     envCreateLat(PropagatorField, getName());
-    envCache(Lattice<iScalar<vInteger>>, tName_, 1, env().getGrid());
+    envCache(Lattice<iScalar<vInteger>>, tName_, 1, envGetGrid(LatticeComplex));
     envTmpLat(LatticeComplex, "eta");
 }
 
diff --git a/Hadrons/Modules/MUtilities/RandomVectors.hpp b/Hadrons/Modules/MUtilities/RandomVectors.hpp
index 2444b850..65126eb7 100644
--- a/Hadrons/Modules/MUtilities/RandomVectors.hpp
+++ b/Hadrons/Modules/MUtilities/RandomVectors.hpp
@@ -99,11 +99,11 @@ void TRandomVectors<Field>::setup(void)
     if (par().Ls > 1)
     {
         envCreate(std::vector<Field>, getName(), par().Ls, par().size, 
-                  env().getGrid(par().Ls));
+                  envGetGrid(Field, par().Ls));
     }
     else
     {
-        envCreate(std::vector<Field>, getName(), 1, par().size, env().getGrid());
+        envCreate(std::vector<Field>, getName(), 1, par().size, envGetGrid(Field));
     }
 }