Save Hartree potential for write in restart file

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)
 {

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();

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);

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)
         {