Index: src/Jobs/InterfaceVMGJob.cpp
===================================================================
--- src/Jobs/InterfaceVMGJob.cpp	(revision 6ff62c8bb90891ea66dfc4882a5f7298de44ba0c)
+++ src/Jobs/InterfaceVMGJob.cpp	(revision 10f0cb73e0cc1105d7dc5360f9c6409f9674a51e)
@@ -71,4 +71,5 @@
     vmg_float _box_end,
     const int& near_field_cells,
+    const bool _DoPrintDebug,
     int coarseningSteps,
     double alpha) :
@@ -78,5 +79,6 @@
   sampled_input(_sampled_input),
   returndata(_returndata),
-  level(levelMax)
+  level(levelMax),
+  DoPrintDebug(_DoPrintDebug)
 {
   for (size_t i=0;i<3;++i) {
@@ -144,6 +146,8 @@
   comm.CommFromGhosts(particle_grid);
 
-  // print nuclei grid to vtk
-  comm.PrintGrid(particle_grid, "Sampled Nuclei Density");
+  if (DoPrintDebug) {
+    // print nuclei grid to vtk
+    comm.PrintGrid(particle_grid, "Sampled Nuclei Density");
+  }
 
   // add sampled electron charge density onto grid
@@ -153,6 +157,8 @@
       -1.);
 
-  // print electron grid to vtk
-  comm.PrintGrid(grid, "Sampled Electron Density");
+  if (DoPrintDebug) {
+    // print electron grid to vtk
+    comm.PrintGrid(grid, "Sampled Electron Density");
+  }
 
   // add particle_grid onto grid
@@ -180,6 +186,8 @@
   comm.PrintStringOnce("Grid charge integral: %e", charge_sum/(4.0 * VMG::Math::pi));
 
-  // print total grid to vtk
-  comm.PrintGrid(grid, "Total Charge Density");
+  if (DoPrintDebug) {
+    // print total grid to vtk
+    comm.PrintGrid(grid, "Total Charge Density");
+  }
 
 //  delete temp_grid;
@@ -197,6 +205,8 @@
   Index i;
 
-  // print output grid to vtk
-  comm.PrintGrid(grid, "Potential Solution");
+  if (DoPrintDebug) {
+    // print output grid to vtk
+    comm.PrintGrid(grid, "Potential Solution");
+  }
 
   // obtain sampled potential from grid
Index: src/Jobs/InterfaceVMGJob.hpp
===================================================================
--- src/Jobs/InterfaceVMGJob.hpp	(revision 6ff62c8bb90891ea66dfc4882a5f7298de44ba0c)
+++ src/Jobs/InterfaceVMGJob.hpp	(revision 10f0cb73e0cc1105d7dc5360f9c6409f9674a51e)
@@ -46,6 +46,8 @@
       vmg_float _box_end,
       const int& near_field_cells,
+      const bool _DoPrintDebug=false,
       int coarseningSteps=0,
-      double alpha=1.6);
+      double alpha=1.6
+      );
 
   ~InterfaceVMGJob() {}
@@ -71,4 +73,7 @@
   //!> contains end of grid intervals
   double box_end[3];
+
+  //!> whether we do print grid for debug visualization or not
+  const bool DoPrintDebug;
 };
 
Index: src/Jobs/VMGJob.cpp
===================================================================
--- src/Jobs/VMGJob.cpp	(revision 6ff62c8bb90891ea66dfc4882a5f7298de44ba0c)
+++ src/Jobs/VMGJob.cpp	(revision 10f0cb73e0cc1105d7dc5360f9c6409f9674a51e)
@@ -82,5 +82,6 @@
     const std::vector< double > &_particle_charges,
     const size_t _near_field_cells,
-    const size_t _interpolation_degree) :
+    const size_t _interpolation_degree,
+    const bool _DoPrintDebug) :
   FragmentJob(_JobId),
   density_grid(_density_grid),
@@ -90,5 +91,6 @@
   interpolation_degree(_interpolation_degree),
   returndata(static_cast<const SamplingGridProperties &>(_density_grid)),
-  particles()
+  particles(),
+  DoPrintDebug(_DoPrintDebug)
 {}
 
@@ -97,5 +99,6 @@
   near_field_cells(3),
   interpolation_degree(3),
-  particles()
+  particles(),
+  DoPrintDebug(false)
 {}
 
@@ -216,5 +219,6 @@
       Vector(density_grid.begin),
       density_grid.end[0]-density_grid.begin[0],
-      near_field_cells);
+      near_field_cells,
+      DoPrintDebug);
   new LevelOperatorCS(Stencils::RestrictionFullWeight, Stencils::InterpolationTrilinear);
   new Givens<SolverSingular>();
Index: src/Jobs/VMGJob.hpp
===================================================================
--- src/Jobs/VMGJob.hpp	(revision 6ff62c8bb90891ea66dfc4882a5f7298de44ba0c)
+++ src/Jobs/VMGJob.hpp	(revision 10f0cb73e0cc1105d7dc5360f9c6409f9674a51e)
@@ -42,4 +42,5 @@
    * @param _interpolation_degree degree of interpolation polynomial for getting nuclei
    *        potential from grid
+   * @param _DoPrintDebug whether we do print grid for debug visualization or not
    */
   VMGJob(const JobId_t _JobId,
@@ -48,5 +49,6 @@
       const std::vector< double > &_particle_charges,
       const size_t _near_field_cells,
-      const size_t _interpolation_degree
+      const size_t _interpolation_degree,
+      const bool _DoPrintDebug=false
       );
   virtual ~VMGJob();
@@ -70,4 +72,6 @@
   //!> interpolation degree used in sampling the potential of the nuclei
   const size_t interpolation_degree;
+  //!> whether we do print grid for debug visualization or not
+  const bool DoPrintDebug;
 
 private:
@@ -120,4 +124,5 @@
     ar & const_cast< size_t &>(near_field_cells);
     ar & const_cast< size_t &>(interpolation_degree);
+    ar & const_cast< bool &>(DoPrintDebug);
     ar & returndata;
   }
