diff --git a/src/Potentials.cc b/src/Potentials.cc
index bed005db..abe277f1 100644
--- a/src/Potentials.cc
+++ b/src/Potentials.cc
@@ -75,10 +75,13 @@ Potentials::Potentials()
 
     dv_.resize(size_);
 
-    memset(&vepsilon_[0], 0, size_ * sizeof(POTDTYPE));
-    memset(&vh_rho_[0], 0, size_ * sizeof(POTDTYPE));
-    memset(&vxc_rho_[0], 0, size_ * sizeof(POTDTYPE));
-    memset(&v_ext_[0], 0, size_ * sizeof(POTDTYPE));
+    vh_rho_backup_.resize(size_);
+
+    memset(vepsilon_.data(), 0, size_ * sizeof(POTDTYPE));
+    memset(vh_rho_.data(), 0, size_ * sizeof(POTDTYPE));
+    memset(vxc_rho_.data(), 0, size_ * sizeof(POTDTYPE));
+    memset(v_ext_.data(), 0, size_ * sizeof(POTDTYPE));
+    memset(vh_rho_backup_.data(), 0, size_ * sizeof(POTDTYPE));
 
 #ifdef HAVE_TRICUBIC
     vext_tricubic_ = NULL;
@@ -596,6 +599,11 @@ 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));
+}
+
 void Potentials::initializeSupersampledRadialDataOnMesh(
     const Vector3D& position, const Species& sp)
 {
diff --git a/src/Potentials.h b/src/Potentials.h
index 970cf95d..71819a75 100644
--- a/src/Potentials.h
+++ b/src/Potentials.h
@@ -76,6 +76,11 @@ class Potentials
 
     std::vector<POTDTYPE> dv_;
 
+    /*!
+     * Backpup copy of Hartree potential to save previous state
+     */
+    std::vector<POTDTYPE> vh_rho_backup_;
+
     int itindex_vxc_;
     int itindex_vh_;
 
@@ -138,14 +143,15 @@ class Potentials
 
     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]; }
-    RHODTYPE* rho_comp() { return &rho_comp_[0]; }
+    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_[0]; }
-    POTDTYPE* vext() { return &v_ext_[0]; }
-    POTDTYPE* vepsilon() { return &vepsilon_[0]; }
+    POTDTYPE* vnuc() { return v_nuc_.data(); }
+    POTDTYPE* vext() { return v_ext_.data(); }
+    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);
@@ -196,6 +202,11 @@ class Potentials
     void initBackground(Ions& ions);
     void addBackgroundToRhoComp();
 
+    /*!
+     * Save current Hartree potential into backup array
+     */
+    void backupVh();
+
 #ifdef HAVE_TRICUBIC
     void readExternalPot(const string filename, const char type);
     void setupVextTricubic();
diff --git a/src/md.cc b/src/md.cc
index b50618c8..dfb4dfe6 100644
--- a/src/md.cc
+++ b/src/md.cc
@@ -54,6 +54,9 @@ void MGmol<OrbitalsType>::moveVnuc(Ions& ions)
 
     Potentials& pot = hamiltonian_->potential();
 
+    // save Hartree potential internally
+    pot.backupVh();
+
     // Update items that change when the ionic coordinates change
     pot.axpVcompToVh(1.);
     initNuc(ions);
diff --git a/src/restart.cc b/src/restart.cc
index 2fa76ecb..fde45b2c 100644
--- a/src/restart.cc
+++ b/src/restart.cc
@@ -141,6 +141,14 @@ int MGmol<OrbitalsType>::write_hdf5(HDFrestart& h5f_file,
             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)
         {