diff --git a/Hadrons/Modules/MAction/DWF.hpp b/Hadrons/Modules/MAction/DWF.hpp
index bac70ae8..21024a4f 100644
--- a/Hadrons/Modules/MAction/DWF.hpp
+++ b/Hadrons/Modules/MAction/DWF.hpp
@@ -112,8 +112,6 @@ void TDWF<FImpl>::setup(void)
                  << par().mass << ", M5= " << par().M5 << " and Ls= "
                  << par().Ls << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                  
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
@@ -121,8 +119,26 @@ void TDWF<FImpl>::setup(void)
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     typename DomainWallFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, DomainWallFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, implParams);
 }
diff --git a/Hadrons/Modules/MAction/MobiusDWF.hpp b/Hadrons/Modules/MAction/MobiusDWF.hpp
index 6b960e2e..3ce02b91 100644
--- a/Hadrons/Modules/MAction/MobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/MobiusDWF.hpp
@@ -112,17 +112,33 @@ void TMobiusDWF<FImpl>::setup(void)
                  << ", b= " << par().b << ", c= " << par().c
                  << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
-
+                 
     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);
     typename MobiusFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, MobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().b, par().c,
                      implParams);
diff --git a/Hadrons/Modules/MAction/ScaledDWF.hpp b/Hadrons/Modules/MAction/ScaledDWF.hpp
index 8742a820..3e052225 100644
--- a/Hadrons/Modules/MAction/ScaledDWF.hpp
+++ b/Hadrons/Modules/MAction/ScaledDWF.hpp
@@ -111,8 +111,6 @@ void TScaledDWF<FImpl>::setup(void)
                  << ", scale= " << par().scale
                  << " using gauge field '" << par().gauge << "'"
                  << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
 
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
@@ -120,8 +118,26 @@ void TScaledDWF<FImpl>::setup(void)
     auto &g5   = *envGetGrid(FermionField, par().Ls);
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     typename ScaledShamirFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, ScaledShamirFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, par().scale,
                      implParams);
diff --git a/Hadrons/Modules/MAction/Wilson.hpp b/Hadrons/Modules/MAction/Wilson.hpp
index 3ce5569e..11df89c3 100644
--- a/Hadrons/Modules/MAction/Wilson.hpp
+++ b/Hadrons/Modules/MAction/Wilson.hpp
@@ -109,15 +109,31 @@ void TWilson<FImpl>::setup(void)
 {
     LOG(Message) << "Setting up Wilson fermion matrix with m= " << par().mass
                  << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
                  
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
     typename WilsonFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, WilsonFermion<FImpl>, getName(), 1, U, grid, gridRb,
                      par().mass, implParams);
 }
diff --git a/Hadrons/Modules/MAction/WilsonClover.hpp b/Hadrons/Modules/MAction/WilsonClover.hpp
index a9ea8de5..d80bccdf 100644
--- a/Hadrons/Modules/MAction/WilsonClover.hpp
+++ b/Hadrons/Modules/MAction/WilsonClover.hpp
@@ -112,17 +112,34 @@ void TWilsonClover<FImpl>::setup(void)
 {
     LOG(Message) << "Setting up Wilson clover fermion matrix with m= " << par().mass
                  << " using gauge field '" << par().gauge << "'" << std::endl;
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary 
-                 << std::endl;
     LOG(Message) << "Clover term csw_r: " << par().csw_r
                  << " csw_t: " << par().csw_t
                  << std::endl;
+                 
     auto &U      = envGet(GaugeField, par().gauge);
     auto &grid   = *envGetGrid(FermionField);
     auto &gridRb = *envGetRbGrid(FermionField);
     typename WilsonCloverFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, WilsonCloverFermion<FImpl>, getName(), 1, U, grid,
                      gridRb, par().mass, par().csw_r, par().csw_t, 
                      par().clover_anisotropy, implParams); 
diff --git a/Hadrons/Modules/MAction/ZMobiusDWF.hpp b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
index e780bd73..56efe1a2 100644
--- a/Hadrons/Modules/MAction/ZMobiusDWF.hpp
+++ b/Hadrons/Modules/MAction/ZMobiusDWF.hpp
@@ -118,10 +118,7 @@ void TZMobiusDWF<FImpl>::setup(void)
     {
         LOG(Message) << "  omega[" << i << "]= " << par().omega[i] << std::endl;
     }
-    LOG(Message) << "Fermion boundary conditions: " << par().boundary
-                 << std::endl;
 
-    env().createGrid(par().Ls);
     auto &U    = envGet(GaugeField, par().gauge);
     auto &g4   = *envGetGrid(FermionField);
     auto &grb4 = *envGetRbGrid(FermionField);
@@ -129,8 +126,26 @@ void TZMobiusDWF<FImpl>::setup(void)
     auto &grb5 = *envGetRbGrid(FermionField, par().Ls);
     auto omega = par().omega;
     typename ZMobiusFermion<FImpl>::ImplParams implParams;
-    implParams.boundary_phases = strToVec<Complex>(par().boundary);
-    implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    if (!par().boundary.empty())
+    {
+        implParams.boundary_phases = strToVec<Complex>(par().boundary);
+    }
+    if (!par().twist.empty())
+    {
+        implParams.twist_n_2pi_L   = strToVec<Real>(par().twist);
+    }
+    LOG(Message) << "Fermion boundary conditions: " << implParams.boundary_phases
+                 << std::endl;
+    LOG(Message) << "Twists: " << implParams.twist_n_2pi_L
+                 << std::endl;
+    if (implParams.boundary_phases.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of boundary phase");
+    }
+    if (implParams.twist_n_2pi_L.size() != env().getNd())
+    {
+        HADRONS_ERROR(Size, "Wrong number of twist");
+    }
     envCreateDerived(FMat, ZMobiusFermion<FImpl>, getName(), par().Ls, U, g5,
                      grb5, g4, grb4, par().mass, par().M5, omega,
                      par().b, par().c, implParams);