diff --git a/src/Control.cc b/src/Control.cc
index 2231813e..f1dcef23 100644
--- a/src/Control.cc
+++ b/src/Control.cc
@@ -70,6 +70,7 @@ Control::Control()
     dm_approx_ndigits      = 1;
     dm_approx_power_maxits = 100;
     wf_extrapolation_      = 1;
+    verbose                = 0;
 
     // undefined values
     dm_algo_                         = -1;
diff --git a/src/DensityMatrix.cc b/src/DensityMatrix.cc
index 6e6ea652..e2c7b473 100644
--- a/src/DensityMatrix.cc
+++ b/src/DensityMatrix.cc
@@ -28,7 +28,7 @@ const double factor_kernel4dot = 10.;
 template <class MatrixType>
 DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
 {
-    assert(ndim >= 0);
+    assert(ndim > 0);
 
     dim_ = ndim;
 
@@ -45,6 +45,7 @@ DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
     kernel4dot_ = new MatrixType("K4dot", ndim, ndim);
     work_       = new MatrixType("work", ndim, ndim);
     occupation_.resize(dim_);
+    setDummyOcc();
 }
 
 template <class MatrixType>
@@ -109,6 +110,7 @@ void DensityMatrix<MatrixType>::build(
     const std::vector<double>& occ, const int new_orbitals_index)
 {
     assert(dm_ != nullptr);
+    assert(!occ.empty());
 
     setOccupations(occ);
 
@@ -149,6 +151,8 @@ template <class MatrixType>
 void DensityMatrix<MatrixType>::setUniform(
     const double nel, const int new_orbitals_index)
 {
+    assert(!occupation_.empty());
+
     MGmol_MPI& mmpi  = *(MGmol_MPI::instance());
     const double occ = (double)((double)nel / (double)dim_);
     if (mmpi.instancePE0())
@@ -314,6 +318,7 @@ void DensityMatrix<MatrixType>::computeOccupations(const MatrixType& ls)
 template <class MatrixType>
 void DensityMatrix<MatrixType>::setOccupations(const std::vector<double>& occ)
 {
+    assert(!occ.empty());
 #ifdef PRINT_OPERATIONS
     MGmol_MPI& mmpi = *(MGmol_MPI::instance());
     if (mmpi.instancePE0())
diff --git a/src/DensityMatrix.h b/src/DensityMatrix.h
index 960ad18e..84804a86 100644
--- a/src/DensityMatrix.h
+++ b/src/DensityMatrix.h
@@ -84,17 +84,24 @@ class DensityMatrix
         *dm_            = mat;
         orbitals_index_ = orbitals_index;
 
-        occupation_.clear();
+        setDummyOcc();
 
         occ_uptodate_ = false;
         uniform_occ_  = false;
         stripped_     = false;
     }
 
+    // set occupations to meaningless values to catch uninitialized use
+    void setDummyOcc()
+    {
+        for (auto& occ : occupation_)
+            occ = -1.;
+    }
+
     void initMatrix(const double* const val)
     {
         dm_->init(val, dim_);
-        occupation_.clear();
+        setDummyOcc();
 
         occ_uptodate_ = false;
         uniform_occ_  = false;
@@ -105,8 +112,10 @@ class DensityMatrix
 
     void getOccupations(std::vector<double>& occ) const
     {
+        assert(!occupation_.empty());
         assert(occ_uptodate_);
         assert((int)occ.size() == dim_);
+
         memcpy(&occ[0], &occupation_[0], dim_ * sizeof(double));
     }
 
diff --git a/src/Electrostatic.cc b/src/Electrostatic.cc
index 80cb754b..6ac4fef3 100644
--- a/src/Electrostatic.cc
+++ b/src/Electrostatic.cc
@@ -11,6 +11,7 @@
 #include "Control.h"
 #include "ExtendedGridOrbitals.h"
 #include "GridFactory.h"
+#include "GridFunc.h"
 #include "Hartree.h"
 #include "Hartree_CG.h"
 #include "Ions.h"
@@ -23,15 +24,6 @@
 #include "ShiftedHartree.h"
 #include "mputils.h"
 
-#include "GridFunc.h"
-#include "Laph2.h"
-#include "Laph4.h"
-#include "Laph4M.h"
-#include "Laph4MP.h"
-#include "Laph6.h"
-#include "Laph8.h"
-#include "ShiftedLaph4M.h"
-
 Timer Electrostatic::solve_tm_("Electrostatic::solve");
 
 Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3],
@@ -49,109 +41,9 @@ Electrostatic::Electrostatic(PoissonFDtype lap_type, const short bcPoisson[3],
     Mesh* mymesh           = Mesh::instance();
     const pb::Grid& myGrid = mymesh->grid();
 
-    Control& ct = *(Control::instance());
-    if (ct.MGPoissonSolver()) // use MG for Poisson Solver
-    {
-        if (screening_const > 0.)
-        {
-            switch (lap_type)
-            {
-                case PoissonFDtype::h4M:
-                    poisson_solver_
-                        = new ShiftedHartree<pb::ShiftedLaph4M<POTDTYPE>>(
-                            myGrid, bc_, screening_const);
-                    break;
-                default:
-                    (*MPIdata::sout)
-                        << "Electrostatic, shifted, Undefined option: "
-                        << static_cast<int>(lap_type) << std::endl;
-            }
-        }
-        else
-        {
-            switch (lap_type)
-            {
-                case PoissonFDtype::h4M:
-                    poisson_solver_
-                        = new Hartree<pb::Laph4M<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h2:
-                    poisson_solver_
-                        = new Hartree<pb::Laph2<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h4:
-                    poisson_solver_
-                        = new Hartree<pb::Laph4<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h6:
-                    poisson_solver_
-                        = new Hartree<pb::Laph6<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h8:
-                    poisson_solver_
-                        = new Hartree<pb::Laph8<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h4MP:
-                    poisson_solver_
-                        = new Hartree<pb::Laph4MP<POTDTYPE>>(myGrid, bc_);
-                    break;
-                default:
-                    (*MPIdata::sout) << "Electrostatic, Undefined option: "
-                                     << static_cast<int>(lap_type) << std::endl;
-            }
-        }
-    }
-    else // use PCG for Poisson Solver
-    {
-        if (screening_const > 0.)
-        {
-            switch (lap_type)
-            {
-                case PoissonFDtype::h4M:
-                    poisson_solver_
-                        = new ShiftedHartree<pb::ShiftedLaph4M<POTDTYPE>>(
-                            myGrid, bc_, screening_const);
-                    break;
-                default:
-                    (*MPIdata::sout)
-                        << "PCG Electrostatic, shifted, Undefined option: "
-                        << static_cast<int>(lap_type) << std::endl;
-            }
-        }
-        else
-        {
-            switch (lap_type)
-            {
-                case PoissonFDtype::h4M:
-                    poisson_solver_
-                        = new Hartree_CG<pb::Laph4M<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h2:
-                    poisson_solver_
-                        = new Hartree_CG<pb::Laph2<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h4:
-                    poisson_solver_
-                        = new Hartree_CG<pb::Laph4<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h6:
-                    poisson_solver_
-                        = new Hartree_CG<pb::Laph6<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h8:
-                    poisson_solver_
-                        = new Hartree_CG<pb::Laph8<POTDTYPE>>(myGrid, bc_);
-                    break;
-                case PoissonFDtype::h4MP:
-                    poisson_solver_
-                        = new Hartree_CG<pb::Laph4MP<POTDTYPE>>(myGrid, bc_);
-                    break;
-                default:
-                    (*MPIdata::sout) << "PCG Electrostatic, Undefined option: "
-                                     << static_cast<int>(lap_type) << std::endl;
-            }
-        }
-    }
+    // create Poisson solver
+    poisson_solver_ = PoissonSolverFactory::create(
+        myGrid, lap_type, bcPoisson, screening_const);
 
     grhoc_     = nullptr;
     diel_flag_ = false;
@@ -244,73 +136,8 @@ void Electrostatic::setupPB(
         ngpts, origin, cell, static_cast<int>(laptype_), true, myPEenv);
     if (poisson_solver_ != nullptr) delete poisson_solver_;
 
-    Control& ct = *(Control::instance());
-    if (ct.MGPoissonSolver()) // use MG for Poisson Solver
-    {
-        switch (laptype_)
-        {
-            case PoissonFDtype::h4M:
-                poisson_solver_ = new PBdiel<pb::PBh4M<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h2:
-                poisson_solver_ = new PBdiel<pb::PBh2<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h4:
-                poisson_solver_ = new PBdiel<pb::PBh4<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h6:
-                poisson_solver_ = new PBdiel<pb::PBh6<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h8:
-                poisson_solver_ = new PBdiel<pb::PBh8<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h4MP:
-                poisson_solver_ = new PBdiel<pb::PBh4MP<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            default:
-                (*MPIdata::sout)
-                    << "Electrostatic, Undefined option" << std::endl;
-        }
-    }
-    else // use PCG for Poisson Solver
-    {
-        switch (laptype_)
-        {
-            case PoissonFDtype::h4M:
-                poisson_solver_ = new PBdiel_CG<pb::PBh4M<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h2:
-                poisson_solver_ = new PBdiel_CG<pb::PBh2<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h4:
-                poisson_solver_ = new PBdiel_CG<pb::PBh4<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h6:
-                poisson_solver_ = new PBdiel_CG<pb::PBh6<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h8:
-                poisson_solver_ = new PBdiel_CG<pb::PBh8<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            case PoissonFDtype::h4MP:
-                poisson_solver_ = new PBdiel_CG<pb::PBh4MP<POTDTYPE>>(
-                    *pbGrid_, bc_, e0, rho0, drho0);
-                break;
-            default:
-                (*MPIdata::sout)
-                    << "Electrostatic, Undefined option" << std::endl;
-        }
-    }
+    poisson_solver_ = PoissonSolverFactory::createDiel(
+        *pbGrid_, laptype_, bc_, e0, rho0, drho0);
 
     if (grhoc_ != nullptr)
     {
@@ -330,6 +157,7 @@ void Electrostatic::setupPB(
     poisson_solver_->set_vh(gf_vh);
 }
 
+// This function is only useful for Hartree problem with dielectric continuum
 void Electrostatic::fillFuncAroundIons(const Ions& ions)
 {
     assert(grhod_ != nullptr);
@@ -352,7 +180,6 @@ void Electrostatic::fillFuncAroundIons(const Ions& ions)
     std::vector<Ion*>::const_iterator ion = rc_ions.begin();
     while (ion != rc_ions.end())
     {
-
         double rc = (*ion)->getRC();
         // Special case: silicon
         if ((*ion)->isMass28()) rc = 2.0;
@@ -373,43 +200,32 @@ void Electrostatic::fillFuncAroundIons(const Ions& ions)
 #endif
             for (unsigned int ix = 0; ix < pbGrid_->dim(0); ix++)
             {
-
                 xc[1]         = pbGrid_->start(1);
                 const int ix1 = (ix + shift) * incx;
 
                 for (unsigned int iy = 0; iy < pbGrid_->dim(1); iy++)
                 {
-
                     xc[2]         = pbGrid_->start(2);
                     const int iy1 = ix1 + (iy + shift) * incy;
 
                     for (unsigned int iz = 0; iz < pbGrid_->dim(2); iz++)
                     {
-
                         const double r = (*ion)->minimage(xc, lattice, bc_);
 
                         if (r < rc)
                         {
                             const double alpha = 0.2 * (1. + cos(r * pi_rc));
-
-                            const int iz1 = iy1 + iz + shift;
+                            const int iz1      = iy1 + iz + shift;
                             vv[iz1] += alpha;
                         }
-
                         xc[2] += pbGrid_->hgrid(2);
                     }
-
                     xc[1] += pbGrid_->hgrid(1);
-
                 } // end for iy
-
                 xc[0] += pbGrid_->hgrid(0);
-
             } // end for ix
         }
-
         ion++;
-
     } // end loop on list of ions
 
     return;
diff --git a/src/Electrostatic.h b/src/Electrostatic.h
index 785bc909..7ae977e1 100644
--- a/src/Electrostatic.h
+++ b/src/Electrostatic.h
@@ -13,6 +13,7 @@
 #include "Control.h"
 #include "GridFunc.h"
 #include "Poisson.h"
+#include "PoissonSolverFactory.h"
 #include "Rho.h"
 #include "Timer.h"
 
diff --git a/src/Hartree.cc b/src/Hartree.cc
index b3984b87..7d0ba4d8 100644
--- a/src/Hartree.cc
+++ b/src/Hartree.cc
@@ -7,10 +7,8 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-// $Id$
 #include <iomanip>
 #include <iostream>
-using namespace std;
 
 #include "Control.h"
 #include "Hartree.h"
@@ -82,6 +80,8 @@ void Hartree<T>::solve(
     //    {
     /* solve with POTDTYPE precision */
     pb::GridFunc<POTDTYPE> rhs(work_rho);
+    // Hartree units
+    rhs *= (4. * M_PI);
     poisson_solver_->solve(*Poisson::vh_, rhs);
     //    }
     //    else
@@ -92,10 +92,11 @@ void Hartree<T>::solve(
     double final_residual     = poisson_solver_->getFinalResidual();
 
     if (onpe0)
-        (*MPIdata::sout) << setprecision(2) << scientific
+        (*MPIdata::sout) << std::setprecision(2) << std::scientific
                          << "Hartree: residual reduction = "
                          << residual_reduction
-                         << ", final residual = " << final_residual << endl;
+                         << ", final residual = " << final_residual
+                         << std::endl;
 
     Poisson::Int_vhrho_  = vel * Poisson::vh_->gdot(rho);
     Poisson::Int_vhrhoc_ = vel * Poisson::vh_->gdot(rhoc);
diff --git a/src/Hartree.h b/src/Hartree.h
index a7050e24..b87ff01e 100644
--- a/src/Hartree.h
+++ b/src/Hartree.h
@@ -7,14 +7,12 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-// $Id$
 #ifndef included_Hartree
 #define included_Hartree
 
 #include "Poisson.h"
 #include "PoissonInterface.h"
 
-// pb
 #include "SolverLap.h"
 
 template <class T>
diff --git a/src/Hartree_CG.cc b/src/Hartree_CG.cc
index f9bf3860..375237f8 100644
--- a/src/Hartree_CG.cc
+++ b/src/Hartree_CG.cc
@@ -81,6 +81,7 @@ void Hartree_CG<T>::solve(
     //    {
     /* solve with POTDTYPE precision */
     pb::GridFunc<POTDTYPE> rhs(work_rho);
+    rhs *= (4. * M_PI);
     poisson_solver_->solve(*Poisson::vh_, rhs);
     //    }
     //    else
diff --git a/src/Ions.cc b/src/Ions.cc
index 940e12e8..852151c7 100644
--- a/src/Ions.cc
+++ b/src/Ions.cc
@@ -584,11 +584,9 @@ void Ions::printPositionsLocal(std::ostream& os, const int root) const
         os.setf(std::ios::right, std::ios::adjustfield);
         os.setf(std::ios::fixed, std::ios::floatfield);
 
-        std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-        while (ion != local_ions_.end())
+        for (auto& ion : local_ions_)
         {
-            (*ion)->printPosition(os);
-            ion++;
+            ion->printPosition(os);
         }
 
         os << std::endl;
@@ -627,14 +625,12 @@ void Ions::writeAtomicNumbers(HDFrestart& h5f_file)
     }
     else
     {
-        std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-        while (ion != local_ions_.end())
+        for (auto& ion : local_ions_)
         {
-            assert((*ion)->atomic_number() > 0);
-            assert((*ion)->atomic_number() < 200);
+            assert(ion->atomic_number() > 0);
+            assert(ion->atomic_number() < 200);
 
-            data.push_back((*ion)->atomic_number());
-            ion++;
+            data.push_back(ion->atomic_number());
         }
     }
 
@@ -685,17 +681,15 @@ void Ions::writeAtomNames(HDFrestart& h5f_file)
 
 void Ions::lockAtom(const std::string& name)
 {
-    std::vector<Ion*>::iterator ion = local_ions_.begin();
-    while (ion != local_ions_.end())
+    for (auto& ion : local_ions_)
     {
-        std::string name_ion((*ion)->name());
+        std::string name_ion(ion->name());
         if (name.compare(name_ion) == 0)
         {
-            (*ion)->lock();
+            ion->lock();
             if (onpe0) (*MPIdata::sout) << "Lock atom " << name << std::endl;
             break;
         }
-        ion++;
     }
 }
 
@@ -975,12 +969,10 @@ void Ions::readRestartPositions(HDFrestart& h5_file)
     std::vector<double> data;
     h5_file.readAtomicPositions(data);
 
-    int i                           = 0;
-    std::vector<Ion*>::iterator ion = local_ions_.begin();
-    while (ion != local_ions_.end())
+    int i = 0;
+    for (auto& ion : local_ions_)
     {
-        (*ion)->setPosition(data[3 * i], data[3 * i + 1], data[3 * i + 2]);
-        ion++;
+        ion->setPosition(data[3 * i], data[3 * i + 1], data[3 * i + 2]);
         i++;
     }
 }
@@ -1109,14 +1101,12 @@ void Ions::removeMassCenterMotion()
 #endif
     }
 
-    ion = local_ions_.begin();
-    i   = 0;
-    while (ion != local_ions_.end())
+    i = 0;
+    for (auto& ion : local_ions_)
     {
         const int threei = 3 * i;
-        (*ion)->setVelocity(velocities[threei] - tmp[0],
+        ion->setVelocity(velocities[threei] - tmp[0],
             velocities[threei + 1] - tmp[1], velocities[threei + 2] - tmp[2]);
-        ion++;
         i++;
     }
 
@@ -1124,16 +1114,15 @@ void Ions::removeMassCenterMotion()
     mv[0] = 0.;
     mv[1] = 0.;
     mv[2] = 0.;
-    while (ion != local_ions_.end())
+    for (auto& ion : local_ions_)
     {
-        velocities[3 * i]     = (*ion)->velocity(0);
-        velocities[3 * i + 1] = (*ion)->velocity(1);
-        velocities[3 * i + 2] = (*ion)->velocity(2);
+        velocities[3 * i]     = ion->velocity(0);
+        velocities[3 * i + 1] = ion->velocity(1);
+        velocities[3 * i + 2] = ion->velocity(2);
         mv[0] += mass[i] * velocities[3 * i];
         mv[1] += mass[i] * velocities[3 * i + 1];
         mv[2] += mass[i] * velocities[3 * i + 2];
 
-        ion++;
         i++;
     }
 
@@ -1155,12 +1144,10 @@ void Ions::readRestartVelocities(HDFrestart& h5_file)
     std::vector<double> data;
     h5_file.readAtomicVelocities(data);
 
-    std::vector<Ion*>::iterator ion = local_ions_.begin();
-    int i                           = 0;
-    while (ion != local_ions_.end())
+    int i = 0;
+    for (auto& ion : local_ions_)
     {
-        (*ion)->setVelocity(data[3 * i], data[3 * i + 1], data[3 * i + 2]);
-        ion++;
+        ion->setVelocity(data[3 * i], data[3 * i + 1], data[3 * i + 2]);
         i++;
     }
 }
@@ -1174,12 +1161,10 @@ void Ions::readRestartRandomStates(HDFrestart& h5f_file)
     std::vector<unsigned short> data;
     h5f_file.readRestartRandomStates(data);
 
-    std::vector<Ion*>::iterator ion = local_ions_.begin();
-    int i                           = 0;
-    while (ion != local_ions_.end())
+    int i = 0;
+    for (auto& ion : local_ions_)
     {
-        (*ion)->setRandomState(data[3 * i], data[3 * i + 1], data[3 * i + 2]);
-        ion++;
+        ion->setRandomState(data[3 * i], data[3 * i + 1], data[3 * i + 2]);
         i++;
     }
 }
@@ -1202,17 +1187,14 @@ void Ions::writeForces(HDFrestart& h5f_file)
     }
     else
     {
-        std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-        while (ion != local_ions().end())
+        for (auto& ion : local_ions_)
         {
             // get position of local ion
             double force[3];
-            (*ion)->getForce(&force[0]);
+            ion->getForce(&force[0]);
             data.push_back(force[0]);
             data.push_back(force[1]);
             data.push_back(force[2]);
-
-            ++ion;
         }
     }
 
@@ -1364,64 +1346,54 @@ void Ions::printForcesLocal(std::ostream& os, const int root) const
            << "FX" << std::setw(10) << "FY" << std::setw(10) << "FZ"
            << std::endl;
 
-        std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-        while (ion != local_ions_.end())
+        for (auto& ion : local_ions_)
         {
+            ion->printPositionAndForce(os);
 
-            (*ion)->printPositionAndForce(os);
-
-            if (!(*ion)->locked())
+            if (!ion->locked())
             {
+                avg_forces[0] += fabs(ion->force(0));
+                avg_forces[1] += fabs(ion->force(1));
+                avg_forces[2] += fabs(ion->force(2));
 
-                avg_forces[0] += fabs((*ion)->force(0));
-                avg_forces[1] += fabs((*ion)->force(1));
-                avg_forces[2] += fabs((*ion)->force(2));
-
-                double ff = (*ion)->norm2F();
+                double ff = ion->norm2F();
                 maxf[0]   = std::max(maxf[0], ff);
 
-                max_forces[0] = std::max(max_forces[0], fabs((*ion)->force(0)));
-                max_forces[1] = std::max(max_forces[1], fabs((*ion)->force(1)));
-                max_forces[2] = std::max(max_forces[2], fabs((*ion)->force(2)));
+                max_forces[0] = std::max(max_forces[0], fabs(ion->force(0)));
+                max_forces[1] = std::max(max_forces[1], fabs(ion->force(1)));
+                max_forces[2] = std::max(max_forces[2], fabs(ion->force(2)));
 
                 num_movable++;
             }
 
             for (short ii = 0; ii < 3; ii++)
-                sum_forces[ii] += (*ion)->force(ii);
+                sum_forces[ii] += ion->force(ii);
             num_atoms++;
-
-            ion++;
         }
     }
     else
     {
-        std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-        while (ion != local_ions_.end())
+        for (auto& ion : local_ions_)
         {
-
-            if (!(*ion)->locked())
+            if (!ion->locked())
             {
+                avg_forces[0] += fabs(ion->force(0));
+                avg_forces[1] += fabs(ion->force(1));
+                avg_forces[2] += fabs(ion->force(2));
 
-                avg_forces[0] += fabs((*ion)->force(0));
-                avg_forces[1] += fabs((*ion)->force(1));
-                avg_forces[2] += fabs((*ion)->force(2));
-
-                double ff = (*ion)->norm2F();
+                double ff = ion->norm2F();
                 maxf[0]   = std::max(maxf[0], ff);
 
-                max_forces[0] = std::max(max_forces[0], fabs((*ion)->force(0)));
-                max_forces[1] = std::max(max_forces[1], fabs((*ion)->force(1)));
-                max_forces[2] = std::max(max_forces[2], fabs((*ion)->force(2)));
+                max_forces[0] = std::max(max_forces[0], fabs(ion->force(0)));
+                max_forces[1] = std::max(max_forces[1], fabs(ion->force(1)));
+                max_forces[2] = std::max(max_forces[2], fabs(ion->force(2)));
 
                 num_movable++;
             }
 
             for (short ii = 0; ii < 3; ii++)
-                sum_forces[ii] += (*ion)->force(ii);
+                sum_forces[ii] += ion->force(ii);
             num_atoms++;
-
-            ion++;
         }
     }
     // global statistics
@@ -1479,12 +1451,10 @@ int Ions::countIonsHere() const { return (int)local_ions_.size(); }
 
 int Ions::countProjectorsHere() const
 {
-    int count                             = 0;
-    std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-    while (ion != local_ions_.end())
+    int count = 0;
+    for (auto& ion : local_ions_)
     {
-        count += (*ion)->nProjectors();
-        ion++;
+        count += ion->nProjectors();
     }
     return count;
 }
@@ -1497,12 +1467,10 @@ int Ions::countProjectors() const
     Mesh* mymesh             = Mesh::instance();
     const pb::PEenv& myPEenv = mymesh->peenv();
 
-    int nproj                             = 0;
-    std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-    while (ion != local_ions_.end())
+    int nproj = 0;
+    for (auto& ion : local_ions_)
     {
-        nproj += (*ion)->nProjectors();
-        ion++;
+        nproj += ion->nProjectors();
     }
     int tmp = nproj;
     MPI_Allreduce(&tmp, &nproj, 1, MPI_INT, MPI_SUM, myPEenv.comm());
@@ -1569,13 +1537,10 @@ void Ions::setLocalPositions(const std::vector<double>& tau)
 {
     assert(tau.size() == 3 * local_ions_.size());
 
-    std::vector<Ion*>::iterator ion = local_ions_.begin();
-    int ia                          = 0;
-    while (ion != local_ions_.end())
+    int ia = 0;
+    for (auto& ion : local_ions_)
     {
-        (*ion)->setPosition(tau[3 * ia + 0], tau[3 * ia + 1], tau[3 * ia + 2]);
-
-        ion++;
+        ion->setPosition(tau[3 * ia + 0], tau[3 * ia + 1], tau[3 * ia + 2]);
         ia++;
     }
 
@@ -1780,45 +1745,8 @@ int Ions::setAtoms(
         if (ia < 1000) aname.append("0");
         if (ia < 10000) aname.append("0");
         aname.append(ss.str());
-        Ion* new_ion
-            = new Ion(species_[isp], aname, &crds[3 * ia], velocity, locked);
-        new_ion->bcast(mmpi.commGlobal());
-
-        // Populate list_ions_ list
-        // std::cout<<"crds: "<<crds[3*ia+0]<<", "<<crds[3*ia+1]<<",
-        // "<<crds[3*ia+2]<<std::endl;
-        if (inListIons(crds[3 * ia + 0], crds[3 * ia + 1], crds[3 * ia + 2]))
-        {
-            list_ions_.push_back(new_ion);
-            if (ct.verbose > 2)
-                (*MPIdata::sout)
-                    << "Ion " << aname << " at position " << crds[3 * ia + 0]
-                    << "," << crds[3 * ia + 1] << "," << crds[3 * ia + 2]
-                    << " added to the list... on PE" << mmpi.mypeGlobal()
-                    << std::endl;
-            // populate local_ions_ list
-            if (inLocalIons(
-                    crds[3 * ia + 0], crds[3 * ia + 1], crds[3 * ia + 2]))
-            {
-                (new_ion)->set_here(true);
-                local_ions_.push_back(new_ion);
-                if (onpe0 && ct.verbose > 2)
-                    (*MPIdata::sout)
-                        << "Ion " << aname << " at position "
-                        << crds[3 * ia + 0] << "," << crds[3 * ia + 1] << ","
-                        << crds[3 * ia + 2]
-                        << " added to the list of local ions... on PE"
-                        << mmpi.mypeGlobal() << std::endl;
-            }
-            else
-                (new_ion)->set_here(false);
-        }
-        else
-        {
-            //(*MPIdata::sout)<<"Ion "<<aname<<" NOT added to the list... on
-            // PE"<<mmpi.mype() <<endl;
-            delete new_ion;
-        }
+
+        addIonToList(species_[isp], aname, &crds[3 * ia], velocity, locked);
     }
     //    std::cout<<mmpi.mype()<<"...list size = "<<list_ions_.size()<<" local
     //    ions size = "<<local_ions_.size()<<endl;
@@ -1826,6 +1754,46 @@ int Ions::setAtoms(
     return natoms;
 }
 
+void Ions::addIonToList(const Species& sp, const std::string& name,
+    const double crds[3], const double velocity[3], const bool locked)
+{
+    MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
+    Control& ct(*(Control::instance()));
+
+    // create a new Ion
+    Ion* new_ion = new Ion(sp, name, &crds[0], velocity, locked);
+    new_ion->bcast(mmpi.commGlobal());
+
+    if (inListIons(crds[0], crds[1], crds[2]))
+    {
+        list_ions_.push_back(new_ion);
+        if (ct.verbose > 2)
+            (*MPIdata::sout)
+                << "Ion " << name << " at position " << crds[0] << ","
+                << crds[1] << "," << crds[2] << " added to the list... on PE"
+                << mmpi.mypeGlobal() << std::endl;
+
+        // populate local_ions_ list
+        if (inLocalIons(crds[0], crds[1], crds[2]))
+        {
+            (new_ion)->set_here(true);
+            local_ions_.push_back(new_ion);
+            if (onpe0 && ct.verbose > 2)
+                (*MPIdata::sout) << "Ion " << name << " at position " << crds[0]
+                                 << "," << crds[1] << "," << crds[2]
+                                 << " added to the list of local ions... on PE"
+                                 << mmpi.mypeGlobal() << std::endl;
+        }
+        else
+            (new_ion)->set_here(false);
+    }
+    else
+    {
+        // delete Ion if not put in list
+        delete new_ion;
+    }
+}
+
 int Ions::readNatoms(const std::string& filename, const bool cell_relative)
 {
     Control& ct(*(Control::instance()));
@@ -2117,9 +2085,8 @@ void Ions::setVelocities(const std::vector<double>& tau0,
     assert(tau0.size() == 3 * local_ions_.size());
     assert(taup.size() == 3 * local_ions_.size());
 
-    int ia                           = 0;
-    std::vector<Ion*>::iterator iion = local_ions_.begin();
-    while (iion != local_ions_.end())
+    int ia = 0;
+    for (auto& ion : local_ions_)
     {
         double v[3];
         for (short i = 0; i < 3; i++)
@@ -2128,8 +2095,7 @@ void Ions::setVelocities(const std::vector<double>& tau0,
             v[i] -= tau0[3 * ia + i];
             v[i] /= dt;
         }
-        (*iion)->setVelocity(v[0], v[1], v[2]);
-        iion++;
+        ion->setVelocity(v[0], v[1], v[2]);
         ia++;
     }
 }
@@ -2138,12 +2104,10 @@ void Ions::getLocalPositions(std::vector<double>& tau) const
 {
     assert(tau.size() == 3 * local_ions_.size());
 
-    int ia                                 = 0;
-    std::vector<Ion*>::const_iterator iion = local_ions_.begin();
-    while (iion != local_ions_.end())
+    int ia = 0;
+    for (auto& ion : local_ions_)
     {
-        (*iion)->getPosition(&tau[3 * ia]);
-        iion++;
+        ion->getPosition(&tau[3 * ia]);
         ia++;
     }
 }
@@ -2188,12 +2152,10 @@ void Ions::setTau0()
 {
     assert(tau0_.size() == 3 * local_ions_.size());
 
-    int ia                                 = 0;
-    std::vector<Ion*>::const_iterator iion = local_ions_.begin();
-    while (iion != local_ions_.end())
+    int ia = 0;
+    for (auto& ion : local_ions_)
     {
-        (*iion)->getPosition(&tau0_[3 * ia]);
-        iion++;
+        ion->getPosition(&tau0_[3 * ia]);
         ia++;
     }
 }
@@ -2202,13 +2164,11 @@ void Ions::setPositionsToTau0()
 {
     assert(tau0_.size() == 3 * local_ions_.size());
 
-    int ia                                 = 0;
-    std::vector<Ion*>::const_iterator iion = local_ions_.begin();
-    while (iion != local_ions_.end())
+    int ia = 0;
+    for (auto& ion : local_ions_)
     {
-        (*iion)->setPosition(
+        ion->setPosition(
             tau0_[3 * ia + 0], tau0_[3 * ia + 1], tau0_[3 * ia + 2]);
-        iion++;
         ia++;
     }
 }
@@ -2244,13 +2204,11 @@ void Ions::getLocalForces(std::vector<double>& tau) const
 {
     assert(tau.size() == 3 * local_ions_.size());
 
-    int ia                                 = 0;
-    std::vector<Ion*>::const_iterator iion = local_ions_.begin();
-    while (iion != local_ions_.end())
+    int ia = 0;
+    for (auto& ion : local_ions_)
     {
         assert(3 * ia + 2 < (int)tau.size());
-        (*iion)->getForce(&tau[3 * ia]);
-        iion++;
+        ion->getForce(&tau[3 * ia]);
         ia++;
     }
 }
@@ -2351,12 +2309,10 @@ double Ions::computeMaxVlRadius() const
 {
     double radius = 0.;
 
-    std::vector<Ion*>::const_iterator iion = local_ions_.begin();
-    while (iion != local_ions_.end())
+    for (auto& ion : local_ions_)
     {
-        double r = (*iion)->computeRadiusVl();
+        double r = ion->computeRadiusVl();
         radius   = r > radius ? r : radius;
-        iion++;
     }
 
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
@@ -2393,13 +2349,11 @@ void Ions::gatherNames(std::map<int, std::string>& names, const int root,
     std::vector<std::string> data;
     std::vector<std::string> local_names;
 
-    std::vector<Ion*>::const_iterator ion = local_ions_.begin();
-    while (ion != local_ions().end())
+    for (auto& ion : local_ions_)
     {
         // get local name and index
-        local_names.push_back((*ion)->name());
-        local_indexes.push_back((*ion)->index());
-        ++ion;
+        local_names.push_back(ion->name());
+        local_indexes.push_back(ion->index());
     }
 
     // gather data to PE root
@@ -2695,15 +2649,13 @@ bool Ions::hasLockedAtoms() const
     return (flag == 1);
 }
 
-double Ions::getMaxNLradius() const
+double Ions::getSpeciesMaxNLradius() const
 {
-    double radius                            = 0;
-    std::vector<Species>::const_iterator spi = species_.begin();
-    while (spi != species_.end())
+    double radius = 0.;
+    for (const auto& spi : species_)
     {
-        const double nlradius = spi->nlradius();
+        const double nlradius = spi.nlradius();
         radius                = radius > nlradius ? radius : nlradius;
-        spi++;
     }
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
     mmpi.allreduce(&radius, 1, MPI_MAX);
@@ -2711,15 +2663,13 @@ double Ions::getMaxNLradius() const
     return radius;
 }
 
-double Ions::getMaxLradius() const
+double Ions::getSpeciesMaxLradius() const
 {
-    double radius                            = 0;
-    std::vector<Species>::const_iterator spi = species_.begin();
-    while (spi != species_.end())
+    double radius = 0.;
+    for (const auto& spi : species_)
     {
-        const double lradius = spi->lradius();
+        const double lradius = spi.lradius();
         radius               = radius > lradius ? radius : lradius;
-        spi++;
     }
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
     mmpi.allreduce(&radius, 1, MPI_MAX);
@@ -2918,8 +2868,8 @@ void Ions::computeNumIons(void)
 double Ions::getMaxListRadius() const
 {
     // get radius of projectors
-    const double nlradius = getMaxNLradius();
-    const double lradius  = getMaxLradius();
+    const double nlradius = getSpeciesMaxNLradius();
+    const double lradius  = getSpeciesMaxLradius();
 
     double rmax = nlradius > lradius ? nlradius : lradius;
 
@@ -3289,21 +3239,15 @@ void Ions::updateIons()
     Control& ct(*(Control::instance()));
 
     // update local_ions data
-    std::vector<Ion*>::iterator ion = local_ions_.begin();
-    int ia                          = 0;
-    while (ion != local_ions_.end())
+    int ia = 0;
+    for (auto& ion : local_ions_)
     {
-        (*ion)->setPosition(
+        ion->setPosition(
             tau0_[3 * ia + 0], tau0_[3 * ia + 1], tau0_[3 * ia + 2]);
 
-        //        (*ion)->setForce(fion_[3*ia+0],
-        //                          fion_[3*ia+1],
-        //                          fion_[3*ia+2]);
-
-        (*ion)->setRandomState(rand_states_[3 * ia + 0],
-            rand_states_[3 * ia + 1], rand_states_[3 * ia + 2]);
+        ion->setRandomState(rand_states_[3 * ia + 0], rand_states_[3 * ia + 1],
+            rand_states_[3 * ia + 2]);
 
-        ion++;
         ia++;
     }
 
diff --git a/src/Ions.h b/src/Ions.h
index 772c51d1..ae4f6adb 100644
--- a/src/Ions.h
+++ b/src/Ions.h
@@ -71,13 +71,14 @@ class Ions
     void setMapVL();
     double computeMaxVlRadius() const;
     double computeMaxNLprojRadius() const;
-    double getMaxNLradius() const;
-    double getMaxLradius() const;
-    void updateListIons();
 
-    // compute boundary for box containing all atoms to be known on local
-    // processor that is values for list_boundary_left_, list_boundary_right_
-    void setupListIonsBoundaries(const double rmax);
+    /*
+     * Evaluate maximum pseudopotential radius among all species in class
+     */
+    double getSpeciesMaxNLradius() const;
+    double getSpeciesMaxLradius() const;
+
+    void updateListIons();
 
     void augmentIonsData(const int nsteps, const int dir, const int disp,
         const int locSize, const int maxLocSize, std::vector<IonData>& data,
@@ -163,6 +164,8 @@ class Ions
     bool hasLockedAtoms() const;
     void clearLists();
 
+    void rescaleVelocities(const double factor);
+
 public:
     Ions(const double lat[3], const std::vector<Species>& sp);
 
@@ -172,6 +175,10 @@ class Ions
 
     void setup();
 
+    // compute boundary for box containing all atoms to be known on local
+    // processor that is values for list_boundary_left_, list_boundary_right_
+    void setupListIonsBoundaries(const double rmax);
+
     std::vector<std::string>& getLocalNames() { return local_names_; }
     std::vector<double>& getTau0() { return tau0_; }
     std::vector<double>& getTaup() { return taup_; }
@@ -334,13 +341,15 @@ class Ions
 
     void updateForcesInteractingIons();
     void updateTaupInteractingIons();
-    void rescaleVelocities(const double factor);
 
     /*!
      * Calculate minimum distance between local pairs
      */
     double computeMinLocalSpacing() const;
 
+    void addIonToList(const Species& sp, const std::string& name,
+        const double crds[3], const double velocity[3], const bool lock);
+
     // void checkUnicityLocalIons();
 };
 
diff --git a/src/KBprojectorSparse.cc b/src/KBprojectorSparse.cc
index afd8871c..937b95d2 100644
--- a/src/KBprojectorSparse.cc
+++ b/src/KBprojectorSparse.cc
@@ -103,7 +103,7 @@ void KBprojectorSparse::setNLindex(
 
     for (int i = 0; i < size_nl; i++)
     {
-        assert(i < lnumpt);
+        // assert(i < lnumpt);
         if ((pvec[i] < iloc * lnumpt) || (pvec[i] >= (iloc + 1) * lnumpt))
         {
             (*MPIdata::sout) << " iloc=" << iloc << ", i=" << i
@@ -960,7 +960,7 @@ bool KBprojectorSparse::setIndexesAndProjectors()
             // get "pvec" and "is_in_domain"
             int icount = get_index_array(pvec, iloc, index_low, index_high);
             assert(icount <= nl3);
-            assert(icount <= mymesh->npointsPatch());
+            //            assert(icount <= mymesh->npointsPatch());
             assert(icount > 0);
 
             setNLindex(iloc, icount, pvec);
diff --git a/src/MGmol.cc b/src/MGmol.cc
index 8da0d331..acc2b944 100644
--- a/src/MGmol.cc
+++ b/src/MGmol.cc
@@ -837,7 +837,7 @@ void MGmol<OrbitalsType>::initNuc(Ions& ions)
     pot.initialize(ions);
 
     // Check compensating charges
-    double comp_rho = pot.getCharge(pot.rho_comp());
+    double comp_rho = getCharge(pot.rho_comp());
 
     if (onpe0 && ct.verbose > 1)
     {
diff --git a/src/MGmol.h b/src/MGmol.h
index 6a4cc830..a9bca4ba 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -182,7 +182,10 @@ class MGmol : public MGmolInterface
 
     /* access functions */
     OrbitalsType* getOrbitals() { return current_orbitals_; }
-    std::shared_ptr<Hamiltonian<OrbitalsType>> getHamiltonian() { return hamiltonian_; }
+    std::shared_ptr<Hamiltonian<OrbitalsType>> getHamiltonian()
+    {
+        return hamiltonian_;
+    }
     std::shared_ptr<Rho<OrbitalsType>> getRho() { return rho_; }
 
     void run() override;
diff --git a/src/PCGSolver.cc b/src/PCGSolver.cc
index 12103aaa..b3bd7ec3 100644
--- a/src/PCGSolver.cc
+++ b/src/PCGSolver.cc
@@ -11,14 +11,13 @@
 
 #include <iomanip>
 #include <iostream>
-using namespace std;
 
-template <class T, typename T2>
-void PCGSolver<T, T2>::clear()
+template <class T, typename ScalarType>
+void PCGSolver<T, ScalarType>::clear()
 {
-    for (short i = 0; i < (short)pc_oper_.size(); i++)
+    for (short i = 0; i < (short)precond_oper_.size(); i++)
     {
-        delete pc_oper_[i];
+        delete precond_oper_[i];
     }
     for (short i = 0; i < (short)gf_work_.size(); i++)
     {
@@ -35,26 +34,26 @@ void PCGSolver<T, T2>::clear()
         assert(gf_newv_[i] != nullptr);
         delete gf_newv_[i];
     }
-    // delete grids after pb::GridFunc<T2> objects since those
+    // delete grids after pb::GridFunc<ScalarType> objects since those
     // have data members references to grids
     for (short i = 0; i < (short)grid_.size(); i++)
     {
         delete grid_[i];
     }
-    pc_oper_.clear();
+    precond_oper_.clear();
     grid_.clear();
     gf_work_.clear();
     gf_rcoarse_.clear();
     gf_newv_.clear();
 }
 
-template <class T, typename T2>
-void PCGSolver<T, T2>::setupPrecon()
+template <class T, typename ScalarType>
+void PCGSolver<T, ScalarType>::setupPrecon()
 {
     // check if precon is already setup
     // Assumes operator does not change, hence
     // a single setup is sufficient
-    if (is_pc_setup_) return;
+    if (is_precond_setup_) return;
 
     // fine level
     pb::Grid* mygrid = new pb::Grid(oper_.grid());
@@ -63,7 +62,7 @@ void PCGSolver<T, T2>::setupPrecon()
 
     pb::Lap<POISSONPRECONDTYPE>* myoper
         = LapFactory<POISSONPRECONDTYPE>::createLap(*grid_[0], lap_type_);
-    pc_oper_.push_back(myoper);
+    precond_oper_.push_back(myoper);
 
     pb::GridFunc<POISSONPRECONDTYPE>* gf_work
         = new pb::GridFunc<POISSONPRECONDTYPE>(
@@ -92,7 +91,7 @@ void PCGSolver<T, T2>::setupPrecon()
 
         pb::Lap<POISSONPRECONDTYPE>* myoper
             = LapFactory<POISSONPRECONDTYPE>::createLap(*coarse_grid, 1);
-        pc_oper_.push_back(myoper);
+        precond_oper_.push_back(myoper);
 
         gf_work = new pb::GridFunc<POISSONPRECONDTYPE>(
             *coarse_grid, bc_[0], bc_[1], bc_[2]);
@@ -109,12 +108,13 @@ void PCGSolver<T, T2>::setupPrecon()
 
         mygrid = coarse_grid;
     }
-    is_pc_setup_ = true;
+    is_precond_setup_ = true;
 }
 
 // MG V-cycle with no mask
-template <class T, typename T2>
-void PCGSolver<T, T2>::preconSolve(pb::GridFunc<POISSONPRECONDTYPE>& gf_v,
+template <class T, typename ScalarType>
+void PCGSolver<T, ScalarType>::preconSolve(
+    pb::GridFunc<POISSONPRECONDTYPE>& gf_v,
     const pb::GridFunc<POISSONPRECONDTYPE>& gf_f, const short level)
 {
     //(*MPIdata::sout)<<"Preconditioning::mg() at level "<<level<<endl;
@@ -124,7 +124,7 @@ void PCGSolver<T, T2>::preconSolve(pb::GridFunc<POISSONPRECONDTYPE>& gf_v,
         ncycl = 4 > (nu1_ + nu2_) ? 4 : (nu1_ + nu2_);
     }
 
-    pb::Lap<POISSONPRECONDTYPE>* myoper = pc_oper_[level];
+    pb::Lap<POISSONPRECONDTYPE>* myoper = precond_oper_[level];
 
     // SMOOTHING
     for (short it = 0; it < ncycl; it++)
@@ -161,25 +161,23 @@ void PCGSolver<T, T2>::preconSolve(pb::GridFunc<POISSONPRECONDTYPE>& gf_v,
 }
 
 // Left Preconditioned CG
-template <class T, typename T2>
-bool PCGSolver<T, T2>::solve(pb::GridFunc<T2>& gf_phi, pb::GridFunc<T2>& gf_rhs)
+template <class T, typename ScalarType>
+bool PCGSolver<T, ScalarType>::solve(
+    pb::GridFunc<ScalarType>& gf_phi, const pb::GridFunc<ScalarType>& gf_rhs)
 {
     bool converged           = false;
     const pb::Grid& finegrid = gf_phi.grid();
 
     // initial data and residual - We assume a nonzero initial guess
-    pb::GridFunc<T2> lhs(finegrid, bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType> lhs(finegrid, bc_[0], bc_[1], bc_[2]);
     // scale initial guess with epsilon
     oper_.inv_transform(gf_phi);
     // compute initial residual: r := b - Ax
     /* compute Ax */
     oper_.apply(gf_phi, lhs);
     /* set r = b */
-    pb::GridFunc<T2> res(gf_rhs);
+    pb::GridFunc<ScalarType> res(gf_rhs);
     oper_.transform(res);
-    // Hartree units
-    const double hu = 4. * M_PI;
-    res *= hu;
     /* compute r = r - Ax */
     res -= lhs;
 
@@ -199,11 +197,11 @@ bool PCGSolver<T, T2>::solve(pb::GridFunc<T2>& gf_phi, pb::GridFunc<T2>& gf_rhs)
     /* preconditioning step */
     prec_z.setValues(0.);
     preconSolve(prec_z, prec_res, 0);
-    pb::GridFunc<T2> z(prec_z);
+    pb::GridFunc<ScalarType> z(prec_z);
 
     // conjugate vectors
-    pb::GridFunc<T2> p(prec_z);
-    pb::GridFunc<T2> ap(p.grid(), bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType> p(prec_z);
+    pb::GridFunc<ScalarType> ap(p.grid(), bc_[0], bc_[1], bc_[2]);
 
     double rtz = res.gdot(z);
 
@@ -251,13 +249,14 @@ bool PCGSolver<T, T2>::solve(pb::GridFunc<T2>& gf_phi, pb::GridFunc<T2>& gf_rhs)
 }
 
 // Left Preconditioned CG
-template <class T, typename T2>
-bool PCGSolver<T, T2>::solve(T2* phi, T2* rhs, const char dis)
+template <class T, typename ScalarType>
+bool PCGSolver<T, ScalarType>::solve(
+    ScalarType* phi, ScalarType* rhs, const char dis)
 {
-    pb::GridFunc<T2> gf_phi(oper_.grid(), bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType> gf_phi(oper_.grid(), bc_[0], bc_[1], bc_[2]);
     gf_phi.assign(phi, dis);
 
-    pb::GridFunc<T2> gf_work(oper_.grid(), bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType> gf_work(oper_.grid(), bc_[0], bc_[1], bc_[2]);
     gf_work.assign(rhs, dis);
 
     bool converged = solve(gf_phi, gf_work);
diff --git a/src/PCGSolver.h b/src/PCGSolver.h
index d38787f1..f6db7cc6 100644
--- a/src/PCGSolver.h
+++ b/src/PCGSolver.h
@@ -7,8 +7,8 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-#ifndef _PCG_SOLVER_H_
-#define _PCG_SOLVER_H_
+#ifndef MGMOL_PCG_SOLVER_H
+#define MGMOL_PCG_SOLVER_H
 
 #include "Control.h"
 #include "Lap.h"
@@ -17,32 +17,41 @@
 
 #include <vector>
 
-template <class T, typename T2>
+template <class T, typename ScalarType>
 class PCGSolver
 {
-
 private:
     std::vector<pb::Grid*> grid_;
     short lap_type_;
     short bc_[3];
     bool fully_periodic_;
-    // operators
+
+    // operator to solve for
     T oper_;
-    std::vector<pb::Lap<POISSONPRECONDTYPE>*> pc_oper_;
+
+    // preconditioner operator for each MG level
+    std::vector<pb::Lap<POISSONPRECONDTYPE>*> precond_oper_;
     std::vector<pb::GridFunc<POISSONPRECONDTYPE>*> gf_work_;
     std::vector<pb::GridFunc<POISSONPRECONDTYPE>*> gf_rcoarse_;
     std::vector<pb::GridFunc<POISSONPRECONDTYPE>*> gf_newv_;
-    // solver params
+
+    // solver parameters
     int maxiters_;
     double tol_;
     double final_residual_;
     double residual_reduction_;
-    // precon params
+
+    // preconditioner parameters
     short nu1_;
     short nu2_;
     short max_nlevels_;
     short nlevels_;
-    bool is_pc_setup_;
+    bool is_precond_setup_;
+
+    void preconSolve(pb::GridFunc<POISSONPRECONDTYPE>& gf_v,
+        const pb::GridFunc<POISSONPRECONDTYPE>& gf_f, const short level = 0);
+    void setupPrecon();
+    void clear();
 
 public:
     PCGSolver(T& oper, const short px, const short py, const short pz)
@@ -61,15 +70,10 @@ class PCGSolver
         bc_[1]          = py;
         bc_[2]          = pz;
         fully_periodic_ = ((bc_[0] == 1) && (bc_[1] == 1) && (bc_[2] == 1));
-        // fine grid info
-        //        pb::Grid* mygrid=new pb::Grid(oper.grid());
-        //        grid_.push_back(mygrid);
-        // fine grid operator
-        Control& ct  = *(Control::instance());
-        lap_type_    = ct.lap_type;
-        is_pc_setup_ = false;
-        //        pb::Lap* myoper = LapFactory::createLap(*grid_[0],lap_type_);
-        //        pc_oper_.push_back(myoper);
+
+        Control& ct       = *(Control::instance());
+        lap_type_         = ct.lap_type;
+        is_precond_setup_ = false;
     };
 
     void setup(const short nu1, const short nu2, const short max_sweeps,
@@ -83,16 +87,10 @@ class PCGSolver
         setupPrecon();
     }
 
-    void clear();
-
-    void setupPrecon();
-
-    void preconSolve(pb::GridFunc<POISSONPRECONDTYPE>& gf_v,
-        const pb::GridFunc<POISSONPRECONDTYPE>& gf_f, const short level = 0);
-
-    bool solve(pb::GridFunc<T2>& gf_phi, pb::GridFunc<T2>& gf_rhs);
+    bool solve(pb::GridFunc<ScalarType>& gf_phi,
+        const pb::GridFunc<ScalarType>& gf_rhs);
 
-    bool solve(T2* phi, T2* rhs, const char dis);
+    bool solve(ScalarType* phi, ScalarType* rhs, const char dis);
 
     double getFinalResidual() const { return final_residual_; }
     double getResidualReduction() const { return residual_reduction_; }
diff --git a/src/PCGSolver_Diel.cc b/src/PCGSolver_Diel.cc
index 0cdc5aed..4e7c88a9 100644
--- a/src/PCGSolver_Diel.cc
+++ b/src/PCGSolver_Diel.cc
@@ -9,10 +9,8 @@
 
 #include "PCGSolver_Diel.h"
 
-using namespace std;
-
-template <class T, typename T2>
-void PCGSolver_Diel<T, T2>::clear()
+template <class T, typename ScalarType>
+void PCGSolver_Diel<T, ScalarType>::clear()
 {
     for (short i = 0; i < (short)pc_oper_.size(); i++)
     {
@@ -33,7 +31,7 @@ void PCGSolver_Diel<T, T2>::clear()
         assert(gf_newv_[i] != nullptr);
         delete gf_newv_[i];
     }
-    // delete grids after pb::GridFunc<T2> objects since those
+    // delete grids after pb::GridFunc<ScalarType> objects since those
     // have data members references to grids
     for (short i = 0; i < (short)grid_.size(); i++)
     {
@@ -46,8 +44,8 @@ void PCGSolver_Diel<T, T2>::clear()
     gf_newv_.clear();
 }
 
-template <class T, typename T2>
-void PCGSolver_Diel<T, T2>::setupPrecon()
+template <class T, typename ScalarType>
+void PCGSolver_Diel<T, ScalarType>::setupPrecon()
 {
     // fine level
     pb::Grid* mygrid = new pb::Grid(oper_.grid());
@@ -57,8 +55,8 @@ void PCGSolver_Diel<T, T2>::setupPrecon()
     T* myoper = new T(oper_);
     pc_oper_.push_back(myoper);
 
-    pb::GridFunc<T2>* gf_work
-        = new pb::GridFunc<T2>(*grid_[0], bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType>* gf_work
+        = new pb::GridFunc<ScalarType>(*grid_[0], bc_[0], bc_[1], bc_[2]);
     gf_work_.push_back(gf_work);
 
     // coarse levels
@@ -86,24 +84,25 @@ void PCGSolver_Diel<T, T2>::setupPrecon()
         T* myoper = new T(pc_oper_[ln - 1]->coarseOp(*mygrid));
         pc_oper_.push_back(myoper);
 
-        gf_work = new pb::GridFunc<T2>(*coarse_grid, bc_[0], bc_[1], bc_[2]);
+        gf_work = new pb::GridFunc<ScalarType>(
+            *coarse_grid, bc_[0], bc_[1], bc_[2]);
         gf_work_.push_back(gf_work);
 
-        pb::GridFunc<T2>* gf_rcoarse
-            = new pb::GridFunc<T2>(*coarse_grid, bc_[0], bc_[1], bc_[2]);
+        pb::GridFunc<ScalarType>* gf_rcoarse = new pb::GridFunc<ScalarType>(
+            *coarse_grid, bc_[0], bc_[1], bc_[2]);
         gf_rcoarse_.push_back(gf_rcoarse);
-        pb::GridFunc<T2>* gf_newv
-            = new pb::GridFunc<T2>(*coarse_grid, bc_[0], bc_[1], bc_[2]);
+        pb::GridFunc<ScalarType>* gf_newv = new pb::GridFunc<ScalarType>(
+            *coarse_grid, bc_[0], bc_[1], bc_[2]);
         gf_newv_.push_back(gf_newv);
 
         mygrid = coarse_grid;
     }
 }
 
-template <class T, typename T2>
+template <class T, typename ScalarType>
 // MG V-cycle with no mask
-void PCGSolver_Diel<T, T2>::preconSolve(
-    pb::GridFunc<T2>& gf_v, const pb::GridFunc<T2>& gf_f, const short level)
+void PCGSolver_Diel<T, ScalarType>::preconSolve(pb::GridFunc<ScalarType>& gf_v,
+    const pb::GridFunc<ScalarType>& gf_f, const short level)
 {
     //(*MPIdata::sout)<<"Preconditioning::mg() at level "<<level<<endl;
     short ncycl = nu1_;
@@ -126,11 +125,11 @@ void PCGSolver_Diel<T, T2>::preconSolve(
     // COARSE GRID CORRECTION
 
     // restrictions
-    pb::GridFunc<T2>* rcoarse = gf_rcoarse_[level];
+    pb::GridFunc<ScalarType>* rcoarse = gf_rcoarse_[level];
     gf_work_[level]->restrict3D(*rcoarse);
 
     // storage functions for coarse grid
-    pb::GridFunc<T2>* newv = gf_newv_[level];
+    pb::GridFunc<ScalarType>* newv = gf_newv_[level];
 
     // call mgrid solver on a coarser level
     newv->resetData();
@@ -149,15 +148,16 @@ void PCGSolver_Diel<T, T2>::preconSolve(
     if (bc_[0] != 1 || bc_[2] != 1 || bc_[2] != 1) gf_v.trade_boundaries();
 }
 
-template <class T, typename T2>
+template <class T, typename ScalarType>
 // Left Preconditioned CG
-bool PCGSolver_Diel<T, T2>::solve(
-    pb::GridFunc<T2>& gf_phi, pb::GridFunc<T2>& gf_rhs)
+bool PCGSolver_Diel<T, ScalarType>::solve(
+    pb::GridFunc<ScalarType>& gf_phi, pb::GridFunc<ScalarType>& gf_rhs)
 {
     if (!oper_.initialized())
     {
-        cout << "Error in PCGSolver_Diel<T>::solve: operator not initialized"
-             << endl;
+        std::cout
+            << "Error in PCGSolver_Diel<T>::solve: operator not initialized"
+            << std::endl;
         return 0.;
     }
 
@@ -165,13 +165,13 @@ bool PCGSolver_Diel<T, T2>::solve(
     const pb::Grid& finegrid = gf_phi.grid();
 
     // initial data and residual - We assume a nonzero initial guess
-    pb::GridFunc<T2> lhs(finegrid, bc_[0], bc_[1], bc_[2]);
-    pb::GridFunc<T2> res(finegrid, bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType> lhs(finegrid, bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType> res(finegrid, bc_[0], bc_[1], bc_[2]);
     // scale initial guess with epsilon
     oper_.inv_transform(gf_phi);
     // compute initial residual
     oper_.apply(gf_phi, lhs);
-    pb::GridFunc<T2> rhs(gf_rhs);
+    pb::GridFunc<ScalarType> rhs(gf_rhs);
     oper_.transform(rhs);
     // Hartree units
     rhs *= (4. * M_PI);
@@ -180,13 +180,13 @@ bool PCGSolver_Diel<T, T2>::solve(
     double rnorm      = init_rnorm;
 
     // preconditioned residual
-    pb::GridFunc<T2> z(finegrid, bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType> z(finegrid, bc_[0], bc_[1], bc_[2]);
     // preconditioning step
     z = 0.;
     preconSolve(z, res, 0);
     // conjugate vectors
-    pb::GridFunc<T2> p(z);
-    pb::GridFunc<T2> ap(p.grid(), bc_[0], bc_[1], bc_[2]);
+    pb::GridFunc<ScalarType> p(z);
+    pb::GridFunc<ScalarType> ap(p.grid(), bc_[0], bc_[1], bc_[2]);
 
     double rtz = res.gdot(z);
 
@@ -225,11 +225,11 @@ bool PCGSolver_Diel<T, T2>::solve(
     return converged;
 }
 
-template <class T, typename T2>
+template <class T, typename ScalarType>
 // Left Preconditioned CG
-bool PCGSolver_Diel<T, T2>::solve(pb::GridFunc<T2>& gf_phi,
-    pb::GridFunc<T2>& gf_rhs, pb::GridFunc<T2>& gf_rhod,
-    pb::GridFunc<T2>& gf_vks)
+bool PCGSolver_Diel<T, ScalarType>::solve(pb::GridFunc<ScalarType>& gf_phi,
+    pb::GridFunc<ScalarType>& gf_rhs, pb::GridFunc<ScalarType>& gf_rhod,
+    pb::GridFunc<ScalarType>& gf_vks)
 {
     // initialize the linear system operator and the preconditioner
     oper_.init(gf_rhod);
diff --git a/src/PCGSolver_Diel.h b/src/PCGSolver_Diel.h
index 9fd97914..f2001a4c 100644
--- a/src/PCGSolver_Diel.h
+++ b/src/PCGSolver_Diel.h
@@ -7,8 +7,8 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-#ifndef _PCG_SOLVER_DIEL_H_
-#define _PCG_SOLVER_DIEL_H_
+#ifndef MGMOL_PCG_SOLVER_DIEL_H_
+#define MGMOL_PCG_SOLVER_DIEL_H_
 
 #include "Control.h"
 #include "PB.h"
@@ -21,7 +21,7 @@
 
 #include <vector>
 
-template <class T, typename T2>
+template <class T, typename ScalarType>
 class PCGSolver_Diel
 {
 
@@ -32,9 +32,9 @@ class PCGSolver_Diel
     // operators
     T oper_;
     std::vector<T*> pc_oper_;
-    std::vector<pb::GridFunc<T2>*> gf_work_;
-    std::vector<pb::GridFunc<T2>*> gf_rcoarse_;
-    std::vector<pb::GridFunc<T2>*> gf_newv_;
+    std::vector<pb::GridFunc<ScalarType>*> gf_work_;
+    std::vector<pb::GridFunc<ScalarType>*> gf_rcoarse_;
+    std::vector<pb::GridFunc<ScalarType>*> gf_newv_;
     // solver params
     int maxiters_;
     double tol_;
@@ -47,6 +47,11 @@ class PCGSolver_Diel
     short nlevels_;
     void setupPrecon();
 
+    void clear();
+
+    void preconSolve(pb::GridFunc<ScalarType>& gf_v,
+        const pb::GridFunc<ScalarType>& gf_f, const short level = 0);
+
 public:
     PCGSolver_Diel(T& oper, const short px, const short py, const short pz)
         : oper_(oper)
@@ -78,15 +83,12 @@ class PCGSolver_Diel
         max_nlevels_ = max_nlevels;
     }
 
-    void clear();
-
-    void preconSolve(pb::GridFunc<T2>& gf_v, const pb::GridFunc<T2>& gf_f,
-        const short level = 0);
-
-    bool solve(pb::GridFunc<T2>& gf_phi, pb::GridFunc<T2>& gf_rhs);
+    bool solve(
+        pb::GridFunc<ScalarType>& gf_phi, pb::GridFunc<ScalarType>& gf_rhs);
 
-    bool solve(pb::GridFunc<T2>& gf_phi, pb::GridFunc<T2>& gf_rhs,
-        pb::GridFunc<T2>& gf_rhod, pb::GridFunc<T2>& gf_vks);
+    bool solve(pb::GridFunc<ScalarType>& gf_phi,
+        pb::GridFunc<ScalarType>& gf_rhs, pb::GridFunc<ScalarType>& gf_rhod,
+        pb::GridFunc<ScalarType>& gf_vks);
 
     double getFinalResidual() const { return final_residual_; }
     double getResidualReduction() const { return residual_reduction_; }
diff --git a/src/PoissonSolverFactory.h b/src/PoissonSolverFactory.h
new file mode 100644
index 00000000..91b99bf2
--- /dev/null
+++ b/src/PoissonSolverFactory.h
@@ -0,0 +1,228 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC and
+// UT-Battelle, LLC.
+// Produced at the Lawrence Livermore National Laboratory and the Oak Ridge
+// National Laboratory.
+// LLNL-CODE-743438
+// All rights reserved.
+// This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
+// Please also read this link https://github.com/llnl/mgmol/LICENSE
+#ifndef MGMOL_PoissonSolverFactory
+#define MGMOL_PoissonSolverFactory
+
+#include "Control.h"
+#include "Hartree.h"
+#include "Hartree_CG.h"
+#include "Mesh.h"
+#include "PBdiel.h"
+#include "PBdiel_CG.h"
+#include "ShiftedHartree.h"
+#include "mputils.h"
+
+#include "GridFunc.h"
+#include "Laph2.h"
+#include "Laph4.h"
+#include "Laph4M.h"
+#include "Laph4MP.h"
+#include "Laph6.h"
+#include "Laph8.h"
+#include "ShiftedLaph4M.h"
+
+class PoissonSolverFactory
+{
+
+public:
+    /*!
+     * return specific Poisson solver needed to solve Hartree problem
+     */
+    static Poisson* create(const pb::Grid& myGrid, PoissonFDtype lap_type,
+        const short bc[3], const double screening_const)
+    {
+        Poisson* poisson_solver = nullptr;
+
+        Control& ct = *(Control::instance());
+        if (ct.MGPoissonSolver()) // use MG for Poisson Solver
+        {
+            if (screening_const > 0.)
+            {
+                switch (lap_type)
+                {
+                    case PoissonFDtype::h4M:
+                        poisson_solver
+                            = new ShiftedHartree<pb::ShiftedLaph4M<POTDTYPE>>(
+                                myGrid, bc, screening_const);
+                        break;
+                    default:
+                        (*MPIdata::sout)
+                            << "Electrostatic, shifted, Undefined option: "
+                            << static_cast<int>(lap_type) << std::endl;
+                }
+            }
+            else
+            {
+                switch (lap_type)
+                {
+                    case PoissonFDtype::h4M:
+                        poisson_solver
+                            = new Hartree<pb::Laph4M<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h2:
+                        poisson_solver
+                            = new Hartree<pb::Laph2<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h4:
+                        poisson_solver
+                            = new Hartree<pb::Laph4<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h6:
+                        poisson_solver
+                            = new Hartree<pb::Laph6<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h8:
+                        poisson_solver
+                            = new Hartree<pb::Laph8<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h4MP:
+                        poisson_solver
+                            = new Hartree<pb::Laph4MP<POTDTYPE>>(myGrid, bc);
+                        break;
+                    default:
+                        (*MPIdata::sout)
+                            << "Electrostatic, Undefined option: "
+                            << static_cast<int>(lap_type) << std::endl;
+                }
+            }
+        }
+        else // use PCG for Poisson Solver
+        {
+            if (screening_const > 0.)
+            {
+                switch (lap_type)
+                {
+                    case PoissonFDtype::h4M:
+                        poisson_solver
+                            = new ShiftedHartree<pb::ShiftedLaph4M<POTDTYPE>>(
+                                myGrid, bc, screening_const);
+                        break;
+                    default:
+                        (*MPIdata::sout)
+                            << "PCG Electrostatic, shifted, Undefined option: "
+                            << static_cast<int>(lap_type) << std::endl;
+                }
+            }
+            else
+            {
+                switch (lap_type)
+                {
+                    case PoissonFDtype::h4M:
+                        poisson_solver
+                            = new Hartree_CG<pb::Laph4M<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h2:
+                        poisson_solver
+                            = new Hartree_CG<pb::Laph2<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h4:
+                        poisson_solver
+                            = new Hartree_CG<pb::Laph4<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h6:
+                        poisson_solver
+                            = new Hartree_CG<pb::Laph6<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h8:
+                        poisson_solver
+                            = new Hartree_CG<pb::Laph8<POTDTYPE>>(myGrid, bc);
+                        break;
+                    case PoissonFDtype::h4MP:
+                        poisson_solver
+                            = new Hartree_CG<pb::Laph4MP<POTDTYPE>>(myGrid, bc);
+                        break;
+                    default:
+                        (*MPIdata::sout)
+                            << "PCG Electrostatic, Undefined option: "
+                            << static_cast<int>(lap_type) << std::endl;
+                }
+            }
+        }
+
+        return poisson_solver;
+    }
+
+    static Poisson* createDiel(pb::Grid& pbGrid, PoissonFDtype lap_type,
+        const short bc[3], const double e0, const double rho0,
+        const double drho0)
+    {
+        Poisson* poisson_solver = nullptr;
+
+        Control& ct = *(Control::instance());
+        if (ct.MGPoissonSolver()) // use MG for Poisson Solver
+        {
+            switch (lap_type)
+            {
+                case PoissonFDtype::h4M:
+                    poisson_solver = new PBdiel<pb::PBh4M<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h2:
+                    poisson_solver = new PBdiel<pb::PBh2<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h4:
+                    poisson_solver = new PBdiel<pb::PBh4<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h6:
+                    poisson_solver = new PBdiel<pb::PBh6<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h8:
+                    poisson_solver = new PBdiel<pb::PBh8<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h4MP:
+                    poisson_solver = new PBdiel<pb::PBh4MP<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                default:
+                    (*MPIdata::sout)
+                        << "Electrostatic, Undefined option" << std::endl;
+            }
+        }
+        else // use PCG for Poisson Solver
+        {
+            switch (lap_type)
+            {
+                case PoissonFDtype::h4M:
+                    poisson_solver = new PBdiel_CG<pb::PBh4M<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h2:
+                    poisson_solver = new PBdiel_CG<pb::PBh2<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h4:
+                    poisson_solver = new PBdiel_CG<pb::PBh4<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h6:
+                    poisson_solver = new PBdiel_CG<pb::PBh6<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h8:
+                    poisson_solver = new PBdiel_CG<pb::PBh8<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                case PoissonFDtype::h4MP:
+                    poisson_solver = new PBdiel_CG<pb::PBh4MP<POTDTYPE>>(
+                        pbGrid, bc, e0, rho0, drho0);
+                    break;
+                default:
+                    (*MPIdata::sout)
+                        << "Electrostatic, Undefined option" << std::endl;
+            }
+        }
+        return poisson_solver;
+    }
+};
+
+#endif
diff --git a/src/Potentials.cc b/src/Potentials.cc
index ee874f10..380f05a0 100644
--- a/src/Potentials.cc
+++ b/src/Potentials.cc
@@ -27,6 +27,7 @@
 #include <fstream>
 using namespace std;
 
+// unit conversion factor Ha -> Ry
 const double ha2ry = 2.;
 
 Potentials::~Potentials()
@@ -57,13 +58,9 @@ Potentials::Potentials(const bool vh_frozen)
 
     size_ = dim_[0] * dim_[1] * dim_[2];
 
-    mix_ = 1.;
-
     scf_dvrho_ = 1000.;
     scf_dv_    = 1000.;
 
-    vh_frozen_ = vh_frozen;
-
     vtot_.resize(size_);
     vtot_old_.resize(size_);
 
@@ -101,7 +98,7 @@ void Potentials::initWithVnuc()
     double one = 1.;
     LinearAlgebraUtils<MemorySpace::Host>::MPaxpy(
         size_, one, &v_ext_[0], &vtot_[0]);
-    // factor 2 to get total potential in [Ry] for calculations
+    // factor ha2ry to get total potential in [Ry] for calculations
     LinearAlgebraUtils<MemorySpace::Host>::MPscal(size_, ha2ry, &vtot_[0]);
 }
 
@@ -568,9 +565,8 @@ template <typename T>
 void Potentials::setVxc(const T* const vxc, const int iterativeIndex)
 {
     assert(iterativeIndex >= 0);
-    //    int ione=1;
+
     itindex_vxc_ = iterativeIndex;
-    //    Tcopy(&size_, vxc,  &ione, &vxc_rho_[0], &ione);
     MPcpy(&vxc_rho_[0], vxc, size_);
 }
 void Potentials::setVh(const POTDTYPE* const vh, const int iterativeIndex)
@@ -849,20 +845,6 @@ void Potentials::rescaleRhoComp()
     if (comp_rho < 0.) mmpi.abort();
 }
 
-double Potentials::getCharge(RHODTYPE* rho)
-{
-    Control& ct            = *(Control::instance());
-    Mesh* mymesh           = Mesh::instance();
-    const pb::Grid& mygrid = mymesh->grid();
-
-    double charge = mygrid.integralOverMesh(rho);
-
-    if (onpe0 && ct.verbose > 0)
-        cout << setprecision(8) << fixed << "Charge: " << charge << endl;
-
-    return charge;
-}
-
 void Potentials::addBackgroundToRhoComp()
 {
     if (fabs(background_charge_) > 0.)
diff --git a/src/Potentials.h b/src/Potentials.h
index 0cd7a8a6..205fb538 100644
--- a/src/Potentials.h
+++ b/src/Potentials.h
@@ -30,7 +30,6 @@ class Potentials
     int gdim_[3];
     int dim_[3];
     bool diel_;
-    double mix_;
 
     double scf_dvrho_;
     double scf_dv_;
@@ -39,8 +38,9 @@ class Potentials
     double charge_in_cell_;
     double ionic_charge_;
 
-    bool vh_frozen_;
-
+    /*!
+     * Total KS potential seen by electrons
+     */
     std::vector<POTDTYPE> vtot_;
     std::vector<POTDTYPE> vtot_old_;
 
@@ -48,16 +48,30 @@ class Potentials
     std::vector<POTDTYPE> vh_rho_;
     std::vector<POTDTYPE> vxc_rho_;
 
-    // nuclei local potential
+    /*
+     * Potential contribution from atomic cores (local pseudopotential)
+     */
     std::vector<POTDTYPE> v_nuc_;
 
-    // external potential (read from input)
+    /*!
+     * Optional external potential (read from input)
+     * Used only in special cases.
+     */
     std::vector<POTDTYPE> v_ext_;
 #ifdef HAVE_TRICUBIC
     pb::TriCubic<POTDTYPE>* vext_tricubic_;
 #endif
 
+    /*!
+     * Potential associated with the sum of Gaussian charge distributions
+     * compensating  the Coulomb potential of each atom
+     */
     std::vector<POTDTYPE> v_comp_;
+
+    /*!
+     * Sum of Gaussian charge distributions compensating the Coulomb potential
+     * of each atom
+     */
     std::vector<RHODTYPE> rho_comp_;
 
     std::vector<POTDTYPE> dv_;
@@ -121,14 +135,11 @@ class Potentials
     void turnOnDiel() { diel_ = true; }
 
     int size() const { return size_; }
-    bool vh_frozen() const { return vh_frozen_; }
-    void freeze_vh() { vh_frozen_ = true; }
 
     double scf_dvrho(void) const { return scf_dvrho_; }
     double scf_dv(void) const { return scf_dv_; }
     POTDTYPE* vtot() { return &vtot_[0]; }
     POTDTYPE* vh_rho() { return &vh_rho_[0]; }
-    POTDTYPE* vxc_rho() { return &vxc_rho_[0]; }
     RHODTYPE* rho_comp() { return &rho_comp_[0]; }
 
     const std::vector<POTDTYPE>& vnuc() const { return v_nuc_; }
@@ -136,12 +147,6 @@ class Potentials
     POTDTYPE* vext() { return &v_ext_[0]; }
     POTDTYPE* vepsilon() { return &vepsilon_[0]; }
 
-    void set_vcomp(const POTDTYPE val)
-    {
-        const int n = (int)v_comp_.size();
-        for (int i = 0; i < n; i++)
-            v_comp_[i] = val;
-    }
     void axpVcompToVh(const double alpha);
     void axpVcomp(POTDTYPE* v, const double alpha);
 
@@ -154,13 +159,29 @@ class Potentials
 
     double getChargeInCell() const { return charge_in_cell_; }
 
+    /*!
+     * initialize total potential as local pseudopotential
+     */
     void initWithVnuc();
 
     void getVofRho(std::vector<POTDTYPE>& vrho) const;
 
-    double delta_v(const std::vector<std::vector<RHODTYPE>>&);
-    double update(const std::vector<std::vector<RHODTYPE>>&);
-    void update(const double);
+    /*!
+     * evaluate potential correction associated with a new rho
+     */
+    double delta_v(const std::vector<std::vector<RHODTYPE>>& rho);
+
+    /*!
+     * update potentials based on argument rho
+     */
+    double update(const std::vector<std::vector<RHODTYPE>>& rho);
+
+    /*!
+     * update potentials based on potential correction delta v and mixing
+     * parameter
+     */
+    void update(const double mix);
+
     double max() const;
     double min() const;
     void readAll(std::vector<Species>& sp);
@@ -172,7 +193,6 @@ class Potentials
 
     void initialize(Ions& ions);
     void rescaleRhoComp();
-    double getCharge(RHODTYPE* rho);
     void initBackground(Ions& ions);
     void addBackgroundToRhoComp();
 
diff --git a/src/ProjectedMatrices.cc b/src/ProjectedMatrices.cc
index e0e679e9..3606aa40 100644
--- a/src/ProjectedMatrices.cc
+++ b/src/ProjectedMatrices.cc
@@ -509,6 +509,7 @@ template <class MatrixType>
 void ProjectedMatrices<MatrixType>::setOccupations(
     const std::vector<DISTMATDTYPE>& occ)
 {
+    assert(!occ.empty());
 #ifdef PRINT_OPERATIONS
     if (mmpi.instancePE0())
         (*MPIdata::sout) << "ProjectedMatrices<MatrixType>::setOccupations()"
diff --git a/src/pb/Mgm.h b/src/pb/Mgm.h
index d16f5610..4b5ad68d 100644
--- a/src/pb/Mgm.h
+++ b/src/pb/Mgm.h
@@ -7,7 +7,6 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-// $Id: Mgm.h,v 1.13 2010/01/28 22:56:47 jeanluc Exp $
 #ifndef PB_MGM_H
 #define PB_MGM_H
 
@@ -41,8 +40,6 @@ bool Mgm(T1& A, T2& vh, const GridFunc<T3>& rho, const short cogr,
     A.rhs(res, rhs);
 
     // Hartree units
-    rhs *= (4. * M_PI);
-
     // work GridFunc<T3>
     GridFunc<T3> lhs(finegrid, bcx, bcy, bcz);
 
diff --git a/src/pb/Solver.h b/src/pb/Solver.h
index 062fbaa2..f82b3b6a 100644
--- a/src/pb/Solver.h
+++ b/src/pb/Solver.h
@@ -7,9 +7,8 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-// $Id: Solver.h,v 1.14 2010/01/28 22:56:47 jeanluc Exp $
-#ifndef SOLVE_H
-#define SOLVE_H
+#ifndef PB_SOLVER_H
+#define PB_SOLVER_H
 
 #include "GridFunc.h"
 
@@ -33,7 +32,7 @@ class Solver
         fully_periodic_ = ((bc_[0] == 1) && (bc_[1] == 1) && (bc_[2] == 1));
     }
 
-    virtual bool solve(GridFunc<T>&, GridFunc<T>&) = 0;
+    virtual bool solve(GridFunc<T>&, const GridFunc<T>&) = 0;
 
     virtual void setup(const short nu1, const short nu2, const short max_sweeps,
         const double tol, const short max_nlevels,
diff --git a/src/pb/SolverLap.cc b/src/pb/SolverLap.cc
index 5a85fefe..062755a8 100644
--- a/src/pb/SolverLap.cc
+++ b/src/pb/SolverLap.cc
@@ -7,8 +7,6 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-// $Id: SolverLap.cc,v 1.15 2010/01/28 22:56:31 jeanluc Exp $
-
 #include "SolverLap.h"
 #include "Laph2.h"
 #include "Laph4.h"
@@ -26,9 +24,8 @@ Timer vcycle_repvcycle_tm("vcycle_repvcycle");
 
 namespace pb
 {
-
 // explicit instantiation declaration
-#if 1
+
 // double
 template class SolverLap<ShiftedLaph4M<double>, double>;
 template class SolverLap<Laph4MP<double>, double>;
@@ -45,77 +42,7 @@ template class SolverLap<Laph4<float>, float>;
 template class SolverLap<Laph2<float>, float>;
 template class SolverLap<Laph6<float>, float>;
 template class SolverLap<Laph8<float>, float>;
-#else
-// double
-template bool SolverLap<Laph4M<double>, double>::solve(
-    double*, double*, const char, const int, const int);
-template bool SolverLap<Laph2<double>, double>::solve(
-    double*, double*, const char, const int, const int);
-template bool SolverLap<Laph4<double>, double>::solve(
-    double*, double*, const char, const int, const int);
-template bool SolverLap<Laph6<double>, double>::solve(
-    double*, double*, const char, const int, const int);
-
-template bool SolverLap<Laph4M<double>, double>::solve(
-    GridFunc<double>&, GridFunc<double>&, const int);
-template bool SolverLap<Laph2<double>, double>::solve(
-    GridFunc<double>&, GridFunc<double>&, const int);
-template bool SolverLap<Laph4<double>, double>::solve(
-    GridFunc<double>&, GridFunc<double>&, const int);
-template bool SolverLap<Laph6<double>, double>::solve(
-    GridFunc<double>&, GridFunc<double>&, const int);
-// float
-template bool SolverLap<Laph4M<float>, float>::solve(
-    float*, float*, const char, const int, const int);
-template bool SolverLap<Laph2<float>, float>::solve(
-    float*, float*, const char, const int, const int);
-template bool SolverLap<Laph4<float>, float>::solve(
-    float*, float*, const char, const int, const int);
-template bool SolverLap<Laph6<float>, float>::solve(
-    float*, float*, const char, const int, const int);
 
-template bool SolverLap<Laph4M<float>, float>::solve(
-    GridFunc<float>&, GridFunc<float>&, const int);
-template bool SolverLap<Laph2<float>, float>::solve(
-    GridFunc<float>&, GridFunc<float>&, const int);
-template bool SolverLap<Laph4<float>, float>::solve(
-    GridFunc<float>&, GridFunc<float>&, const int);
-template bool SolverLap<Laph6<float>, float>::solve(
-    GridFunc<float>&, GridFunc<float>&, const int);
-#endif
-/*
-template bool Mgm(Laph4MP&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE>&, const
-short, const short, const double, const short, const short, const bool,
-                    double&,double&,double&,short&);
-template bool Mgm(Laph4M&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE>&, const
-short, const short, const double, const short, const short, const bool,
-                    double&,double&,double&,short&);
-template bool Mgm(Laph2&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE>&, const
-short, const short, const double, const short, const short, const bool,
-                    double&,double&,double&,short&);
-template bool Mgm(Laph4&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE>&, const
-short, const short, const double, const short, const short, const bool,
-                    double&,double&,double&,short&);
-template bool Mgm(Laph6&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE>&, const
-short, const short, const double, const short, const short, const bool,
-                    double&,double&,double&,short&);
-template bool Mgm(ShiftedLaph4M&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE>&,
-const short, const short, const double, const short, const short, const bool,
-                    double&,double&,double&,short&);
-*/
-/*
-template int Vcycle(ShiftedLaph4M&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE>
-&, const short, const short, const short, const bool); template int
-Vcycle(Laph4MP&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE> &, const short,
-const short, const short, const bool); template int Vcycle(Laph4M&,
-GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE> &, const short, const short, const
-short, const bool); template int Vcycle(Laph2&, GridFunc<GFDTYPE>&, const
-GridFunc<GFDTYPE> &, const short, const short, const short, const bool);
-template int Vcycle(Laph4&, GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE> &, const
-short, const short, const short, const bool); template int Vcycle(Laph6&,
-GridFunc<GFDTYPE>&, const GridFunc<GFDTYPE> &, const short, const short, const
-short, const bool);
-*/
 template <class T, typename T2>
 bool SolverLap<T, T2>::solve(T2* phi, T2* rhs, const char dis)
 {
@@ -139,7 +66,7 @@ bool SolverLap<T, T2>::solve(T2* phi, T2* rhs, const char dis)
 }
 
 template <class T, typename T2>
-bool SolverLap<T, T2>::solve(GridFunc<T2>& gf_phi, GridFunc<T2>& gf_rhs)
+bool SolverLap<T, T2>::solve(GridFunc<T2>& gf_phi, const GridFunc<T2>& gf_rhs)
 {
     bool conv = Mgm(oper_, gf_phi, gf_rhs, max_nlevels_, max_sweeps_, tol_,
         nu1_, nu2_, gather_coarse_level_, final_residual_,
diff --git a/src/pb/SolverLap.h b/src/pb/SolverLap.h
index 939e15bc..567bed05 100644
--- a/src/pb/SolverLap.h
+++ b/src/pb/SolverLap.h
@@ -64,7 +64,7 @@ class SolverLap : public Solver<T2>
 
     bool solve(T2* phi, T2* rhs, const char dis);
 
-    bool solve(GridFunc<T2>& gf_phi, GridFunc<T2>& gf_rhs) override;
+    bool solve(GridFunc<T2>& gf_phi, const GridFunc<T2>& gf_rhs) override;
 
     short getNbSweeps() const override { return nb_sweeps_; }
     double getFinalResidual() const override { return final_residual_; }
diff --git a/src/pb/SolverPB.cc b/src/pb/SolverPB.cc
index 53252c26..756b49d8 100644
--- a/src/pb/SolverPB.cc
+++ b/src/pb/SolverPB.cc
@@ -7,7 +7,6 @@
 // This file is part of MGmol. For details, see https://github.com/llnl/mgmol.
 // Please also read this link https://github.com/llnl/mgmol/LICENSE
 
-// $Id: SolverPB.cc,v 1.13 2010/01/28 22:56:32 jeanluc Exp $
 #include "SolverPB.h"
 #include "Mgm.h"
 #include "PBh2.h"
@@ -85,7 +84,7 @@ bool SolverPB<T, T2>::solve(T2* phi, T2* rhs, T2* rhod, T2* vks, const char dis)
 }
 
 template <class T, typename T2>
-bool SolverPB<T, T2>::solve(GridFunc<T2>& gf_phi, GridFunc<T2>& gf_rhs,
+bool SolverPB<T, T2>::solve(GridFunc<T2>& gf_phi, const GridFunc<T2>& gf_rhs,
     GridFunc<T2>& gf_rhod, GridFunc<T2>& gf_vks)
 {
     oper_.init(gf_rhod);
@@ -101,7 +100,7 @@ bool SolverPB<T, T2>::solve(GridFunc<T2>& gf_phi, GridFunc<T2>& gf_rhs,
 }
 
 template <class T, typename T2>
-bool SolverPB<T, T2>::solve(GridFunc<T2>& gf_phi, GridFunc<T2>& gf_rhs)
+bool SolverPB<T, T2>::solve(GridFunc<T2>& gf_phi, const GridFunc<T2>& gf_rhs)
 {
     if (!oper_.initialized())
     {
diff --git a/src/pb/SolverPB.h b/src/pb/SolverPB.h
index a05e1d49..99eae0e3 100644
--- a/src/pb/SolverPB.h
+++ b/src/pb/SolverPB.h
@@ -64,9 +64,9 @@ class SolverPB : public Solver<T2>
     }
 
     bool solve(T2* phi, T2* rhs, T2* rhod, T2* vks, const char dis);
-    bool solve(GridFunc<T2>& gf_phi, GridFunc<T2>& gf_rhs,
+    bool solve(GridFunc<T2>& gf_phi, const GridFunc<T2>& gf_rhs,
         GridFunc<T2>& gf_rhod, GridFunc<T2>& gf_vks);
-    bool solve(GridFunc<T2>& gf_phi, GridFunc<T2>& gf_rhs) override;
+    bool solve(GridFunc<T2>& gf_phi, const GridFunc<T2>& gf_rhs) override;
 
     ~SolverPB() override{};
 
diff --git a/src/tools.cc b/src/tools.cc
index 4e1bd314..95790c59 100644
--- a/src/tools.cc
+++ b/src/tools.cc
@@ -228,7 +228,7 @@ void printWithTimeStamp(const std::string& string2print, std::ostream& os)
     int mpierr=MPI_Reduce(&s, &r, 1, MPI_INT, MPI_SUM, 0, mmpi.commGlobal());
     if( mpierr!=MPI_SUCCESS )
     {
-        cerr << " Error in MPI!!!" << endl;
+        cerr << " Error in MPI!!!" << std::endl;
         MPI_Abort(mmpi.commGlobal(),1);
     }
     if( r!=mmpi.size()*s && onpe0 )
@@ -388,3 +388,18 @@ void getkvector(const int index, const int kmax, int kvector[3])
     //    std::cout << " k=(" << kvector[0] << "," << kvector[1] << ","
     //              << kvector[2] << ")" << std::endl;
 }
+
+double getCharge(double* rho)
+{
+    Control& ct            = *(Control::instance());
+    Mesh* mymesh           = Mesh::instance();
+    const pb::Grid& mygrid = mymesh->grid();
+
+    double charge = mygrid.integralOverMesh(rho);
+
+    if (onpe0 && ct.verbose > 0)
+        std::cout << std::setprecision(8) << std::fixed << "Charge: " << charge
+                  << std::endl;
+
+    return charge;
+}
diff --git a/src/tools.h b/src/tools.h
index 4e8374da..be3d13e0 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -37,5 +37,6 @@ double minQuadPolynomial(const double e0, const double e1, const double de0,
 double minQuadPolynomialFrom3values(const double e0, const double e1,
     const double e12, const bool print_flag, std::ostream& os);
 void getkvector(const int index, const int kmax, int kvector[3]);
+double getCharge(double* rho);
 
 #endif
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index cc205167..6de34971 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -209,6 +209,8 @@ add_executable(testMGkernels
                ${CMAKE_SOURCE_DIR}/src/pb/FDkernels.cc
                ${CMAKE_SOURCE_DIR}/src/tools/Timer.cc
                ${CMAKE_SOURCE_DIR}/tests/ut_main.cc)
+add_executable(testIons
+               ${CMAKE_SOURCE_DIR}/tests/testIons.cc)
 add_executable(testGramMatrix
                ${CMAKE_SOURCE_DIR}/tests/testGramMatrix.cc
                ${CMAKE_SOURCE_DIR}/src/GramMatrix.cc
@@ -334,6 +336,10 @@ add_test(NAME testBatchLaph4
 add_test(NAME testtMGkernels
          COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
                  ${CMAKE_CURRENT_BINARY_DIR}/testMGkernels)
+add_test(NAME testIons
+         COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+                 ${CMAKE_CURRENT_BINARY_DIR}/testIons
+                 ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
 add_test(NAME testGramMatrix
          COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
                  ${CMAKE_CURRENT_BINARY_DIR}/testGramMatrix)
@@ -541,6 +547,7 @@ target_include_directories(testPowerDistMatrix PRIVATE ${Boost_INCLUDE_DIRS})
 target_include_directories(testDensityMatrix PRIVATE ${Boost_INCLUDE_DIRS})
 target_include_directories(testGramMatrix PRIVATE ${Boost_INCLUDE_DIRS})
 target_include_directories(testAndersonMix PRIVATE ${Boost_INCLUDE_DIRS})
+target_include_directories(testIons PRIVATE ${Boost_INCLUDE_DIRS} ${HDF5_INCLUDE_DIRS})
 
 target_link_libraries(testMPI PRIVATE MPI::MPI_CXX)
 target_link_libraries(testBlacsContext PRIVATE ${SCALAPACK_LIBRARIES}
@@ -550,6 +557,7 @@ target_link_libraries(testDirectionalReduce PRIVATE MPI::MPI_CXX)
 target_link_libraries(testEnergyAndForces PRIVATE mgmol_src)
 target_link_libraries(testWFEnergyAndForces PRIVATE mgmol_src)
 target_link_libraries(testDMandEnergyAndForces PRIVATE mgmol_src)
+target_link_libraries(testIons PRIVATE mgmol_src)
 
 if(${MAGMA_FOUND})
   target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES}
diff --git a/tests/WFEnergyAndForces/mgmol.cfg b/tests/WFEnergyAndForces/mgmol.cfg
index d543b626..12002703 100644
--- a/tests/WFEnergyAndForces/mgmol.cfg
+++ b/tests/WFEnergyAndForces/mgmol.cfg
@@ -18,7 +18,7 @@ pseudopotential=pseudo.H
 type=QUENCH
 [Quench]
 max_steps=50
-atol=1.e-9
+atol=1.e-8
 num_lin_iterations=2
 [Orbitals]
 initial_type=Gaussian
diff --git a/tests/testIons.cc b/tests/testIons.cc
new file mode 100644
index 00000000..30b2ad0c
--- /dev/null
+++ b/tests/testIons.cc
@@ -0,0 +1,104 @@
+#include "Control.h"
+#include "Ions.h"
+#include "MGmol_MPI.h"
+#include "Mesh.h"
+#include "Species.h"
+
+#include <random>
+
+int main(int argc, char** argv)
+{
+    int mpirc = MPI_Init(&argc, &argv);
+
+    MPI_Comm comm = MPI_COMM_WORLD;
+
+    int myrank;
+    MPI_Comm_rank(MPI_COMM_WORLD, &myrank);
+
+    MGmol_MPI::setup(comm, std::cout);
+    Control::setup(comm, false, 0.);
+
+    // create a domain [0.10.]^3
+    const double origin[3]  = { 0., 0., 0. };
+    const double ll         = 10.;
+    const double lattice[3] = { ll, ll, ll };
+    const unsigned ngpts[3] = { 32, 24, 20 };
+    short lap_type          = 0;
+
+    Mesh::setup(comm, ngpts, origin, lattice, lap_type);
+
+    const double h[3] = { ll / (double(ngpts[0])), ll / (double(ngpts[1])),
+        ll / (double(ngpts[2])) };
+
+    // random number generator
+    static std::random_device rd;
+    static std::mt19937 gen(rd());
+    static std::uniform_real_distribution<> dis(0.0, 1.0);
+
+    // create one species
+    Species sp(MPI_COMM_WORLD);
+
+    // read species info from pseudopotential file
+    std::string file_path = argv[1];
+    std::string filename(file_path + "/pseudo.C_ONCV_PBE_SG15");
+    std::cout << "Potential = " << filename << std::endl;
+
+    sp.read_1species(filename);
+    sp.set_dim_nl(h[0]);
+    sp.set_dim_l(h[0]);
+    sp.initPotentials('f', h[0], true);
+
+    // put species into a vector
+    std::vector<Species> vsp;
+    vsp.push_back(sp);
+
+    Ions ions(lattice, vsp);
+    ions.setupListIonsBoundaries(10000.);
+
+    double velocity[3] = { 0., 0., 0. };
+
+    // set "na" atoms coordinates and add them to "ions"
+    const int na = 10;
+    for (int i = 0; i < na; i++)
+    {
+        double x[3] = { origin[0] + lattice[0] * dis(gen),
+            origin[1] + lattice[1] * dis(gen),
+            origin[2] + lattice[2] * dis(gen) };
+        if (myrank == 0)
+            std::cout << "x,y,z = " << x[0] << ", " << x[1] << ", " << x[2]
+                      << std::endl;
+
+        // set all x to the values of PE0
+        MPI_Bcast(&x[0], 3, MPI_DOUBLE, 0, comm);
+
+        // make a name for atom based on species and order of reading in
+        std::string stri = std::to_string(i);
+        std::string aname("C" + stri);
+
+        ions.addIonToList(sp, aname, &x[0], velocity, false);
+    }
+
+    ions.setup();
+
+    std::vector<Ion*>& new_local_ions(ions.local_ions());
+
+    int nlocal = new_local_ions.size();
+    std::cout << "PE " << myrank << ", nlocal = " << nlocal << std::endl;
+
+    int ntotal = 0;
+    MPI_Allreduce(&nlocal, &ntotal, 1, MPI_INT, MPI_SUM, comm);
+    mpirc = MPI_Finalize();
+    if (mpirc != MPI_SUCCESS)
+    {
+        std::cerr << "MPI Finalize failed!!!" << std::endl;
+        return 1;
+    }
+
+    if (ntotal != na)
+    {
+        std::cout << "ntotal = " << ntotal << std::endl;
+        return 1;
+    }
+
+    return 0;
+}