diff --git a/src/DavidsonSolver.cc b/src/DavidsonSolver.cc
index f21cf7fa..c53e11c0 100644
--- a/src/DavidsonSolver.cc
+++ b/src/DavidsonSolver.cc
@@ -598,7 +598,6 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
 
                 rho_->computeRho(
                     orbitals, work_orbitals, dm11, dm12, dm21, dm22);
-                rho_->rescaleTotalCharge();
 
                 mgmol_strategy_->update_pot(vh_init, ions_);
 
@@ -663,7 +662,6 @@ int DavidsonSolver<OrbitalsType, MatrixType>::solve(
 
                 rho_->computeRho(
                     orbitals, work_orbitals, dm11, dm12, dm21, dm22);
-                rho_->rescaleTotalCharge();
             }
 
         } // inner iterations
diff --git a/src/MGmol.h b/src/MGmol.h
index 2a0d0783..b9d7c854 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -305,6 +305,7 @@ class MGmol : public MGmolInterface
     }
 
     OrbitalsType* loadOrbitalFromRestartFile(const std::string filename);
+
 #ifdef MGMOL_HAS_LIBROM
     int save_orbital_snapshot(std::string snapshot_dir, OrbitalsType& orbitals);
 #endif
diff --git a/src/MVPSolver.cc b/src/MVPSolver.cc
index 2f8060b7..8ac6a0be 100644
--- a/src/MVPSolver.cc
+++ b/src/MVPSolver.cc
@@ -52,18 +52,22 @@ MVPSolver<OrbitalsType, MatrixType>::MVPSolver(MPI_Comm comm, std::ostream& os,
       os_(os),
       n_inner_steps_(n_inner_steps),
       use_old_dm_(use_old_dm),
-      ions_(ions)
+      ions_(ions),
+      numst_(numst)
 {
-    history_length_ = 2;
-    eks_history_.resize(history_length_, 100000.);
+    Control& ct = *(Control::instance());
+    if (onpe0 && ct.verbose > 0)
+    {
+        os_ << "MVPSolver..." << std::endl;
+        if (use_old_dm_) os_ << "MVPSolver uses old DM..." << std::endl;
+    }
 
     rho_            = rho;
     energy_         = energy;
     electrostat_    = electrostat;
     mgmol_strategy_ = mgmol_strategy;
 
-    numst_ = numst;
-    work_  = new MatrixType("workMVP", numst_, numst_);
+    work_ = new MatrixType("workMVP", numst_, numst_);
 
     proj_mat_work_ = new ProjectedMatrices<MatrixType>(numst_, false, kbT);
     proj_mat_work_->setup(global_indexes);
@@ -123,9 +127,9 @@ void MVPSolver<OrbitalsType, MatrixType>::buildTarget_MVP(
 {
     target_tm_.start();
 
-    if (onpe0) std::cout << "buildTarget_MVP()..." << std::endl;
-
     Control& ct = *(Control::instance());
+    if (onpe0 && ct.verbose > 1)
+        std::cout << "buildTarget_MVP()..." << std::endl;
 
     proj_mat_work_->assignH(h11);
 
@@ -135,10 +139,8 @@ void MVPSolver<OrbitalsType, MatrixType>::buildTarget_MVP(
     //
     // compute target density matrix, with occupations>0 in numst only
     //
-    // if( onpe0 )os_<<"Compute XN..."<<endl;
     proj_mat_work_->setHB2H();
 
-    // if( onpe0 )os_<<"Compute DM..."<<endl;
     proj_mat_work_->updateDM(orbitals_index);
 
     target = proj_mat_work_->dm();
@@ -168,18 +170,12 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
         os_ << "---------------------------------------------------------------"
                "-"
             << std::endl;
-        os_ << "Update DM functions using MVP Solver..." << std::endl;
+        os_ << "Update DM using MVP Solver..." << std::endl;
         os_ << "---------------------------------------------------------------"
                "-"
             << std::endl;
     }
 
-    // step for numerical derivative
-
-    // double nel=orbitals.projMatrices()->getNel();
-    // if( onpe0 )
-    //    os_<<"NEW algorithm: Number of electrons at start = "<<nel<<endl;
-
     KBPsiMatrixSparse kbpsi(nullptr);
     kbpsi.setup(ions_);
 
@@ -227,11 +223,12 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                 energy_->saveVofRho();
             }
 
-            ProjectedMatricesInterface* current_proj_mat
-                = (inner_it == 0) ? orbitals.getProjMatrices() : proj_mat_work_;
+            ProjectedMatrices<MatrixType>* current_proj_mat
+                = (inner_it == 0)
+                      ? dynamic_cast<ProjectedMatrices<MatrixType>*>(
+                            orbitals.getProjMatrices())
+                      : proj_mat_work_;
 
-            double ts0;
-            double e0        = 0.;
             const int printE = (ct.verbose > 1) ? 1 : 0;
 
             // compute h11 for the current potential by adding local part to
@@ -239,49 +236,32 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
             MatrixType h11(h11_nl);
             mgmol_strategy_->addHlocal2matrix(orbitals, orbitals, h11);
 
+            current_proj_mat->assignH(h11);
+            current_proj_mat->setHB2H();
+
             if (inner_it == 0)
             {
-                // orbitals are new, so a few things need to recomputed
-                ProjectedMatrices<MatrixType>* projmatrices
-                    = dynamic_cast<ProjectedMatrices<MatrixType>*>(
-                        orbitals.getProjMatrices());
-
-                projmatrices->assignH(h11);
-                projmatrices->setHB2H();
-
                 if (use_old_dm_)
                 {
-                    dmInit = projmatrices->dm();
+                    dmInit = current_proj_mat->dm();
                 }
-
-                ts0 = evalEntropyMVP(projmatrices, true, os_);
-                e0  = energy_->evaluateTotal(
-                    ts0, projmatrices, orbitals, printE, os_);
             }
             else
             {
-
                 dmInit = proj_mat_work_->dm();
-
-                proj_mat_work_->assignH(h11);
-                proj_mat_work_->setHB2H();
-
-                ts0 = evalEntropyMVP(proj_mat_work_, true, os_);
-                e0  = energy_->evaluateTotal(
-                    ts0, proj_mat_work_, orbitals, printE, os_);
             }
 
-            // N x N target...
-            // proj_mat_work_->setHiterativeIndex(orbitals.getIterativeIndex(),
-            // pot.getIterativeIndex());
+            const double ts0 = evalEntropyMVP(current_proj_mat, true, os_);
+            const double e0  = energy_->evaluateTotal(
+                ts0, current_proj_mat, orbitals, printE, os_);
 
             MatrixType target("target", numst_, numst_);
-            MatrixType delta_dm("delta_dm", numst_, numst_);
 
             buildTarget_MVP(h11, s11, target);
 
             if (use_old_dm_ || inner_it > 0)
             {
+                MatrixType delta_dm("delta_dm", numst_, numst_);
                 delta_dm = target;
                 delta_dm -= dmInit;
 
@@ -291,18 +271,16 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                 // evaluate free energy at beta=1
                 //
                 if (onpe0 && ct.verbose > 2)
-                    std::cout << "Target energy..." << std::endl;
+                    std::cout << "MVP --- Target energy..." << std::endl;
                 proj_mat_work_->setDM(target, orbitals.getIterativeIndex());
                 proj_mat_work_->computeOccupationsFromDM();
                 if (ct.verbose > 2) current_proj_mat->printOccupations(os_);
-                double nel = proj_mat_work_->getNel();
+                const double nel = proj_mat_work_->getNel();
                 if (onpe0 && ct.verbose > 1)
-                    os_ << "Number of electrons at beta=1 : " << nel
+                    os_ << "MVP --- Number of electrons at beta=1 : " << nel
                         << std::endl;
 
-                // if( onpe0 )os_<<"Rho..."<<endl;
                 rho_->computeRho(orbitals, target);
-                rho_->rescaleTotalCharge();
 
                 mgmol_strategy_->update_pot(vh_init, ions_);
 
@@ -323,13 +301,13 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                     ts1, proj_mat_work_, orbitals, ct.verbose - 1, os_);
 
                 // line minimization
-                double beta
+                const double beta
                     = minQuadPolynomial(e0, e1, de0, (ct.verbose > 2), os_);
 
                 if (onpe0 && ct.verbose > 0)
                 {
                     os_ << std::setprecision(12);
-                    os_ << std::fixed << "Inner iteration " << inner_it
+                    os_ << std::fixed << "MVP inner iteration " << inner_it
                         << ", E0=" << e0 << ", E1=" << e1;
                     os_ << std::scientific << ", E0'=" << de0
                         << " -> beta=" << beta;
@@ -355,14 +333,13 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
                         os_ << "Number of electrons for interpolated DM = "
                             << pnel << std::endl;
                 }
-
                 // if( onpe0 )os_<<"Rho..."<<endl;
                 rho_->computeRho(orbitals, *work_);
-                rho_->rescaleTotalCharge();
             }
 
         } // inner iterations
 
+        // set DM to latest iteration value
         ProjectedMatrices<MatrixType>* projmatrices
             = dynamic_cast<ProjectedMatrices<MatrixType>*>(
                 orbitals.getProjMatrices());
@@ -375,11 +352,9 @@ int MVPSolver<OrbitalsType, MatrixType>::solve(OrbitalsType& orbitals)
     if (onpe0 && ct.verbose > 1)
     {
         os_ << "---------------------------------------------------------------"
-               "-"
             << std::endl;
         os_ << "End MVP Solver..." << std::endl;
         os_ << "---------------------------------------------------------------"
-               "-"
             << std::endl;
     }
     solve_tm_.stop();
diff --git a/src/MVPSolver.h b/src/MVPSolver.h
index a3c29d1d..8e7a7c70 100644
--- a/src/MVPSolver.h
+++ b/src/MVPSolver.h
@@ -24,27 +24,22 @@ template <class OrbitalsType, class MatrixType>
 class MVPSolver
 {
 private:
-    MPI_Comm comm_;
+    const MPI_Comm comm_;
     std::ostream& os_;
 
-    short n_inner_steps_;
+    const short n_inner_steps_;
 
-    bool use_old_dm_;
+    const bool use_old_dm_;
     Ions& ions_;
 
+    int numst_;
+
     Rho<OrbitalsType>* rho_;
     Energy<OrbitalsType>* energy_;
     Electrostatic* electrostat_;
 
-    int history_length_;
-    std::vector<double> eks_history_;
-
     MGmol<OrbitalsType>* mgmol_strategy_;
 
-    double de_old_;
-    double de_;
-
-    int numst_;
     MatrixType* work_;
     ProjectedMatrices<MatrixType>* proj_mat_work_;
 
diff --git a/src/Rho.cc b/src/Rho.cc
index 823e4478..b456f6c1 100644
--- a/src/Rho.cc
+++ b/src/Rho.cc
@@ -28,8 +28,6 @@ template <class OrbitalsType>
 Timer Rho<OrbitalsType>::compute_tm_("Rho::compute");
 template <class OrbitalsType>
 Timer Rho<OrbitalsType>::compute_blas_tm_("Rho::compute_usingBlas");
-// template <class OrbitalsType>
-// Timer Rho<OrbitalsType>::compute_offdiag_tm_("Rho::compute_offdiag");
 
 #ifdef HAVE_MAGMA
 template <typename ScalarType>
@@ -46,10 +44,6 @@ Rho<OrbitalsType>::Rho()
     myspin_         = mmpi.myspin();
     nspin_          = mmpi.nspin();
 
-    // default values for block sizes
-    //    block_functions_ = 8;
-    //    block_space_     = 256;
-
     Mesh* mymesh           = Mesh::instance();
     const pb::Grid& mygrid = mymesh->grid();
     np_                    = mygrid.size();
@@ -77,31 +71,6 @@ void Rho<OrbitalsType>::setup(const OrthoType orbitals_type,
     orbitals_indexes_ = orbitals_indexes;
 }
 
-template <class OrbitalsType>
-void Rho<OrbitalsType>::extrapolate()
-{
-    double minus = -1;
-    double two   = 2.;
-    if (rho_minus1_.empty())
-    {
-        rho_minus1_.resize(nspin_);
-        rho_minus1_[myspin_].resize(np_);
-        memcpy(&rho_minus1_[myspin_][0], &rho_[myspin_][0],
-            np_ * sizeof(RHODTYPE));
-        return;
-    }
-    RHODTYPE* tmp = new RHODTYPE[np_];
-    memcpy(tmp, &rho_[myspin_][0], np_ * sizeof(RHODTYPE));
-
-    LinearAlgebraUtils<MemorySpace::Host>::MPscal(np_, two, &rho_[myspin_][0]);
-    LinearAlgebraUtils<MemorySpace::Host>::MPaxpy(
-        np_, minus, &rho_minus1_[myspin_][0], &rho_[myspin_][0]);
-
-    memcpy(&rho_minus1_[myspin_][0], tmp, np_ * sizeof(RHODTYPE));
-
-    delete[] tmp;
-}
-
 template <class OrbitalsType>
 void Rho<OrbitalsType>::axpyRhoc(const double alpha, RHODTYPE* rhoc)
 {
@@ -144,8 +113,6 @@ void Rho<OrbitalsType>::update(OrbitalsType& current_orbitals)
 
     computeRho(current_orbitals);
 
-    rescaleTotalCharge();
-
     assert(iterative_index_ >= 0);
 
     update_tm_.stop();
@@ -234,68 +201,6 @@ void Rho<OrbitalsType>::rescaleTotalCharge()
 #endif
 }
 
-// template <class OrbitalsType>
-// int Rho<OrbitalsType>::setupSubdomainData(const int iloc,
-//    const vector<const T*>& vorbitals,
-//    const ProjectedMatricesInterface* const projmatrices,
-//    vector<MATDTYPE>& melements, vector<vector<const ORBDTYPE*>>& vmpsi)
-//{
-//    // printWithTimeStamp("Rho<OrbitalsType>::setupSubdomainData()...",cout);
-//
-//    const short norb = (short)vorbitals.size();
-//    vmpsi.resize(norb);
-//
-//    const vector<int>& loc_indexes(orbitals_indexes_[iloc]);
-//    const int n_colors = (int)loc_indexes.size();
-//
-//    if (n_colors == 0) return 0;
-//
-//    vector<int> mycolors;
-//    mycolors.reserve(n_colors);
-//    for (int icolor = 0; icolor < n_colors; icolor++)
-//        if (loc_indexes[icolor] != -1) mycolors.push_back(icolor);
-//    const int nmycolors = (int)mycolors.size();
-//
-//    for (short j = 0; j < norb; j++)
-//        vmpsi[j].resize(nmycolors);
-//
-//    short j = 0;
-//    for (typename vector<const T*>::const_iterator it = vorbitals.begin();
-//         it != vorbitals.end(); ++it)
-//    {
-//        for (int color = 0; color < nmycolors; color++)
-//        {
-//            vmpsi[j][color] = (*it)->getPsi(mycolors[color], iloc);
-//            assert(vmpsi[j][color] != NULL);
-//        }
-//        j++;
-//    }
-//
-//    SquareLocalMatrices<MATDTYPE,MemorySpace::Host>&
-//    localX(projmatrices->getLocalX()); const MATDTYPE* const localX_iloc =
-//    localX.getSubMatrix(iloc); melements.clear(); melements.resize(nmycolors *
-//    nmycolors);
-//
-//    for (int i = 0; i < nmycolors; i++)
-//    {
-//        const int icolor = mycolors[i];
-//        if (norb == 1)
-//        {
-//            melements[i * nmycolors + i]
-//                = localX_iloc[icolor + n_colors * icolor];
-//        }
-//        const int jmax = (norb == 1) ? i : nmycolors;
-//        for (int j = 0; j < jmax; j++)
-//        {
-//            int jcolor = mycolors[j];
-//            melements[j * nmycolors + i]
-//                = localX_iloc[icolor + n_colors * jcolor];
-//        }
-//    }
-//
-//    return nmycolors;
-//}
-
 template <class OrbitalsType>
 void Rho<OrbitalsType>::accumulateCharge(const double alpha, const short ix_max,
     const ORBDTYPE* const psii, const ORBDTYPE* const psij,
@@ -305,235 +210,6 @@ void Rho<OrbitalsType>::accumulateCharge(const double alpha, const short ix_max,
         plrho[ix] += (RHODTYPE)(alpha * (double)psii[ix] * (double)psij[ix]);
 }
 
-// template <class OrbitalsType>
-// void Rho<OrbitalsType>::computeRhoSubdomain(
-//    const int iloc_init, const int iloc_end, const OrbitalsType& orbitals)
-//{
-//    assert(orbitals_type_ == OrthoType::Eigenfunctions
-//        || orbitals_type_ == OrthoType::Nonorthogonal);
-//
-//    compute_tm_.start();
-//
-//    Mesh* mymesh = Mesh::instance();
-//
-//    const int loc_numpt = mymesh->locNumpt();
-//
-//    RHODTYPE* const prho = &rho_[myspin_][0];
-//
-//    vector<const T*> vorbitals;
-//    vorbitals.push_back(&orbitals);
-//
-//    const double nondiagfactor = 2.;
-//
-//    for (int iloc = iloc_init; iloc < iloc_end; iloc++)
-//    {
-//        const int istart = iloc * loc_numpt;
-//
-//        RHODTYPE* const lrho = &prho[istart];
-//
-//        vector<MATDTYPE> melements;
-//        vector<vector<const ORBDTYPE*>> vmpsi;
-//
-//        const int nmycolors = setupSubdomainData(
-//            iloc, vorbitals, orbitals.projMatrices(), melements, vmpsi);
-//        assert(vmpsi.size() == 1);
-//        vector<const ORBDTYPE*> mpsi = vmpsi[0];
-//
-//        const int nblocks_color = nmycolors / block_functions_;
-//        const int max_icolor    = (nmycolors % block_functions_ == 0)
-//                                   ? nmycolors
-//                                   : nblocks_color * block_functions_;
-//        const int missed_rows = nmycolors - max_icolor;
-//        //(*MPIdata::sout)<<"max_icolor="<<max_icolor<<endl;
-//        //(*MPIdata::sout)<<"Rho<OrbitalsType>::computeRhoSubdomain:
-//        //missed_rows="<<missed_rows<<endl;
-//        //(*MPIdata::sout)<<"nblocks_color="<<nblocks_color<<endl;
-//
-//#ifdef _OPENMP
-//#pragma omp parallel for
-//#endif
-//        for (int idx = 0; idx < loc_numpt; idx += block_space_)
-//        {
-//            const short ix_max = min(block_space_, loc_numpt - idx);
-//            // RHODTYPE* const plrho=lrho+idx;
-//
-//            // non-diagonal blocks
-//            for (int icolor = 0; icolor < max_icolor;
-//                 icolor += block_functions_)
-//                for (int jcolor = 0; jcolor < icolor;
-//                     jcolor += block_functions_)
-//                {
-//
-//                    nonOrthoRhoKernel(icolor, block_functions_, jcolor,
-//                        block_functions_, idx, ix_max, &melements[0],
-//                        nmycolors, mpsi, nondiagfactor, lrho);
-//                }
-//            // finish missing rows ()
-//            if (missed_rows)
-//                for (int jcolor = 0; jcolor < max_icolor;
-//                     jcolor += block_functions_)
-//                {
-//                    nonOrthoRhoKernel(max_icolor, missed_rows, jcolor,
-//                        block_functions_, idx, ix_max, &melements[0],
-//                        nmycolors, mpsi, nondiagfactor, lrho);
-//                }
-//
-//            // diagonal blocks (jcolor=icolor)
-//            for (int icolor = 0; icolor < max_icolor;
-//                 icolor += block_functions_)
-//            {
-//                nonOrthoRhoKernelDiagonalBlock(icolor, block_functions_, idx,
-//                    ix_max, &melements[0], nmycolors, mpsi, lrho);
-//            }
-//            // finish missing rows () for diagonal blocks
-//            if (missed_rows)
-//            {
-//                nonOrthoRhoKernelDiagonalBlock(max_icolor, missed_rows, idx,
-//                    ix_max, &melements[0], nmycolors, mpsi, lrho);
-//            }
-//        }
-//    }
-//
-//    compute_tm_.stop();
-//}
-
-// template <class OrbitalsType>
-// void Rho<OrbitalsType>::computeRhoSubdomainOffDiagBlock(const int iloc_init,
-//    const int iloc_end, const vector<const T*>& vorbitals,
-//    const ProjectedMatricesInterface* const projmatrices)
-//{
-//    assert(orbitals_type_ == OrthoType::Eigenfunctions
-//        || orbitals_type_ == OrthoType::Nonorthogonal);
-//    assert(vorbitals.size() == 2);
-//
-//    compute_offdiag_tm_.start();
-//
-//    //
-//    printWithTimeStamp("Rho<OrbitalsType>::computeRhoSubdomainOffDiagBlock()...",cout);
-//
-//    Mesh* mymesh = Mesh::instance();
-//
-//    const int loc_numpt = mymesh->locNumpt();
-//
-//    RHODTYPE* const prho = &rho_[myspin_][0];
-//
-//    for (int iloc = iloc_init; iloc < iloc_end; iloc++)
-//    {
-//        const int istart = iloc * loc_numpt;
-//
-//        RHODTYPE* const lrho = &prho[istart];
-//
-//        vector<MATDTYPE> melements;
-//        vector<vector<const ORBDTYPE*>> vmpsi;
-//
-//        const int nmycolors = setupSubdomainData(
-//            iloc, vorbitals, projmatrices, melements, vmpsi);
-//        const int nblocks_color = nmycolors / block_functions_;
-//        const int max_icolor    = (nmycolors % block_functions_ == 0)
-//                                   ? nmycolors
-//                                   : nblocks_color * block_functions_;
-//        const int missed_rows = nmycolors - max_icolor;
-//        //(*MPIdata::sout)<<"max_icolor="<<max_icolor<<endl;
-//        //(*MPIdata::sout)<<"Rho<OrbitalsType>::computeRhoSubdomainOffDiagBlock:
-//        //missed_rows="<<missed_rows<<endl;
-//        //(*MPIdata::sout)<<"nblocks_color="<<nblocks_color<<endl;
-//
-//        for (int idx = 0; idx < loc_numpt; idx += block_space_)
-//        {
-//            const short ix_max    = min(block_space_, loc_numpt - idx);
-//            RHODTYPE* const plrho = lrho + idx;
-//
-//            for (int istart = 0; istart < max_icolor;
-//                 istart += block_functions_)
-//                for (int jstart = 0; jstart < max_icolor;
-//                     jstart += block_functions_)
-//                {
-//                    for (short jcolor = jstart;
-//                         jcolor < jstart + block_functions_; jcolor++)
-//                    {
-//                        const int jld              = jcolor * nmycolors;
-//                        const ORBDTYPE* const psij = &vmpsi[1][jcolor][idx];
-//                        for (short icolor = istart;
-//                             icolor < istart + block_functions_; icolor++)
-//                        {
-//                            const ORBDTYPE* const psii =
-//                            &vmpsi[0][icolor][idx]; const double alpha
-//                                = (double)melements[jld + icolor];
-//                            accumulateCharge(alpha, ix_max, psii, psij,
-//                            plrho);
-//                        }
-//                    }
-//                }
-//            // finish missing rows and cols
-//            if (missed_rows)
-//            {
-//                int icolor = max_icolor;
-//                int imax   = nmycolors % block_functions_;
-//                for (int jcolor = 0; jcolor < max_icolor;
-//                     jcolor += block_functions_)
-//                {
-//                    for (short j = 0; j < block_functions_; j++)
-//                    {
-//                        const int jstart = (jcolor + j) * nmycolors + icolor;
-//                        const ORBDTYPE* const psij = &vmpsi[1][jcolor +
-//                        j][idx]; for (short i = 0; i < imax; i++)
-//                        {
-//                            const ORBDTYPE* const psii
-//                                = &vmpsi[0][icolor + i][idx];
-//                            const double alpha = (double)melements[jstart +
-//                            i]; accumulateCharge(alpha, ix_max, psii, psij,
-//                            plrho);
-//                        }
-//                    }
-//                }
-//            }
-//            if (missed_rows)
-//            {
-//                int jcolor = max_icolor;
-//                int jmax   = nmycolors % block_functions_;
-//                for (int icolor = 0; icolor < max_icolor;
-//                     icolor += block_functions_)
-//                {
-//                    for (short j = 0; j < jmax; j++)
-//                    {
-//                        const int jstart = (jcolor + j) * nmycolors + icolor;
-//                        const ORBDTYPE* const psij = &vmpsi[1][jcolor +
-//                        j][idx]; for (short i = 0; i < block_functions_; i++)
-//                        {
-//                            const ORBDTYPE* const psii
-//                                = &vmpsi[0][icolor + i][idx];
-//                            const double alpha = (double)melements[jstart +
-//                            i]; accumulateCharge(alpha, ix_max, psii, psij,
-//                            plrho);
-//                        }
-//                    }
-//                }
-//            }
-//            if (missed_rows)
-//            {
-//                int icolor = max_icolor;
-//                int imax   = nmycolors % block_functions_;
-//                int jcolor = max_icolor;
-//                int jmax   = nmycolors % block_functions_;
-//                for (short j = 0; j < jmax; j++)
-//                {
-//                    const int jstart = (jcolor + j) * nmycolors + icolor;
-//                    const ORBDTYPE* const psij = &vmpsi[1][jcolor + j][idx];
-//                    for (short i = 0; i < imax; i++)
-//                    {
-//                        const ORBDTYPE* const psii = &vmpsi[0][icolor +
-//                        i][idx]; const double alpha = (double)melements[jstart
-//                        + i]; accumulateCharge(alpha, ix_max, psii, psij,
-//                        plrho);
-//                    }
-//                }
-//            }
-//        }
-//    }
-//
-//    compute_offdiag_tm_.stop();
-//}
-
 template <class OrbitalsType>
 void Rho<OrbitalsType>::computeRhoSubdomain(const int iloc_init,
     const int iloc_end, const OrbitalsType& orbitals,
@@ -586,6 +262,16 @@ void Rho<OrbitalsType>::computeRho(OrbitalsType& orbitals)
     computeRho(orbitals, proj_matrices);
 }
 
+template <class OrbitalsType>
+void Rho<OrbitalsType>::resetValues()
+{
+    Mesh* mymesh        = Mesh::instance();
+    const int subdivx   = mymesh->subdivx();
+    const int loc_numpt = mymesh->locNumpt();
+
+    memset(&rho_[myspin_][0], 0, subdivx * loc_numpt * sizeof(RHODTYPE));
+}
+
 template <class OrbitalsType>
 void Rho<OrbitalsType>::computeRho(
     OrbitalsType& orbitals, ProjectedMatricesInterface& proj_matrices)
@@ -593,12 +279,11 @@ void Rho<OrbitalsType>::computeRho(
     assert(rho_.size() > 0);
     assert(rho_[myspin_].size() > 0);
 
-    Mesh* mymesh        = Mesh::instance();
-    const int subdivx   = mymesh->subdivx();
-    const int loc_numpt = mymesh->locNumpt();
-    Control& ct         = *(Control::instance());
+    Mesh* mymesh      = Mesh::instance();
+    const int subdivx = mymesh->subdivx();
+    Control& ct       = *(Control::instance());
 
-    memset(&rho_[myspin_][0], 0, subdivx * loc_numpt * sizeof(RHODTYPE));
+    resetValues();
 
     if (orbitals_type_ == OrthoType::Eigenfunctions
         || (orbitals_type_ == OrthoType::Orthonormal && ct.fullyOccupied()))
@@ -611,7 +296,6 @@ void Rho<OrbitalsType>::computeRho(
     {
         proj_matrices.updateSubMatX();
 
-        // if (dynamic_cast<LocGridOrbitals*>(&orbitals)) but it
         if (std::is_same<OrbitalsType, LocGridOrbitals>::value)
         {
             SquareLocalMatrices<MATDTYPE, memory_space_type>& localX(
@@ -625,6 +309,8 @@ void Rho<OrbitalsType>::computeRho(
     }
 
     gatherSpin();
+
+    rescaleTotalCharge();
 }
 
 template <class OrbitalsType>
@@ -635,11 +321,10 @@ void Rho<OrbitalsType>::computeRho(OrbitalsType& orbitals1,
 {
     assert(orbitals_type_ == OrthoType::Nonorthogonal);
 
-    Mesh* mymesh        = Mesh::instance();
-    const int subdivx   = mymesh->subdivx();
-    const int loc_numpt = mymesh->locNumpt();
+    Mesh* mymesh      = Mesh::instance();
+    const int subdivx = mymesh->subdivx();
 
-    memset(&rho_[myspin_][0], 0, subdivx * loc_numpt * sizeof(RHODTYPE));
+    resetValues();
 
     // 11 diagonal block
     ProjectedMatrices<MatrixType>* projmatrices1
@@ -655,18 +340,14 @@ void Rho<OrbitalsType>::computeRho(OrbitalsType& orbitals1,
     projmatrices2->updateSubMatX(dm22);
     computeRhoSubdomainUsingBlas3(0, subdivx, orbitals2);
 
-    // non diagonal blocks
-    // vector<const T*> vorbitals;
-    // vorbitals.push_back(&orbitals1);
-    // vorbitals.push_back(&orbitals2);
-
     MatrixType dm(dm12);
     dm.scal(2.); // use symmetry to reduce work
     projmatrices1->updateSubMatX(dm);
-    //    computeRhoSubdomainOffDiagBlock(0, subdivx, vorbitals, projmatrices1);
     computeRhoSubdomainUsingBlas3(0, subdivx, orbitals1, orbitals2);
 
     gatherSpin();
+
+    rescaleTotalCharge();
 }
 
 template <class OrbitalsType>
@@ -769,11 +450,10 @@ void Rho<OrbitalsType>::computeRho(OrbitalsType& orbitals, const MatrixType& dm)
 
     iterative_index_++;
 
-    Mesh* mymesh        = Mesh::instance();
-    const int subdivx   = mymesh->subdivx();
-    const int loc_numpt = mymesh->locNumpt();
+    Mesh* mymesh      = Mesh::instance();
+    const int subdivx = mymesh->subdivx();
 
-    memset(&rho_[myspin_][0], 0, subdivx * loc_numpt * sizeof(RHODTYPE));
+    resetValues();
 
     ProjectedMatrices<MatrixType>* projmatrices
         = dynamic_cast<ProjectedMatrices<MatrixType>*>(
@@ -783,6 +463,8 @@ void Rho<OrbitalsType>::computeRho(OrbitalsType& orbitals, const MatrixType& dm)
     computeRhoSubdomainUsingBlas3(0, subdivx, orbitals);
 
     gatherSpin();
+
+    rescaleTotalCharge();
 }
 
 template <class OrbitalsType>
@@ -879,7 +561,6 @@ void Rho<OrbitalsType>::printTimers(std::ostream& os)
     update_tm_.print(os);
     compute_tm_.print(os);
     compute_blas_tm_.print(os);
-    // compute_offdiag_tm_.print(os);
 }
 
 template class Rho<LocGridOrbitals>;
diff --git a/src/Rho.h b/src/Rho.h
index 31130c7e..ea5069da 100644
--- a/src/Rho.h
+++ b/src/Rho.h
@@ -35,24 +35,13 @@ class Rho
 
     int verbosity_level_;
 
-    // parameters for "blocking" in loops over functions and space
-    // int block_functions_;
-    // int block_space_;
-
     static Timer update_tm_;
     static Timer compute_tm_;
     static Timer compute_blas_tm_;
-    // static Timer compute_offdiag_tm_;
 
     double computeTotalCharge();
-    // void computeRhoSubdomain(const int iloc_init, const int iloc_end,
-    //    const T& orbitals);
     void computeRhoSubdomain(const int iloc_init, const int iloc_end,
         const OrbitalsType& orbitals, const std::vector<double>& occ);
-    // void computeRhoSubdomainOffDiagBlock(const int iloc_init,
-    //    const int iloc_end,
-    //    const std::vector<const T*>& vorbitals,
-    //    const ProjectedMatricesInterface* const);
     void computeRhoSubdomainUsingBlas3(const int iloc_init, const int iloc_end,
         const OrbitalsType& orbitals1, const OrbitalsType& orbitals2);
     void computeRhoSubdomainUsingBlas3(
@@ -64,11 +53,6 @@ class Rho
     void accumulateCharge(const double alpha, const short ix_max,
         const ORBDTYPE* const psii, const ORBDTYPE* const psij,
         RHODTYPE* const plrho);
-    // int setupSubdomainData(const int iloc,
-    //    const std::vector<const T*>& vorbitals,
-    //    const ProjectedMatricesInterface* const projmatrices,
-    //    std::vector<MATDTYPE>& melements,
-    //    std::vector<std::vector<const ORBDTYPE*>>& mpsi);
 
     void computeRho(OrbitalsType& orbitals);
     void computeRho(
@@ -76,15 +60,17 @@ class Rho
 
     void gatherSpin();
 
+    void resetValues();
+
+    void rescaleTotalCharge();
+
 public:
     // electronic density on grid
     std::vector<std::vector<RHODTYPE>> rho_;
-    std::vector<std::vector<RHODTYPE>> rho_minus1_;
 
     Rho();
     ~Rho(){};
 
-    void rescaleTotalCharge();
     void setup(
         const OrthoType orbitals_type, const std::vector<std::vector<int>>&);
     void setVerbosityLevel(const int vlevel) { verbosity_level_ = vlevel; }
@@ -105,18 +91,11 @@ class Rho
     void initUniform();
     int getIterativeIndex() const { return iterative_index_; }
     int readRestart(HDFrestart& file);
-    void extrapolate();
     void axpyRhoc(const double alpha, RHODTYPE* rhoc);
 
     template <typename T2>
     double dotWithRho(const T2* const func) const;
 
-    // void setupBlockSizes(const int block_functions, const int block_space)
-    //{
-    //    block_functions_ = block_functions;
-    //    block_space_     = block_space;
-    //}
-
     static void printTimers(std::ostream& os);
 };