diff --git a/examples/PinnedH2O/job.basis_1_50 b/examples/PinnedH2O/job.basis_1_50
index 5291a7b7..67b0d61f 100644
--- a/examples/PinnedH2O/job.basis_1_50
+++ b/examples/PinnedH2O/job.basis_1_50
@@ -8,7 +8,7 @@ date
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 8
+set ncpus = 1
 
 set maindir = /p/lustre2/cheung26/mgmol
 
diff --git a/examples/PinnedH2O/job.offline b/examples/PinnedH2O/job.offline
index 9d7150f8..b6b55e64 100644
--- a/examples/PinnedH2O/job.offline
+++ b/examples/PinnedH2O/job.offline
@@ -8,7 +8,7 @@ date
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 8
+set ncpus = 1
 
 set maindir = /p/lustre2/cheung26/mgmol
 
diff --git a/examples/PinnedH2O/job.rom b/examples/PinnedH2O/job.rom
index 45285136..75a97bf4 100644
--- a/examples/PinnedH2O/job.rom
+++ b/examples/PinnedH2O/job.rom
@@ -8,7 +8,7 @@ date
 setenv OMP_NUM_THREADS 1
 #setenv KMP_DETERMINISTIC_REDUCTION 1
 
-set ncpus = 8
+set ncpus = 1
 
 set maindir = /p/lustre2/cheung26/mgmol
 
diff --git a/src/Control.cc b/src/Control.cc
index 94e2ff30..01843fa9 100644
--- a/src/Control.cc
+++ b/src/Control.cc
@@ -269,7 +269,9 @@ void Control::print(std::ostream& os)
     if (precond_type_ == 10)
     {
         os << " Multigrid preconditioning for wave functions:" << std::endl;
-        os << " # of Multigrid levels   : " << mg_levels_ << std::endl;
+        os << " # of Multigrid levels     : " << mg_levels_ << std::endl;
+        os << " # of pre-smoothing steps  : " << mg_npresmoothing_ << std::endl;
+        os << " # of post-smoothing steps : " << mg_npostsmoothing_ << std::endl;
     }
     else
     {
@@ -334,7 +336,7 @@ void Control::sync(void)
     if (onpe0 && verbose > 0)
         (*MPIdata::sout) << "Control::sync()" << std::endl;
     // pack
-    const short size_short_buffer = 93;
+    const short size_short_buffer = 95;
     short* short_buffer           = new short[size_short_buffer];
     if (mype_ == 0)
     {
@@ -426,6 +428,8 @@ void Control::sync(void)
         short_buffer[90] = (short)static_cast<int>(poisson_lap_type_);
         short_buffer[91] = poisson_pc_data_;
         short_buffer[92] = precond_precision_;
+        short_buffer[93] = mg_npresmoothing_;
+        short_buffer[94] = mg_npostsmoothing_;
     }
     else
     {
@@ -642,6 +646,8 @@ void Control::sync(void)
     poisson_lap_type_  = static_cast<PoissonFDtype>(short_buffer[90]);
     poisson_pc_data_   = short_buffer[91];
     precond_precision_ = short_buffer[92];
+    mg_npresmoothing_  = short_buffer[93];
+    mg_npostsmoothing_ = short_buffer[94];
 
     numst    = int_buffer[0];
     nel_     = int_buffer[1];
@@ -1486,8 +1492,13 @@ void Control::setOptions(const boost::program_options::variables_map& vm)
         std::cout << "Outer solver type: " << str << std::endl;
         assert(it_algo_type_ >= 0);
 
-        mg_levels_ = vm["Quench.preconditioner_num_levels"].as<short>() - 1;
-        precond_precision_ = vm["Quench.preconditioner_precision"].as<short>();
+        // Preconditioner parameters
+        mg_levels_         = vm["Preconditioner.num_levels"].as<short>() - 1;
+        mg_npresmoothing_  = vm["Preconditioner.npresmoothing"].as<short>();
+        mg_npostsmoothing_ = vm["Preconditioner.npostsmoothing"].as<short>();
+        precond_precision_ = vm["Preconditioner.precision"].as<short>();
+	assert(precond_precision_==32 ||precond_precision_==64);
+
         precond_factor     = vm["Quench.step_length"].as<float>();
         if (precond_factor < 0.)
         {
diff --git a/src/Control.h b/src/Control.h
index 0956eb07..6ae93692 100644
--- a/src/Control.h
+++ b/src/Control.h
@@ -176,9 +176,6 @@ class Control
     // 11=local greedy
     short coloring_algo_;
 
-    // Number of MG levels for preconditioning
-    short mg_levels_;
-
     // preconditioning type
     // 10 = MG, block implementation
     short precond_type_;
@@ -404,6 +401,12 @@ class Control
 
     float betaAnderson;
 
+    // Number of MG levels for preconditioning
+    short mg_levels_;
+    short mg_npresmoothing_;
+    short mg_npostsmoothing_;
+
+
     // dielectric model for solvation
     short diel;
     // Parameters for MG solver/ preconditioner for Poisson problem
diff --git a/src/MGOrbitalsPreconditioning.cc b/src/MGOrbitalsPreconditioning.cc
index 957495d9..cdd6c44c 100644
--- a/src/MGOrbitalsPreconditioning.cc
+++ b/src/MGOrbitalsPreconditioning.cc
@@ -22,7 +22,15 @@
 template <class OrbitalsType, typename PDataType>
 MGOrbitalsPreconditioning<OrbitalsType, PDataType>::MGOrbitalsPreconditioning(
     const short mg_levels, const short lap_type)
-    : mg_levels_(mg_levels), lap_type_(lap_type), is_set_(false){};
+    : mg_levels_(mg_levels), lap_type_(lap_type), is_set_(false)
+{
+    Control& ct(*(Control::instance()));
+    Mesh* mymesh = Mesh::instance();
+    const pb::Grid& mygrid(mymesh->grid());
+
+    precond_ = std::make_shared<Preconditioning<PDataType>>(
+        lap_type_, mg_levels_, ct.mg_npresmoothing_, ct.mg_npostsmoothing_, mygrid, ct.bcWF);
+}
 
 template <class OrbitalsType, typename PDataType>
 MGOrbitalsPreconditioning<OrbitalsType, PDataType>::~MGOrbitalsPreconditioning()
@@ -42,9 +50,6 @@ void MGOrbitalsPreconditioning<OrbitalsType, PDataType>::setup(
     Mesh* mymesh = Mesh::instance();
     const pb::Grid& mygrid(mymesh->grid());
 
-    precond_ = std::make_shared<Preconditioning<PDataType>>(
-        lap_type_, mg_levels_, mygrid, ct.bcWF);
-
     if (currentMasks != nullptr)
     {
         // set masks in GridFuncVector class
diff --git a/src/PCGSolver.cc b/src/PCGSolver.cc
index 254b1a23..cb22ab24 100644
--- a/src/PCGSolver.cc
+++ b/src/PCGSolver.cc
@@ -202,7 +202,7 @@ bool PCGSolver<OperatorType, ScalarDataType, PrecondDataType>::solve(
     pb::GridFunc<PrecondDataType> prec_z(finegrid, bc_[0], bc_[1], bc_[2]);
     pb::GridFunc<PrecondDataType> prec_res(res);
     /* preconditioning step */
-    prec_z.setValues(0.);
+    prec_z.setZero();
     preconSolve(prec_z, prec_res, 0);
     pb::GridFunc<ScalarDataType> z(prec_z);
 
@@ -236,7 +236,7 @@ bool PCGSolver<OperatorType, ScalarDataType, PrecondDataType>::solve(
             converged = true;
             break;
         }
-        prec_z.setValues(0.);
+        prec_z.setZero();
         prec_res.setValues(res);
         preconSolve(prec_z, prec_res, 0);
         z.setValues(prec_z);
diff --git a/src/Preconditioning.cc b/src/Preconditioning.cc
index 495d9703..491b9974 100644
--- a/src/Preconditioning.cc
+++ b/src/Preconditioning.cc
@@ -10,13 +10,19 @@
 #include "Preconditioning.h"
 #include "LapFactory.h"
 
-using namespace std;
-
 template <typename T>
 Preconditioning<T>::Preconditioning(const short lap_type, const short maxlevels,
-    const pb::Grid& grid, const short bcWF[3])
+    const short npresmooth, const short npostsmooth,
+    const pb::Grid& grid, const short bcWF[3]):
+	max_levels_(maxlevels),
+	npresmooth_(npresmooth),
+	npostsmooth_(npostsmooth)
 {
-    max_levels_ = maxlevels;
+    assert(npresmooth_>=0);
+    assert(npostsmooth_>=0);
+    assert(npresmooth_<100);
+    assert(npostsmooth_<100);
+
     for (short i = 0; i < 3; i++)
         bc_[i] = bcWF[i];
 
@@ -27,16 +33,6 @@ Preconditioning<T>::Preconditioning(const short lap_type, const short maxlevels,
     jacobi_factor_.push_back(myoper->jacobiFactor());
 }
 
-template <typename T>
-Preconditioning<T>::Preconditioning(const Preconditioning& precond)
-{
-    max_levels_      = precond.max_levels_;
-    pb::Grid* mygrid = new pb::Grid(*(precond.grid_[0]));
-    grid_.push_back(mygrid);
-    for (short i = 0; i < 3; i++)
-        bc_[i] = precond.bc_[i];
-}
-
 template <typename T>
 Preconditioning<T>::~Preconditioning()
 {
@@ -150,8 +146,7 @@ void Preconditioning<T>::setup(
     }
 }
 
-// MG V-cycle with mask corresponding to state istate
-// (no mask if istate==-1)
+// MG V-cycle
 template <typename T>
 void Preconditioning<T>::mg(pb::GridFuncVector<T, memory_space_type>& gfv_v,
     const pb::GridFuncVector<T, memory_space_type>& gfv_f, const short lap_type,
@@ -164,8 +159,8 @@ void Preconditioning<T>::mg(pb::GridFuncVector<T, memory_space_type>& gfv_v,
     assert(static_cast<int>(gfv_work_.size()) > level);
     assert(gfv_work_[level] != nullptr);
 
-    short ncycl = 2;
-    if (level == max_levels_) ncycl = 4;
+    short ncycl = npresmooth_;
+    if (level == max_levels_) ncycl = npresmooth_+npostsmooth_;
 
     const double jacobi_factor = jacobi_factor_[level];
 
@@ -206,7 +201,7 @@ void Preconditioning<T>::mg(pb::GridFuncVector<T, memory_space_type>& gfv_v,
     gfv_v -= (*gfv_work_[level]);
 
     // post-smoothing
-    for (short it = 0; it < 2; it++)
+    for (short it = 0; it < npostsmooth_; it++)
     {
         gfv_v.jacobi(lap_type, gfv_f, *gfv_work_[level], jacobi_factor);
         gfv_v.app_mask(level);
diff --git a/src/Preconditioning.h b/src/Preconditioning.h
index 9f28748f..33aad1fb 100644
--- a/src/Preconditioning.h
+++ b/src/Preconditioning.h
@@ -30,8 +30,12 @@ class Preconditioning
 #endif
 
 private:
-    short max_levels_;
-    short lap_type_;
+
+    // V-cycle parameters
+    const short max_levels_;
+    const short npresmooth_;
+    const short npostsmooth_;
+
     short bc_[3];
 
     // Jacobi factor at each level
@@ -48,9 +52,10 @@ class Preconditioning
 
 public:
     Preconditioning(const short lap_type, const short maxlevels,
+        const short npresmooth, const short npostsmooth,
         const pb::Grid& grid, const short bc[3]);
 
-    Preconditioning(const Preconditioning&);
+    Preconditioning(const Preconditioning&) = delete;
 
     ~Preconditioning();
     void clear();
diff --git a/src/pb/GridFunc.cc b/src/pb/GridFunc.cc
index f883d0d7..d80bad01 100644
--- a/src/pb/GridFunc.cc
+++ b/src/pb/GridFunc.cc
@@ -78,9 +78,7 @@ double GridFunc<double>::fmax()
     int n    = grid_.sizeg();
     int imax = IDAMAX(&n, &uu_[0], &ione) - 1;
 
-    double vmax = mype_env().double_max_all(uu_[imax]);
-
-    return vmax;
+    return mype_env().double_max_all(uu_[imax]);
 }
 
 template <>
@@ -90,9 +88,7 @@ double GridFunc<float>::fmax()
     int n    = grid_.sizeg();
     int imax = ISAMAX(&n, &uu_[0], &ione) - 1;
 
-    double vmax = mype_env().double_max_all((double)uu_[imax]);
-
-    return vmax;
+    return mype_env().double_max_all((double)uu_[imax]);
 }
 
 template <typename T>
@@ -219,8 +215,7 @@ GridFunc<T>::GridFunc(const GridFunc<double>& A) : grid_(A.grid())
 
     alloc();
 
-    int n = grid_.sizeg();
-    MPcpy(uu_, A.uu(), n);
+    MPcpy(uu_, A.uu(), grid_.sizeg());
 
     updated_boundaries_ = A.updated_boundaries();
 }
@@ -276,13 +271,11 @@ GridFunc<T>::GridFunc(const GridFunc<T>& A, const Grid& new_grid)
 
     for (int ix = 0; ix < dim_[0]; ix++)
     {
-
         int ix1 = (ix + shift1) * incx_ + shift1;
         int ix2 = (ix + shift2) * incx2 + shift2;
 
         for (int iy = 0; iy < dim_[1]; iy++)
         {
-
             int iy1 = ix1 + (iy + shift1) * incy_;
             int iy2 = ix2 + (iy + shift2) * incy2;
 
@@ -331,13 +324,11 @@ void GridFunc<T>::set_max(const T val)
 }
 
 template <typename T>
-void GridFunc<T>::setValues(const int n, const T* src, const int pos)
+void GridFunc<T>::setValues(const int n, const T* src)
 {
-    assert((pos + n) <= static_cast<int>(grid_.sizeg()));
+    assert(n <= static_cast<int>(grid_.sizeg()));
 
-    size_t ssize = n * sizeof(T);
-    // int ione=1;
-    memcpy(&uu_[pos], src, ssize);
+    memcpy(uu_, src, n * sizeof(T));
 }
 
 template <typename T>
@@ -365,6 +356,14 @@ GridFunc<T>& GridFunc<T>::operator=(const GridFunc<T>& func)
     return *this;
 }
 
+template <typename T>
+void GridFunc<T>::setZero()
+{
+    memset(uu_, 0, grid_.sizeg()*sizeof(T));
+
+    updated_boundaries_ = true;
+}
+
 template <typename T>
 void GridFunc<T>::setValues(const T val)
 {
@@ -373,6 +372,8 @@ void GridFunc<T>::setValues(const T val)
 
     for (int i = 0; i < n; i++)
         pu[i] = val;
+
+    updated_boundaries_ = true;
 }
 
 template <typename T>
@@ -380,7 +381,6 @@ GridFunc<T>& GridFunc<T>::operator=(const T val)
 {
     setValues(val);
 
-    updated_boundaries_ = true;
     return *this;
 }
 
@@ -516,7 +516,6 @@ GridFunc<T>& GridFunc<T>::operator*=(const double alpha)
 template <typename T>
 void GridFunc<T>::scal(const double alpha)
 {
-
     LinearAlgebraUtils<MemorySpace::Host>::MPscal(grid_.sizeg(), alpha, uu_);
 }
 
diff --git a/src/pb/GridFunc.h b/src/pb/GridFunc.h
index ab72c86d..1e0930a3 100644
--- a/src/pb/GridFunc.h
+++ b/src/pb/GridFunc.h
@@ -95,7 +95,7 @@ class GridFunc : public GridFuncInterface
     static std::vector<T> buf3_;
     static std::vector<T> buf4_;
 
-    void setValues(const int n, const T* src, const int pos = 0);
+    void setValues(const int n, const T* src);
 
 public:
     // Constructors
@@ -114,6 +114,7 @@ class GridFunc : public GridFuncInterface
 
     void setValues(const GridFunc<T>& src);
     void setValues(const T val);
+    void setZero();
 
     int inc(const short dir) const { return grid_.inc(dir); }
 
diff --git a/src/read_config.cc b/src/read_config.cc
index e3272b3b..2444fd6a 100644
--- a/src/read_config.cc
+++ b/src/read_config.cc
@@ -118,10 +118,16 @@ int read_config(int argc, char** argv, po::variables_map& vm,
             "Compute MLWF (apply rotation) in quench")(
             "Quench.num_lin_iterations", po::value<short>()->default_value(0),
             "Number of iterations without potential update in quench")(
-            "Quench.preconditioner_num_levels",
+            "Preconditioner.num_levels",
             po::value<short>()->default_value(2),
             "Number of levels for MG preconditioner")(
-            "Quench.preconditioner_precision",
+            "Preconditioner.npresmoothing",
+            po::value<short>()->default_value(2),
+            "Number of presmoothing steps i preconditioner")(
+            "Preconditioner.npostsmoothing",
+            po::value<short>()->default_value(2),
+            "Number of postsmoothing steps i preconditioner")(
+            "Preconditioner.precision",
             po::value<short>()->default_value(32),
             "Precision for MG preconditioner")("Quench.spread_penalty_damping",
             po::value<float>()->default_value(0.),
diff --git a/tests/LBFGS/lbfgs1.cfg b/tests/LBFGS/lbfgs1.cfg
index c64cb712..0de05091 100644
--- a/tests/LBFGS/lbfgs1.cfg
+++ b/tests/LBFGS/lbfgs1.cfg
@@ -24,7 +24,8 @@ dt=8.
 [Quench]
 max_steps=50
 atol=1.e-8
-preconditioner_num_levels=1
+[Preconditioner]
+num_levels=1
 [Orbitals]
 initial_type=Gaussian
 initial_width=2.
diff --git a/tests/LBFGS/lbfgs2.cfg b/tests/LBFGS/lbfgs2.cfg
index 78a7e737..36064518 100644
--- a/tests/LBFGS/lbfgs2.cfg
+++ b/tests/LBFGS/lbfgs2.cfg
@@ -24,7 +24,8 @@ dt=8.
 [Quench]
 max_steps=50
 atol=1.e-8
-preconditioner_num_levels=1
+[Preconditioner]
+num_levels=1
 [Orbitals]
 initial_type=Gaussian
 initial_width=2.
diff --git a/tests/MLWF/mlwf.cfg b/tests/MLWF/mlwf.cfg
index cc5c52d1..d03cfad9 100644
--- a/tests/MLWF/mlwf.cfg
+++ b/tests/MLWF/mlwf.cfg
@@ -20,8 +20,9 @@ type=QUENCH
 max_steps=300
 atol=1.e-8
 MLWF=true
-preconditioner_num_levels=1
 step_length=1.5
+[Preconditioner]
+num_levels=1
 [Poisson]
 FDtype=4th
 [Orbitals]
diff --git a/tests/SiH4/mgmol.cfg b/tests/SiH4/mgmol.cfg
index aea52881..ebfecf93 100644
--- a/tests/SiH4/mgmol.cfg
+++ b/tests/SiH4/mgmol.cfg
@@ -20,7 +20,8 @@ type=QUENCH
 max_steps=45
 atol=1.e-9
 num_lin_iterations=2
-preconditioner_precision=64
+[Preconditioner]
+precision=64
 [Orbitals]
 initial_type=Gaussian
 initial_width=2.