diff --git a/src/Electrostatic.cc b/src/Electrostatic.cc
index 6ac4fef3..aba09d07 100644
--- a/src/Electrostatic.cc
+++ b/src/Electrostatic.cc
@@ -65,7 +65,7 @@ Electrostatic::~Electrostatic()
     if (grhoc_ != nullptr) delete grhoc_;
 }
 
-void Electrostatic::setupInitialVh(const POTDTYPE* const vh_init)
+void Electrostatic::setupInitialVh(const std::vector<POTDTYPE>& vh_init)
 {
     poisson_solver_->set_vh(vh_init);
 
@@ -153,7 +153,7 @@ void Electrostatic::setupPB(
 
     // initialize vh with last trial solution
     pb::GridFunc<POTDTYPE> gf_vh(*pbGrid_, bc_[0], bc_[1], bc_[2]);
-    gf_vh.assign(pot.vh_rho());
+    gf_vh.assign((pot.vh_rho()).data());
     poisson_solver_->set_vh(gf_vh);
 }
 
diff --git a/src/Electrostatic.h b/src/Electrostatic.h
index 072ce0cc..9beee899 100644
--- a/src/Electrostatic.h
+++ b/src/Electrostatic.h
@@ -58,7 +58,7 @@ class Electrostatic
     template <class T>
     void computeVh(const pb::GridFunc<POTDTYPE>& vhinit, const Ions& ions,
         Rho<T>& rho, Potentials& pot);
-    void setupInitialVh(const POTDTYPE* const);
+    void setupInitialVh(const std::vector<POTDTYPE>&);
     void setupInitialVh(const pb::GridFunc<POTDTYPE>&);
     template <class T>
     void computeVhRho(Rho<T>& rho);
diff --git a/src/Forces.cc b/src/Forces.cc
index 9ab19f6a..04e6ffd7 100644
--- a/src/Forces.cc
+++ b/src/Forces.cc
@@ -272,9 +272,10 @@ void Forces<T>::get_loc_proj(RHODTYPE* rho,
     const int numpt = mymesh->numpt();
 
     Potentials& pot = hamiltonian_->potential();
+    const std::vector<POTDTYPE>& vh_rho(pot.vh_rho());
     for (int idx = 0; idx < numpt; idx++)
     {
-        const double vhrho = pot.vh_rho(idx);
+        const double vhrho = vh_rho[idx];
         for (short dir = 0; dir < 3; dir++)
         {
             double* lproj = &(loc_proj[dir * NPTS]);
diff --git a/src/MGmol.cc b/src/MGmol.cc
index 681efacc..649f2c39 100644
--- a/src/MGmol.cc
+++ b/src/MGmol.cc
@@ -319,8 +319,6 @@ int MGmol<OrbitalsType>::initial()
     mmpi.barrier();
     if (ct.verbose > 0) current_orbitals_->printChromaticNumber(os_);
 
-    pot.initBackground(*ions_);
-
     // Random initialization of the wavefunctions
     if (ct.restart_info <= 2)
     {
@@ -834,23 +832,9 @@ void MGmol<OrbitalsType>::initNuc(Ions& ions)
 
     Potentials& pot = hamiltonian_->potential();
 
+    // initialize poentials based on ionic positions and their species
     pot.initialize(ions);
 
-    // Check compensating charges
-    double comp_rho = getCharge(pot.rho_comp());
-
-    if (onpe0 && ct.verbose > 1)
-    {
-        os_ << std::setprecision(8) << std::fixed
-            << " Charge of rhoc: " << comp_rho << std::endl;
-    }
-
-#if 1
-    pot.rescaleRhoComp();
-#endif
-
-    pot.addBackgroundToRhoComp();
-
     electrostat_->setupRhoc(pot.rho_comp());
 
     if (onpe0 && ct.verbose > 3) os_ << " initNuc done" << std::endl;
diff --git a/src/Poisson.h b/src/Poisson.h
index 87107542..ba5daa39 100644
--- a/src/Poisson.h
+++ b/src/Poisson.h
@@ -71,7 +71,10 @@ class Poisson : public PoissonInterface
 
     virtual void set_rhod(pb::GridFunc<RHODTYPE>* /*rhod*/){};
     void set_vh(const pb::GridFunc<POTDTYPE>& vh) { (*vh_) = vh; };
-    void set_vh(const POTDTYPE* const vh) { vh_->assign(vh, 'd'); };
+    void set_vh(const std::vector<POTDTYPE>& vh)
+    {
+        vh_->assign(vh.data(), 'd');
+    };
     void resetVh() { vh_->resetData(); }
     void set_vepsilon(const POTDTYPE* const vepsilon)
     {
diff --git a/src/Potentials.cc b/src/Potentials.cc
index abe277f1..a0ff3f1a 100644
--- a/src/Potentials.cc
+++ b/src/Potentials.cc
@@ -25,7 +25,6 @@
 #include "mputils.h"
 
 #include <fstream>
-using namespace std;
 
 // unit conversion factor Ha -> Ry
 const double ha2ry = 2.;
@@ -92,7 +91,7 @@ void Potentials::initWithVnuc()
 {
     assert(size_ > 0);
     if (verbosity_level_ > 2 && onpe0)
-        (*MPIdata::sout) << "Potentials::initWithVnuc()" << endl;
+        (*MPIdata::sout) << "Potentials::initWithVnuc()" << std::endl;
     itindex_vxc_ = 0;
     itindex_vh_  = 0;
     int ione     = 1;
@@ -122,7 +121,8 @@ double Potentials::min() const
     return vmin;
 }
 
-void Potentials::evalNormDeltaVtotRho(const vector<vector<RHODTYPE>>& rho)
+void Potentials::evalNormDeltaVtotRho(
+    const std::vector<std::vector<RHODTYPE>>& rho)
 {
     Mesh* mymesh           = Mesh::instance();
     const pb::Grid& mygrid = mymesh->grid();
@@ -141,14 +141,14 @@ void Potentials::evalNormDeltaVtotRho(const vector<vector<RHODTYPE>>& rho)
     mmpi.allreduce(&scf_dvrho_, 1, MPI_SUM);
 }
 
-double Potentials::update(const vector<vector<RHODTYPE>>& rho)
+double Potentials::update(const std::vector<std::vector<RHODTYPE>>& rho)
 {
     assert(itindex_vxc_ >= 0);
     assert(itindex_vh_ >= 0);
     assert(itindex_vxc_ == itindex_vh_);
 
     if (verbosity_level_ > 2 && onpe0)
-        (*MPIdata::sout) << "Potentials::update(rho)" << endl;
+        (*MPIdata::sout) << "Potentials::update(rho)" << std::endl;
     int ione                 = 1;
     Mesh* mymesh             = Mesh::instance();
     const pb::PEenv& myPEenv = mymesh->peenv();
@@ -185,7 +185,7 @@ double Potentials::update(const vector<vector<RHODTYPE>>& rho)
         = MPI_Allreduce(&dvdot, &sum, 1, MPI_DOUBLE, MPI_SUM, myPEenv.comm());
     if (rc != MPI_SUCCESS)
     {
-        cout << "MPI_Allreduce double sum failed!!!" << endl;
+        std::cout << "MPI_Allreduce double sum failed!!!" << std::endl;
         MPI_Abort(myPEenv.comm(), 2);
     }
     dvdot = sum;
@@ -202,7 +202,7 @@ void Potentials::update(const double mix)
     assert(itindex_vxc_ == itindex_vh_);
 
 #ifdef DEBUG
-    if (onpe0) (*MPIdata::sout) << "Potentials::update(mix)" << endl;
+    if (onpe0) (*MPIdata::sout) << "Potentials::update(mix)" << std::endl;
 #endif
     //    int ione=1;
     double potmix = mix;
@@ -210,13 +210,13 @@ void Potentials::update(const double mix)
         size_, potmix, &dv_[0], &vtot_[0]);
 }
 
-double Potentials::delta_v(const vector<vector<RHODTYPE>>& rho)
+double Potentials::delta_v(const std::vector<std::vector<RHODTYPE>>& rho)
 {
     assert(itindex_vxc_ == itindex_vh_);
     assert(size_ > 0);
 
     if (verbosity_level_ > 2 && onpe0)
-        (*MPIdata::sout) << "Potentials::delta_v()" << endl;
+        (*MPIdata::sout) << "Potentials::delta_v()" << std::endl;
 
     int ione                 = 1;
     Mesh* mymesh             = Mesh::instance();
@@ -253,7 +253,7 @@ double Potentials::delta_v(const vector<vector<RHODTYPE>>& rho)
         = MPI_Allreduce(&dvdot, &sum, 1, MPI_DOUBLE, MPI_SUM, myPEenv.comm());
     if (rc != MPI_SUCCESS)
     {
-        cout << "MPI_Allreduce double sum failed!!!" << endl;
+        std::cout << "MPI_Allreduce double sum failed!!!" << std::endl;
         MPI_Abort(myPEenv.comm(), 2);
     }
     dvdot = sum;
@@ -266,7 +266,7 @@ double Potentials::delta_v(const vector<vector<RHODTYPE>>& rho)
 }
 
 // in Ry
-void Potentials::getVofRho(vector<POTDTYPE>& vrho) const
+void Potentials::getVofRho(std::vector<POTDTYPE>& vrho) const
 {
     vrho.resize(size_);
     int ione        = 1;
@@ -283,7 +283,7 @@ void Potentials::getVofRho(vector<POTDTYPE>& vrho) const
 // type:
 // 2->text
 // 3->binary
-void Potentials::readExternalPot(const string filename, const char type)
+void Potentials::readExternalPot(const std::string filename, const char type)
 {
     assert(type == 2 || type == 3);
 
@@ -304,27 +304,27 @@ void Potentials::readExternalPot(const string filename, const char type)
         {
             from->seekg(0, ios::end);
             const int length = from->tellg();
-            (*MPIdata::sout) << "Length file = " << length << endl;
+            (*MPIdata::sout) << "Length file = " << length << std::endl;
             from->seekg(0, ios::beg);
             if (length <= 0)
             {
                 (*MPIdata::serr)
-                    << "ERROR Potential: file length <=0!!!!" << endl;
+                    << "ERROR Potential: file length <=0!!!!" << std::endl;
                 mmpi.abort();
             }
         }
     }
     if (!from)
     {
-        (*MPIdata::serr) << " Cannot open file " << filename << endl;
+        (*MPIdata::serr) << " Cannot open file " << filename << std::endl;
         mmpi.abort();
     }
     if (onpe0)
     {
         (*MPIdata::sout) << "Potentials::read_ExternalPot(), filename="
-                         << filename << endl;
-        if (type == 2) (*MPIdata::sout) << "text file..." << endl;
-        if (type == 3) (*MPIdata::sout) << "binary file..." << endl;
+                         << filename << std::endl;
+        if (type == 2) (*MPIdata::sout) << "text file..." << std::endl;
+        if (type == 3) (*MPIdata::sout) << "binary file..." << std::endl;
     }
 
     // read origin and end of cell (to check compatibility)
@@ -348,29 +348,32 @@ void Potentials::readExternalPot(const string filename, const char type)
     if (onpe0)
         for (short d = 0; d < 3; d++)
         {
-            (*MPIdata::sout) << setprecision(8);
+            (*MPIdata::sout) << std::setprecision(8);
             if (fabs(origin[d] - mygrid.origin(d)) > 1.e-3)
             {
                 (*MPIdata::serr)
                     << "ERROR Potential: Incompatible cell origin in direction "
-                    << d << endl;
-                (*MPIdata::serr) << "Potential origin=" << origin[d] << endl;
-                (*MPIdata::serr) << "MGmol origin=" << mygrid.origin(d) << endl;
+                    << d << std::endl;
+                (*MPIdata::serr)
+                    << "Potential origin=" << origin[d] << std::endl;
+                (*MPIdata::serr)
+                    << "MGmol origin=" << mygrid.origin(d) << std::endl;
                 (*MPIdata::serr)
                     << "Difference=" << fabs(origin[d] - mygrid.origin(d))
-                    << endl;
+                    << std::endl;
                 mmpi.abort();
             }
             if (fabs(ll[d] - mygrid.ll(d)) > 1.e-3)
             {
                 (*MPIdata::serr) << "ERROR Potential: Incompatible cell "
                                     "dimension in direction "
-                                 << d << endl;
-                (*MPIdata::serr) << "Potential cell end=" << end[d] << endl;
+                                 << d << std::endl;
+                (*MPIdata::serr)
+                    << "Potential cell end=" << end[d] << std::endl;
                 (*MPIdata::serr)
-                    << "Potential cell dimension=" << ll[d] << endl;
+                    << "Potential cell dimension=" << ll[d] << std::endl;
                 (*MPIdata::serr)
-                    << "MGmol cell dimension=" << mygrid.ll(d) << endl;
+                    << "MGmol cell dimension=" << mygrid.ll(d) << std::endl;
                 mmpi.abort();
             }
         }
@@ -392,9 +395,9 @@ void Potentials::readExternalPot(const string filename, const char type)
         if (nxyz[i] != gdim_[i])
         {
             (*MPIdata::serr) << "Potentials::read_ExternalPot(): dimension "
-                             << i << " incompatible with Grid!!!" << endl;
+                             << i << " incompatible with Grid!!!" << std::endl;
             (*MPIdata::serr)
-                << "n=" << nxyz[i] << ", gdim_=" << gdim_[i] << endl;
+                << "n=" << nxyz[i] << ", gdim_=" << gdim_[i] << std::endl;
             mmpi.abort();
         }
 
@@ -454,7 +457,7 @@ void Potentials::readExternalPot(const string filename, const char type)
 
     if (type == 3)
     {
-        vector<float> tmp(dim_[2]);
+        std::vector<float> tmp(dim_[2]);
         for (int i = 0; i < dim_[0]; i++)
         {
             // advance (start-file_index) positions
@@ -512,7 +515,8 @@ void Potentials::getGradVext(const double r[3], double dfdr[3]) const
     vext_tricubic_->getGradient(r, dfdr, comm);
 }
 
-void Potentials::getValVext(const vector<double>& r, vector<double>& val) const
+void Potentials::getValVext(
+    const std::vector<double>& r, std::vector<double>& val) const
 {
     assert(vext_tricubic_ != NULL);
 
@@ -522,21 +526,22 @@ void Potentials::getValVext(const vector<double>& r, vector<double>& val) const
 }
 #endif
 
-void Potentials::readAll(vector<Species>& sp)
+void Potentials::readAll(std::vector<Species>& sp)
 {
     assert(sp.size() <= pot_filenames_.size());
 
     if (verbosity_level_ > 2 && onpe0)
         (*MPIdata::sout) << "Potentials::readAll() for " << pot_types_.size()
-                         << " potentials" << endl;
+                         << " potentials" << std::endl;
     Mesh* mymesh           = Mesh::instance();
     const pb::Grid& mygrid = mymesh->grid();
     double hmin            = mygrid.hmin();
     if (verbosity_level_ > 2 && onpe0)
-        (*MPIdata::sout) << "hmin= " << hmin << endl;
+        (*MPIdata::sout) << "hmin= " << hmin << std::endl;
 
-    vector<string>::const_iterator it_filename = pot_filenames_.begin();
-    int isp                                    = 0;
+    std::vector<std::string>::const_iterator it_filename
+        = pot_filenames_.begin();
+    int isp = 0;
     while (it_filename != pot_filenames_.end())
     {
         if (pot_types_[isp] == 'n' || pot_types_[isp] == 's'
@@ -556,7 +561,7 @@ void Potentials::readAll(vector<Species>& sp)
 #else
             (*MPIdata::sout)
                 << "ERROR: cannot read external potential "
-                << " -> need to compile with Tricubic library" << endl;
+                << " -> need to compile with Tricubic library" << std::endl;
 #endif
         }
         it_filename++;
@@ -594,11 +599,6 @@ void Potentials::axpVcompToVh(const double alpha)
         size_, alpha, &v_comp_[0], &vh_rho_[0]);
 }
 
-void Potentials::axpVcomp(POTDTYPE* v, const double alpha)
-{
-    LinearAlgebraUtils<MemorySpace::Host>::MPaxpy(size_, alpha, &v_comp_[0], v);
-}
-
 void Potentials::backupVh()
 {
     memcpy(vh_rho_backup_.data(), vh_rho_.data(), size_ * sizeof(POTDTYPE));
@@ -792,9 +792,12 @@ void Potentials::initialize(Ions& ions)
     const pb::Grid& mygrid = mymesh->grid();
     const int numpt        = mygrid.size();
 
-    memset(&v_comp_[0], 0, numpt * sizeof(POTDTYPE));
-    memset(&rho_comp_[0], 0, numpt * sizeof(RHODTYPE));
-    memset(&v_nuc_[0], 0, numpt * sizeof(RHODTYPE));
+    memset(v_comp_.data(), 0, numpt * sizeof(POTDTYPE));
+    memset(rho_comp_.data(), 0, numpt * sizeof(RHODTYPE));
+    memset(v_nuc_.data(), 0, numpt * sizeof(RHODTYPE));
+
+    // Count up the total ionic charge
+    ionic_charge_ = ions.computeIonicCharge();
 
     char flag_filter = pot_type(0);
 
@@ -805,6 +808,7 @@ void Potentials::initialize(Ions& ions)
 
         Vector3D position(ion->position(0), ion->position(1), ion->position(2));
 
+        // initialize rho_comp_, v_comp_, v_nuc_
         if (flag_filter == 's')
         {
             const int sampleRate  = 3;
@@ -820,6 +824,13 @@ void Potentials::initialize(Ions& ions)
             initializeRadialDataOnMesh(position, sp);
         }
     }
+
+    // rescale rho_comp_ due to finite mesh effects
+    rescaleRhoComp();
+
+    initBackground();
+
+    addBackgroundToRhoComp();
 }
 
 void Potentials::rescaleRhoComp()
@@ -831,24 +842,33 @@ void Potentials::rescaleRhoComp()
     const pb::Grid& mygrid = mymesh->grid();
 
     // Check compensating charges
-    double comp_rho = getCharge(&rho_comp_[0]);
+    double comp_rho = getCharge(rho_comp_.data());
+    if (onpe0 && ct.verbose > 1)
+    {
+        std::cout << std::setprecision(8) << std::fixed
+                  << " Charge of rhoc: " << comp_rho << std::endl;
+    }
 
     if (onpe0 && ct.verbose > 1)
     {
-        cout << " Rescaling rhoc" << endl;
+        std::cout << " Rescaling rhoc" << std::endl;
     }
+
+    // rescale rho_comp_ (initialized by sampling on mesh)
+    // so that its integral exactly matches ionic_charge_
     if (ionic_charge_ > 0.)
     {
         const int numpt = mygrid.size();
         double t        = ionic_charge_ / comp_rho;
-        LinearAlgebraUtils<MemorySpace::Host>::MPscal(numpt, t, &rho_comp_[0]);
+        LinearAlgebraUtils<MemorySpace::Host>::MPscal(
+            numpt, t, rho_comp_.data());
 
         // Check new compensating charges
-        comp_rho = getCharge(&rho_comp_[0]);
+        comp_rho = getCharge(rho_comp_.data());
     }
     if (onpe0 && ct.verbose > 1)
-        cout << " Rescaled compensating charges: " << setprecision(8) << fixed
-             << comp_rho << endl;
+        std::cout << " Rescaled compensating charges: " << std::setprecision(8)
+                  << std::fixed << comp_rho << std::endl;
     if (comp_rho < 0.) mmpi.abort();
 }
 
@@ -868,26 +888,23 @@ void Potentials::addBackgroundToRhoComp()
         {
             if (onpe0)
             {
-                cout << setprecision(12) << scientific
-                     << "Add background charge " << background << " to rhoc "
-                     << endl;
+                std::cout << std::setprecision(12) << std::scientific
+                          << "Add background charge " << background
+                          << " to rhoc " << std::endl;
             }
             for (int i = 0; i < numpt; i++)
                 rho_comp_[i] += background;
 
             // Check new compensating charges
-            getCharge(&rho_comp_[0]);
+            getCharge(rho_comp_.data());
         }
     }
 }
 
-void Potentials::initBackground(Ions& ions)
+void Potentials::initBackground()
 {
     Control& ct = *(Control::instance());
 
-    // Count up the total ionic charge
-    ionic_charge_ = ions.computeIonicCharge();
-
     // calculation the compensating background charge
     //   for charged supercell calculations
     background_charge_ = 0.;
@@ -898,14 +915,76 @@ void Potentials::initBackground(Ions& ions)
     }
     if (onpe0 && ct.verbose > 0)
     {
-        cout << "N electrons=      " << ct.getNel() << endl;
-        cout << "ionic charge=     " << ionic_charge_ << endl;
-        cout << "background charge=" << background_charge_ << endl;
+        std::cout << "N electrons=      " << ct.getNel() << std::endl;
+        std::cout << "ionic charge=     " << ionic_charge_ << std::endl;
+        std::cout << "background charge=" << background_charge_ << std::endl;
     }
 
     if (fabs(background_charge_) < 1.e-10) background_charge_ = 0.;
 }
 
+int Potentials::read(HDFrestart& h5f_file)
+{
+    Control& ct = *(Control::instance());
+
+    // Read total potential
+    h5f_file.read_1func_hdf5(vtot_.data(), "Vtotal");
+
+    // Read the hartree potential
+    h5f_file.read_1func_hdf5(vh_rho_.data(), "Hartree");
+
+    // Read dielectric potential
+    if (ct.diel)
+    {
+        h5f_file.read_1func_hdf5(vepsilon_.data(), "VDielectric");
+    }
+
+    return 0;
+}
+
+int Potentials::write(HDFrestart& h5f_file)
+{
+    Control& ct = *(Control::instance());
+
+    Mesh* mymesh           = Mesh::instance();
+    const pb::Grid& mygrid = mymesh->grid();
+
+    double ll[3]     = { mygrid.ll(0), mygrid.ll(1), mygrid.ll(2) };
+    double origin[3] = { mygrid.origin(0), mygrid.origin(1), mygrid.origin(2) };
+
+    // Write total potential
+    int ierr
+        = h5f_file.write_1func_hdf5(vtot_.data(), "Vtotal", &ll[0], &origin[0]);
+    if (ierr < 0) return ierr;
+
+    // Write the hartree potential
+    ierr = h5f_file.write_1func_hdf5(
+        vh_rho_.data(), "Hartree", &ll[0], &origin[0]);
+    if (ierr < 0) return ierr;
+
+    if (ct.AtomsDynamic() == AtomsDynamicType::MD)
+    {
+        // Write hartree potential before extrapolation
+        ierr = h5f_file.write_1func_hdf5(
+            vh_rho_backup_.data(), "Preceding_Hartree", &ll[0], &origin[0]);
+        if (ierr < 0) return ierr;
+    }
+
+    // Write
+    if (ct.diel)
+    {
+        ierr = h5f_file.write_1func_hdf5(
+            vepsilon_.data(), "VDielectric", &ll[0], &origin[0]);
+    }
+    if (ierr < 0) return ierr;
+
+    // Write external potential
+    ierr = h5f_file.write_1func_hdf5(v_ext_.data(), "Vext", &ll[0], &origin[0]);
+    if (ierr < 0) return ierr;
+
+    return ierr;
+}
+
 template void Potentials::setVxc<double>(
     const double* const vxc, const int iterativeIndex);
 template void Potentials::setVxc<float>(
diff --git a/src/Potentials.h b/src/Potentials.h
index 71819a75..cb0f8e05 100644
--- a/src/Potentials.h
+++ b/src/Potentials.h
@@ -10,6 +10,7 @@
 #ifndef MGMOL_POTENTIALS_H
 #define MGMOL_POTENTIALS_H
 
+#include "HDFrestart.h"
 #include "Rho.h"
 #include "TriCubic.h"
 
@@ -100,6 +101,12 @@ class Potentials
     void initializeSupersampledRadialDataOnMesh(
         const Vector3D& position, const Species& sp);
 
+    void rescaleRhoComp();
+
+    void addBackgroundToRhoComp();
+
+    void initBackground();
+
 public:
     Potentials();
 
@@ -139,27 +146,22 @@ class Potentials
 
     void turnOnDiel() { diel_ = true; }
 
+    int write(HDFrestart& h5f_file);
+    int read(HDFrestart& h5f_file);
+
     int size() const { return size_; }
 
     double scf_dvrho(void) const { return scf_dvrho_; }
     double scf_dv(void) const { return scf_dv_; }
     POTDTYPE* vtot() { return vtot_.data(); }
-    POTDTYPE* vh_rho() { return vh_rho_.data(); }
     RHODTYPE* rho_comp() { return rho_comp_.data(); }
 
     const std::vector<POTDTYPE>& vnuc() const { return v_nuc_; }
-    POTDTYPE* vnuc() { return v_nuc_.data(); }
-    POTDTYPE* vext() { return v_ext_.data(); }
+    const std::vector<POTDTYPE>& vh_rho() const { return vh_rho_; }
+
     POTDTYPE* vepsilon() { return vepsilon_.data(); }
-    POTDTYPE* vh_rho_backup() { return vh_rho_backup_.data(); }
 
     void axpVcompToVh(const double alpha);
-    void axpVcomp(POTDTYPE* v, const double alpha);
-
-    POTDTYPE vtot(const int i) { return vtot_[i]; }
-    POTDTYPE vh_rho(const int i) { return vh_rho_[i]; }
-    POTDTYPE vxc_rho(const int i) { return vxc_rho_[i]; }
-    POTDTYPE vepsilon(const int i) { return vepsilon_[i]; }
 
     bool diel() const { return diel_; }
 
@@ -198,9 +200,6 @@ class Potentials
     void setVh(const pb::GridFunc<POTDTYPE>& vh, const int iterativeIndex);
 
     void initialize(Ions& ions);
-    void rescaleRhoComp();
-    void initBackground(Ions& ions);
-    void addBackgroundToRhoComp();
 
     /*!
      * Save current Hartree potential into backup array
diff --git a/src/restart.cc b/src/restart.cc
index 30bdad73..ed7f478c 100644
--- a/src/restart.cc
+++ b/src/restart.cc
@@ -43,27 +43,11 @@ int MGmol<OrbitalsType>::read_rho_and_pot_hdf5(
         os_ << "Try to read density and potentials" << std::endl;
 
     Potentials& pot = hamiltonian_->potential();
+    pot.read(file);
 
-    Mesh* mymesh           = Mesh::instance();
-    const pb::Grid& mygrid = mymesh->grid();
-    POTDTYPE* tmp          = new POTDTYPE[mygrid.size()];
-
-    // Read total potential
-    file.read_1func_hdf5(pot.vtot(), "Vtotal");
-
-    // Read the hartree potential
-    file.read_1func_hdf5(tmp, "Hartree");
-    pot.setVh(tmp, 0);
-
-    // Read dielectric potential
-    if (ct.diel)
-    {
-        file.read_1func_hdf5(pot.vepsilon(), "VDielectric");
-    }
     // Read the Density
     rho.readRestart(file);
 
-    delete[] tmp;
     return 0;
 }
 
@@ -131,41 +115,12 @@ int MGmol<OrbitalsType>::write_hdf5(HDFrestart& h5f_file,
 
     if (ct.out_restart_info > 1)
     {
-        // Write total potential
-        int ierr = h5f_file.write_1func_hdf5(
-            pot.vtot(), "Vtotal", &ll[0], &origin[0]);
-        if (ierr < 0) return ierr;
-
-        // Write the hartree potential
-        ierr = h5f_file.write_1func_hdf5(
-            pot.vh_rho(), "Hartree", &ll[0], &origin[0]);
-        if (ierr < 0) return ierr;
-
-        if (ct.AtomsDynamic() == AtomsDynamicType::MD)
-        {
-            // Write hartree potential before extrapolation
-            ierr = h5f_file.write_1func_hdf5(
-                pot.vh_rho_backup(), "Preceding_Hartree", &ll[0], &origin[0]);
-            if (ierr < 0) return ierr;
-        }
-
-        // Write
-        if (ct.diel)
-        {
-            ierr = h5f_file.write_1func_hdf5(
-                pot.vepsilon(), "VDielectric", &ll[0], &origin[0]);
-        }
-        if (ierr < 0) return ierr;
+        int ierr = pot.write(h5f_file);
 
         // Write the Density
         ierr = h5f_file.write_1func_hdf5(
             &rho[0][0], "Density", &ll[0], &origin[0]);
         if (ierr < 0) return ierr;
-
-        // Write external potential
-        ierr
-            = h5f_file.write_1func_hdf5(pot.vext(), "Vext", &ll[0], &origin[0]);
-        if (ierr < 0) return ierr;
     }
 
     // Write wavefunctions and old centers.