diff --git a/src/DensityMatrix.cc b/src/DensityMatrix.cc
index e2c7b473..2f5e7c1c 100644
--- a/src/DensityMatrix.cc
+++ b/src/DensityMatrix.cc
@@ -12,6 +12,8 @@
 #include "DistMatrix.h"
 #include "MGmol_MPI.h"
 #include "ReplicatedMatrix.h"
+#include "ReplicatedWorkSpace.h"
+#include "hdf_tools.h"
 
 #include <cmath>
 #include <iomanip>
@@ -22,16 +24,25 @@ const double factor_kernel4dot = 10.;
 
 #define PROCRUSTES 0
 
+#define MGMOL_DENSITYMATRIX_FAIL(X)                                            \
+    {                                                                          \
+        std::cerr << "DensityMatrix failure:" << std::endl;                    \
+        std::cerr << "Error Message: " << X << std::endl;                      \
+    }
+
 // occupations in [0,1]
 // DM eigenvalues in [0,orbital_occupation]
 
 template <class MatrixType>
 DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
+    : dim_(ndim),
+      orbitals_index_(-1),
+      occ_uptodate_(false),
+      uniform_occ_(false),
+      stripped_(false)
 {
     assert(ndim > 0);
 
-    dim_ = ndim;
-
     occ_uptodate_ = false;
     stripped_     = false;
     uniform_occ_  = false;
@@ -39,8 +50,6 @@ DensityMatrix<MatrixType>::DensityMatrix(const int ndim)
     MGmol_MPI& mmpi     = *(MGmol_MPI::instance());
     orbital_occupation_ = mmpi.nspin() > 1 ? 1. : 2.;
 
-    orbitals_index_ = -1;
-
     dm_         = new MatrixType("DM", ndim, ndim);
     kernel4dot_ = new MatrixType("K4dot", ndim, ndim);
     work_       = new MatrixType("work", ndim, ndim);
@@ -438,6 +447,43 @@ void DensityMatrix<MatrixType>::mix(
     orbitals_index_ = new_orbitals_index;
 }
 
+template <class MatrixType>
+int DensityMatrix<MatrixType>::write(HDFrestart& h5f_file, std::string& name)
+{
+    ReplicatedWorkSpace<double>& wspace(
+        ReplicatedWorkSpace<double>::instance());
+
+    wspace.initSquareMatrix(*dm_);
+
+    DISTMATDTYPE* work_matrix = wspace.square_matrix();
+
+    hid_t file_id = h5f_file.file_id();
+    return mgmol_tools::write_matrix(file_id, name, work_matrix, dim_);
+}
+
+template <class MatrixType>
+int DensityMatrix<MatrixType>::read(HDFrestart& h5f_file, std::string& name)
+{
+    ReplicatedWorkSpace<double>& wspace(
+        ReplicatedWorkSpace<double>::instance());
+    DISTMATDTYPE* work_matrix = wspace.square_matrix();
+
+    MGmol_MPI& mmpi = *(MGmol_MPI::instance());
+
+    int ierr = 0;
+    if (mmpi.instancePE0())
+    {
+        hid_t file_id = h5f_file.file_id();
+        ierr          = mgmol_tools::read_matrix(file_id, name, work_matrix);
+    }
+    mmpi.bcast(&ierr, 1);
+
+    if (ierr >= 0) wspace.mpiBcastSquareMatrix();
+    if (ierr >= 0) initMatrix(work_matrix);
+
+    return ierr;
+}
+
 template class DensityMatrix<dist_matrix::DistMatrix<double>>;
 #ifdef HAVE_MAGMA
 template class DensityMatrix<ReplicatedMatrix>;
diff --git a/src/DensityMatrix.h b/src/DensityMatrix.h
index 84804a86..1d09e79e 100644
--- a/src/DensityMatrix.h
+++ b/src/DensityMatrix.h
@@ -10,6 +10,7 @@
 #ifndef MGMOL_DENSITYMATRIX_H
 #define MGMOL_DENSITYMATRIX_H
 
+#include "HDFrestart.h"
 #include "MGmol_MPI.h"
 #include "global.h"
 
@@ -23,7 +24,7 @@
 template <class MatrixType>
 class DensityMatrix
 {
-    int dim_;
+    const int dim_;
     std::vector<double> occupation_;
 
     MatrixType* dm_;
@@ -35,6 +36,10 @@ class DensityMatrix
     bool occ_uptodate_;
     bool uniform_occ_;
     bool stripped_;
+
+    /*!
+     * Max. occupation of an orbital: 1 with spin, 2 otherwise
+     */
     double orbital_occupation_;
 
     DensityMatrix();
@@ -146,6 +151,9 @@ class DensityMatrix
     double getExpectation(const MatrixType& A);
     void mix(
         const double mix, const MatrixType& matA, const int new_orbitals_index);
+
+    int write(HDFrestart& h5f_file, std::string& name);
+    int read(HDFrestart& h5f_file, std::string& name);
 };
 
 #endif
diff --git a/src/ExtendedGridOrbitals.cc b/src/ExtendedGridOrbitals.cc
index 12e5b04c..19c0d26b 100644
--- a/src/ExtendedGridOrbitals.cc
+++ b/src/ExtendedGridOrbitals.cc
@@ -584,7 +584,6 @@ int ExtendedGridOrbitals::read_hdf5(HDFrestart& h5f_file)
 
     Control& ct = *(Control::instance());
 
-    hid_t file_id    = h5f_file.file_id();
     std::string name = "Function";
     int ierr         = read_func_hdf5(h5f_file, name);
     if (ierr < 0)
@@ -603,7 +602,7 @@ int ExtendedGridOrbitals::read_hdf5(HDFrestart& h5f_file)
     // Read DM
     if (!ct.fullyOccupied())
     {
-        ierr = proj_matrices_->read_dm_hdf5(file_id);
+        ierr = proj_matrices_->readDM(h5f_file);
         if (ierr < 0)
         {
             (*MPIdata::serr)
@@ -618,28 +617,7 @@ int ExtendedGridOrbitals::read_hdf5(HDFrestart& h5f_file)
     return ierr;
 }
 
-int ExtendedGridOrbitals::write_hdf5(
-    HDFrestart& h5f_file, const std::string& name)
-{
-    assert(proj_matrices_ != nullptr);
-    Control& ct = *(Control::instance());
-
-    if (!ct.fullyOccupied())
-    {
-        MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
-        mmpi.barrier();
-
-        int ierr = proj_matrices_->writeDM_hdf5(h5f_file);
-        if (ierr < 0) return ierr;
-    }
-
-    int ierr = write_func_hdf5(h5f_file, name);
-
-    return ierr;
-}
-
-int ExtendedGridOrbitals::write_func_hdf5(
-    HDFrestart& h5f_file, const std::string& name)
+int ExtendedGridOrbitals::write(HDFrestart& h5f_file, const std::string& name)
 {
     if (onpe0)
         (*MPIdata::sout) << "ExtendedGridOrbitals::write_func_hdf5()...\n";
@@ -829,7 +807,7 @@ int ExtendedGridOrbitals::read_func_hdf5(
         const std::string datasetname(getDatasetName(name, icolor));
 
         // check if dataset exists...
-        int err_id = h5f_file.dset_exists(datasetname);
+        int err_id = h5f_file.checkDataExists(datasetname);
         if (h5f_file.gatherDataX()) mmpi.bcast(&err_id, 1);
         if (err_id == 0) break; // dataset does not exists
 
diff --git a/src/ExtendedGridOrbitals.h b/src/ExtendedGridOrbitals.h
index fb7a4b42..c1a718d0 100644
--- a/src/ExtendedGridOrbitals.h
+++ b/src/ExtendedGridOrbitals.h
@@ -374,8 +374,7 @@ class ExtendedGridOrbitals : public Orbitals
     void multiplyByMatrix2states(const int st1, const int st2,
         const double* mat, ExtendedGridOrbitals& product);
 
-    int write_hdf5(HDFrestart& h5f_file, const std::string& name = "Function");
-    int write_func_hdf5(HDFrestart&, const std::string& name = "Function");
+    int write(HDFrestart&, const std::string& name = "Function");
     int read_hdf5(HDFrestart& h5f_file);
     int read_func_hdf5(HDFrestart&, const std::string& name = "Function");
 
diff --git a/src/HDFrestart.cc b/src/HDFrestart.cc
index f95292ac..08018d47 100644
--- a/src/HDFrestart.cc
+++ b/src/HDFrestart.cc
@@ -933,12 +933,12 @@ int HDFrestart::getLRCenters(std::multimap<std::string, Vector3D>& centers,
 
         std::string datasetname(getDatasetName(name, color));
 
-        int err_id = dset_exists(datasetname);
+        int err_id = checkDataExistsLocal(datasetname);
         if (err_id == 0)
         { // dataset does not exists
             // try older version
             datasetname = getDatasetName_old(name, color);
-            err_id      = dset_exists(datasetname);
+            err_id      = checkDataExistsLocal(datasetname);
         }
 
         if (err_id == 0)
@@ -1052,12 +1052,12 @@ int HDFrestart::getLRs(std::shared_ptr<LocalizationRegions> lrs,
 
         std::string datasetname(getDatasetName(name, color));
 
-        int err_id = dset_exists(datasetname);
+        int err_id = checkDataExistsLocal(datasetname);
         if (err_id == 0)
         { // dataset does not exists
             // try older version
             datasetname = getDatasetName_old(name, color);
-            err_id      = dset_exists(datasetname);
+            err_id      = checkDataExistsLocal(datasetname);
         }
         if (err_id == 0)
         { // dataset does not exists
diff --git a/src/HDFrestart.h b/src/HDFrestart.h
index fef634c2..7ffadff7 100644
--- a/src/HDFrestart.h
+++ b/src/HDFrestart.h
@@ -155,28 +155,30 @@ class HDFrestart
         return 0;
     }
 
-    hid_t dset_exists(const std::string& datasetname) const
+    hid_t checkDataExists(const std::string& datasetname) const
     {
-        if (active_)
-        {
-            return dset_exists(datasetname.c_str());
-        }
-        return 0;
+        return checkDataExists(datasetname.c_str());
     }
-    herr_t dset_exists(const char* const datasetname) const
+    herr_t checkDataExists(const char* const datasetname) const
+    {
+        herr_t err_id = checkDataExistsLocal(datasetname);
+
+        short id = (short)err_id;
+        MPI_Bcast(&id, 1, MPI_SHORT, 0, comm_data_);
+        return (herr_t)id;
+    }
+
+    hid_t checkDataExistsLocal(const std::string& datasetname) const
+    {
+        return checkDataExistsLocal(datasetname.c_str());
+    }
+
+    herr_t checkDataExistsLocal(const char* const datasetname) const
     {
         herr_t err_id = 0;
         if (active_)
         {
             err_id = H5LTfind_dataset(file_id_, datasetname);
-
-            if (err_id < 0)
-            {
-                if (onpe0)
-                    (*MPIdata::sout)
-                        << "HDFrestart::dset_exists() failed for dataset "
-                        << datasetname << std::endl;
-            }
         }
         return err_id;
     }
diff --git a/src/LBFGS_IonicStepper.cc b/src/LBFGS_IonicStepper.cc
index caea5f84..93ce03f5 100644
--- a/src/LBFGS_IonicStepper.cc
+++ b/src/LBFGS_IonicStepper.cc
@@ -328,11 +328,11 @@ int LBFGS_IonicStepper::read_lbfgs(HDFrestart& h5f_file)
 
     short check_data = 0;
     std::vector<double> attr_d(16);
+    // Open an existing dataset.
+    std::string datasetname("/LBFGS");
+    int err_id = h5f_file.checkDataExists(datasetname);
     if (onpe0)
     {
-        // Open an existing dataset.
-        std::string datasetname("/LBFGS");
-        int err_id = h5f_file.dset_exists(datasetname);
         if (err_id < 0)
         { // dataset does not exists
             (*MPIdata::sout) << "Warning: no dataset " << datasetname
@@ -343,10 +343,9 @@ int LBFGS_IonicStepper::read_lbfgs(HDFrestart& h5f_file)
             dataset_id = H5Dopen2(file_id, "/LBFGS", H5P_DEFAULT);
             if (dataset_id < 0)
             {
-                if (onpe0)
-                    (*MPIdata::sout) << "Warning: H5Dopen failed for /LBFGS-> "
-                                        "no restart info for LBFGS"
-                                     << std::endl;
+                (*MPIdata::sout) << "Warning: H5Dopen failed for /LBFGS-> "
+                                    "no restart info for LBFGS"
+                                 << std::endl;
             }
             else
             {
diff --git a/src/LocGridOrbitals.cc b/src/LocGridOrbitals.cc
index ae8a4d1d..4a1b9442 100644
--- a/src/LocGridOrbitals.cc
+++ b/src/LocGridOrbitals.cc
@@ -966,7 +966,6 @@ int LocGridOrbitals::read_hdf5(HDFrestart& h5f_file)
 
     Control& ct = *(Control::instance());
 
-    hid_t file_id    = h5f_file.file_id();
     std::string name = "Function";
     int ierr         = read_func_hdf5(h5f_file, name);
     if (ierr < 0)
@@ -984,7 +983,7 @@ int LocGridOrbitals::read_hdf5(HDFrestart& h5f_file)
     // Read DM
     if (!ct.fullyOccupied())
     {
-        ierr = proj_matrices_->read_dm_hdf5(file_id);
+        ierr = proj_matrices_->readDM(h5f_file);
         if (ierr < 0)
         {
             (*MPIdata::serr)
@@ -997,27 +996,7 @@ int LocGridOrbitals::read_hdf5(HDFrestart& h5f_file)
     return ierr;
 }
 
-int LocGridOrbitals::write_hdf5(HDFrestart& h5f_file, const std::string& name)
-{
-    assert(proj_matrices_ != nullptr);
-    Control& ct = *(Control::instance());
-
-    if (!ct.fullyOccupied())
-    {
-        MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
-        mmpi.barrier();
-
-        int ierr = proj_matrices_->writeDM_hdf5(h5f_file);
-        if (ierr < 0) return ierr;
-    }
-
-    int ierr = write_func_hdf5(h5f_file, name);
-
-    return ierr;
-}
-
-int LocGridOrbitals::write_func_hdf5(
-    HDFrestart& h5f_file, const std::string& name)
+int LocGridOrbitals::write(HDFrestart& h5f_file, const std::string& name)
 {
     Control& ct   = *(Control::instance());
     hid_t file_id = h5f_file.file_id();
@@ -1249,7 +1228,7 @@ int LocGridOrbitals::read_func_hdf5(
         const std::string key(itcenter->first);
 
         // checkif dataset exists...
-        int err_id = h5f_file.dset_exists(key);
+        int err_id = h5f_file.checkDataExistsLocal(key);
         if (h5f_file.gatherDataX()) mmpi.bcast(&err_id, 1);
         if (err_id == 0) break; // dataset does not exists
 
diff --git a/src/LocGridOrbitals.h b/src/LocGridOrbitals.h
index 5616c6b6..951e418f 100644
--- a/src/LocGridOrbitals.h
+++ b/src/LocGridOrbitals.h
@@ -410,8 +410,7 @@ class LocGridOrbitals : public Orbitals
     void multiplyByMatrix2states(const int st1, const int st2,
         const double* mat, LocGridOrbitals& product);
 
-    int write_hdf5(HDFrestart& h5f_file, const std::string& name = "Function");
-    int write_func_hdf5(HDFrestart&, const std::string& name = "Function");
+    int write(HDFrestart&, const std::string& name = "Function");
     int read_hdf5(HDFrestart& h5f_file);
     int read_func_hdf5(HDFrestart&, const std::string& name = "Function");
 
diff --git a/src/MGmol.h b/src/MGmol.h
index f7525600..1eb8b2d3 100644
--- a/src/MGmol.h
+++ b/src/MGmol.h
@@ -133,7 +133,7 @@ class MGmol : public MGmolInterface
         KBPsiMatrixSparse* kbpsi, dist_matrix::DistMatrix<DISTMATDTYPE>& hij);
     void computeHnlPhiAndAdd2HPhi(Ions& ions, OrbitalsType& phi,
         OrbitalsType& hphi, const KBPsiMatrixSparse* const kbpsi);
-    int dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
+    int dumpMDrestartFile(OrbitalsType& orbitals, Ions& ions,
         Rho<OrbitalsType>& rho, const bool write_extrapolated_wf,
         const short count);
 
diff --git a/src/ProjectedMatrices.cc b/src/ProjectedMatrices.cc
index 3606aa40..74ff01ed 100644
--- a/src/ProjectedMatrices.cc
+++ b/src/ProjectedMatrices.cc
@@ -26,6 +26,7 @@
 #include "SparseDistMatrix.h"
 #include "SquareSubMatrix2DistMatrix.h"
 #include "fermi.h"
+#include "hdf_tools.h"
 
 #include <fstream>
 #include <iomanip>
@@ -732,114 +733,50 @@ double ProjectedMatrices<MatrixType>::checkCond(
             (*MPIdata::sout)
                 << " CONDITION NUMBER OF THE OVERLAP MATRIX EXCEEDS TOL: "
                 << rcond << "!!!" << std::endl;
-        Control& ct = *(Control::instance());
         if (flag) mmpi.abort();
     }
     return rcond;
 }
 
-////// TEMPLATE THIS FOR FLOAT OPTION ??
 template <class MatrixType>
-int ProjectedMatrices<MatrixType>::writeDM_hdf5(HDFrestart& h5f_file)
+int ProjectedMatrices<MatrixType>::writeDM(HDFrestart& h5f_file)
 {
-    hid_t file_id = h5f_file.file_id();
+    std::string name("/Density_Matrix");
+    return dm_->write(h5f_file, name);
+}
+
+template <class MatrixType>
+int ProjectedMatrices<MatrixType>::writeSavedDM(HDFrestart& h5f_file)
+{
+    std::string name("/Density_Matrix_WF");
 
     ReplicatedWorkSpace<double>& wspace(
         ReplicatedWorkSpace<double>::instance());
 
-    wspace.initSquareMatrix(dm_->getMatrix());
-
-    if (file_id < 0) return 0;
-
-    hsize_t dims[2] = { dim_, dim_ };
-
-    // filespace identifier
-    hid_t dataspace = H5Screate_simple(2, dims, nullptr);
-
-    hid_t dset_id = H5Dcreate2(file_id, "/Density_Matrix", H5T_NATIVE_DOUBLE,
-        dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
-    if (dset_id < 0)
-    {
-        (*MPIdata::serr) << "ProjectedMatrices<MatrixType>::write_dm_hdf5: "
-                            "H5Dcreate2 failed!!!"
-                         << std::endl;
-        return -1;
-    }
-
-    hid_t memspace  = dataspace;
-    hid_t filespace = dataspace;
+    const MatrixType* matrix = mat_X_old_.get();
+    wspace.initSquareMatrix(*matrix);
 
     DISTMATDTYPE* work_matrix = wspace.square_matrix();
-    herr_t status = H5Dwrite(dset_id, H5T_NATIVE_DOUBLE, memspace, filespace,
-        H5P_DEFAULT, work_matrix);
-    if (status < 0)
-    {
-        (*MPIdata::serr) << "Orbitals: H5Dwrite failed!!!" << std::endl;
-        return -1;
-    }
-
-    status = H5Dclose(dset_id);
-    if (status < 0)
-    {
-        (*MPIdata::serr) << "ProjectedMatrices<MatrixType>::write_dm_hdf5(), "
-                            "H5Dclose failed!!!"
-                         << std::endl;
-        return -1;
-    }
-    status = H5Sclose(dataspace);
-    if (status < 0)
-    {
-        (*MPIdata::serr) << "ProjectedMatrices<MatrixType>::write_dm_hdf5(), "
-                            "H5Sclose failed!!!"
-                         << std::endl;
-        return -1;
-    }
 
-    return 0;
+    hid_t file_id = h5f_file.file_id();
+    return mgmol_tools::write_matrix(file_id, name, work_matrix, dim_);
 }
-////// TEMPLATE THIS FOR FLOAT OPTION ??
+
 template <class MatrixType>
-int ProjectedMatrices<MatrixType>::read_dm_hdf5(hid_t file_id)
+int ProjectedMatrices<MatrixType>::readDM(HDFrestart& h5f_file)
 {
-    ReplicatedWorkSpace<double>& wspace(
-        ReplicatedWorkSpace<double>::instance());
-    DISTMATDTYPE* work_matrix = wspace.square_matrix();
+    std::string name("/Density_Matrix");
+    return dm_->read(h5f_file, name);
+}
 
-    int ierr        = 0;
+template <class MatrixType>
+int ProjectedMatrices<MatrixType>::readWFDM(HDFrestart& h5f_file)
+{
     MGmol_MPI& mmpi = *(MGmol_MPI::instance());
-    if (mmpi.instancePE0())
-    {
-        hid_t dset_id = H5Dopen2(file_id, "/Density_Matrix", H5P_DEFAULT);
-        if (dset_id < 0)
-        {
-            (*MPIdata::serr)
-                << "H5Dopen failed for /Density_Matrix!!!" << std::endl;
-        }
-        else
-        {
-            ierr          = 1;
-            herr_t status = H5Dread(dset_id, H5T_NATIVE_DOUBLE, H5S_ALL,
-                H5S_ALL, H5P_DEFAULT, work_matrix);
-            if (status < 0)
-            {
-                (*MPIdata::serr)
-                    << "H5Dread failed for /Density_Matrix!!!" << std::endl;
-                return -1;
-            }
-
-            status = H5Dclose(dset_id);
-            if (status < 0)
-            {
-                (*MPIdata::serr) << "H5Dclose failed!!!" << std::endl;
-                return -1;
-            }
-        }
-    }
-    mmpi.bcast(&ierr, 1);
-    if (ierr >= 0) wspace.mpiBcastSquareMatrix();
-    if (ierr >= 0) dm_->initMatrix(work_matrix);
-
-    return ierr;
+    mmpi.barrier();
+    if (mmpi.PE0()) std::cout << "ProjectedMatrices::readWFDM..." << std::endl;
+    std::string name("/Density_Matrix_WF");
+    return dm_->read(h5f_file, name);
 }
 
 template <class MatrixType>
diff --git a/src/ProjectedMatrices.h b/src/ProjectedMatrices.h
index 6725b827..f2881763 100644
--- a/src/ProjectedMatrices.h
+++ b/src/ProjectedMatrices.h
@@ -93,7 +93,14 @@ class ProjectedMatrices : public ProjectedMatricesInterface
     std::unique_ptr<MatrixType> matHB_;
     std::unique_ptr<MatrixType> matH_;
 
+    /*!
+     * Density Matrix
+     */
     std::unique_ptr<DensityMatrix<MatrixType>> dm_;
+
+    /*!
+     * Gram matrix of orbitals overlaps
+     */
     std::unique_ptr<GramMatrix<MatrixType>> gm_;
 
     // work matrix for tmp usage
@@ -311,8 +318,10 @@ class ProjectedMatrices : public ProjectedMatricesInterface
     double computeEntropy() override;
     double computeEntropyWithCheb(const double kbt);
     double checkCond(const double tol, const bool flag = true) override;
-    int writeDM_hdf5(HDFrestart& h5f_file) override;
-    int read_dm_hdf5(hid_t file_id) override;
+    int writeDM(HDFrestart& h5f_file) override;
+    int writeSavedDM(HDFrestart& h5f_file);
+    int readDM(HDFrestart& h5f_file) override;
+    int readWFDM(HDFrestart& h5f_file);
     void printEigenvalues(std::ostream& os) const;
     void updateDM(const int iterative_index) override;
     void updateDMwithEigenstates(const int iterative_index);
diff --git a/src/ProjectedMatricesInterface.h b/src/ProjectedMatricesInterface.h
index 653b386d..c52595e1 100644
--- a/src/ProjectedMatricesInterface.h
+++ b/src/ProjectedMatricesInterface.h
@@ -265,19 +265,35 @@ class ProjectedMatricesInterface : public ChebyshevApproximationFunction
 
         exitWithErrorMessage("updateDMwithRelax");
     }
-    virtual int read_dm_hdf5(hid_t file_id)
+    virtual int readDM(HDFrestart& h5f_file)
     {
-        (void)file_id;
+        (void)h5f_file;
+
+        exitWithErrorMessage("readDM");
+
+        return 0;
+    }
+    virtual int readWFDM(HDFrestart& h5f_file)
+    {
+        (void)h5f_file;
+
+        exitWithErrorMessage("readWFDM");
+
+        return 0;
+    }
+    virtual int writeDM(HDFrestart& h5f_file)
+    {
+        (void)h5f_file;
 
-        exitWithErrorMessage("read_dm_hdf5");
+        exitWithErrorMessage("writeDM");
 
         return 0;
     }
-    virtual int writeDM_hdf5(HDFrestart& h5f_file)
+    virtual int writeSavedDM(HDFrestart& h5f_file)
     {
         (void)h5f_file;
 
-        exitWithErrorMessage("writeDM_hdf5");
+        exitWithErrorMessage("writeSavedDM");
 
         return 0;
     }
diff --git a/src/hdf_tools.cc b/src/hdf_tools.cc
index 69a105dc..e56c0478 100644
--- a/src/hdf_tools.cc
+++ b/src/hdf_tools.cc
@@ -14,12 +14,16 @@
 #include <iostream>
 #include <string.h>
 
-using namespace std;
+#define MGMOL_HDF5_FAIL(X)                                                     \
+    {                                                                          \
+        std::cerr << "MGMOL_HDF5 failure:" << std::endl;                       \
+        std::cerr << "Error Message: " << X << std::endl;                      \
+    }
 
 namespace mgmol_tools
 {
-void write1d(
-    hid_t file_id, const string& datasetname, vector<int>& data, size_t length)
+void write1d(hid_t file_id, const std::string& datasetname,
+    std::vector<int>& data, size_t length)
 {
     assert(file_id >= 0);
 
@@ -31,7 +35,7 @@ void write1d(
     hid_t dataspace_id = H5Screate_simple(1, &dim, nullptr);
     if (dataspace_id < 0)
     {
-        cerr << "write1d(), H5Screate_simple failed!!!" << endl;
+        std::cerr << "write1d(), H5Screate_simple failed!!!" << std::endl;
         return;
     }
 
@@ -40,7 +44,7 @@ void write1d(
         dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
     if (dataset_id < 0)
     {
-        cerr << "write1d(), H5Dcreate2 failed!!!" << endl;
+        std::cerr << "write1d(), H5Dcreate2 failed!!!" << std::endl;
         return;
     }
     H5Sclose(dataspace_id);
@@ -49,24 +53,24 @@ void write1d(
         dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &data[0]);
     if (status < 0)
     {
-        cerr << "write1d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "write1d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
     status = H5Dclose(dataset_id);
     if (status < 0)
     {
-        cerr << "write1d(), H5Dclose failed!!!" << endl;
+        std::cerr << "write1d(), H5Dclose failed!!!" << std::endl;
         return;
     }
 }
 
-void write2d(
-    hid_t file_id, const string& datasetname, vector<int>& data, size_t* dims)
+void write2d(hid_t file_id, const std::string& datasetname,
+    std::vector<int>& data, size_t* dims)
 {
     assert(file_id >= 0);
 
-    // cout<<"Write "<<dims[0]<<" atomic numbers..."<<endl;
+    // std::cout<<"Write "<<dims[0]<<" atomic numbers..."<<endl;
 
     if (dims[0] == 0) return;
 
@@ -76,7 +80,7 @@ void write2d(
     hid_t dataspace_id = H5Screate_simple(2, dimsm, nullptr);
     if (dataspace_id < 0)
     {
-        cerr << "write2d(), H5Screate_simple failed!!!" << endl;
+        std::cerr << "write2d(), H5Screate_simple failed!!!" << std::endl;
         return;
     }
 
@@ -85,7 +89,7 @@ void write2d(
         dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
     if (dataset_id < 0)
     {
-        cerr << "write2d(), H5Dcreate2 failed!!!" << endl;
+        std::cerr << "write2d(), H5Dcreate2 failed!!!" << std::endl;
         return;
     }
     H5Sclose(dataspace_id);
@@ -94,24 +98,24 @@ void write2d(
         dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &data[0]);
     if (status < 0)
     {
-        cerr << "write2d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "write2d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
     status = H5Dclose(dataset_id);
     if (status < 0)
     {
-        cerr << "write2d(), H5Dclose failed!!!" << endl;
+        std::cerr << "write2d(), H5Dclose failed!!!" << std::endl;
         return;
     }
 }
 
-void write2d(hid_t file_id, const string& datasetname,
-    vector<unsigned short>& data, size_t* dims)
+void write2d(hid_t file_id, const std::string& datasetname,
+    std::vector<unsigned short>& data, size_t* dims)
 {
     assert(file_id >= 0);
 
-    // cout<<"Write "<<dims[0]<<" atomic numbers..."<<endl;
+    // std::cout<<"Write "<<dims[0]<<" atomic numbers..."<<endl;
 
     if (dims[0] == 0) return;
 
@@ -121,7 +125,7 @@ void write2d(hid_t file_id, const string& datasetname,
     hid_t dataspace_id = H5Screate_simple(2, dimsm, nullptr);
     if (dataspace_id < 0)
     {
-        cerr << "write2d(), H5Screate_simple failed!!!" << endl;
+        std::cerr << "write2d(), H5Screate_simple failed!!!" << std::endl;
         return;
     }
 
@@ -130,7 +134,7 @@ void write2d(hid_t file_id, const string& datasetname,
         H5T_NATIVE_USHORT, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
     if (dataset_id < 0)
     {
-        cerr << "write2d(), H5Dcreate2 failed!!!" << endl;
+        std::cerr << "write2d(), H5Dcreate2 failed!!!" << std::endl;
         return;
     }
     H5Sclose(dataspace_id);
@@ -139,24 +143,24 @@ void write2d(hid_t file_id, const string& datasetname,
         dataset_id, H5T_NATIVE_USHORT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &data[0]);
     if (status < 0)
     {
-        cerr << "write2d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "write2d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
     status = H5Dclose(dataset_id);
     if (status < 0)
     {
-        cerr << "write2d(), H5Dclose failed!!!" << endl;
+        std::cerr << "write2d(), H5Dclose failed!!!" << std::endl;
         return;
     }
 }
 
-void write2d(hid_t file_id, const string& datasetname, vector<double>& data,
-    size_t* dims)
+void write2d(hid_t file_id, const std::string& datasetname,
+    std::vector<double>& data, size_t* dims)
 {
     assert(file_id >= 0);
 
-    // cout<<"Write "<<dims[0]<<" atomic numbers..."<<endl;
+    // std::cout<<"Write "<<dims[0]<<" atomic numbers..."<<endl;
 
     if (dims[0] == 0) return;
 
@@ -166,7 +170,7 @@ void write2d(hid_t file_id, const string& datasetname, vector<double>& data,
     hid_t dataspace_id = H5Screate_simple(2, dimsm, nullptr);
     if (dataspace_id < 0)
     {
-        cerr << "write2d(), H5Screate_simple failed!!!" << endl;
+        std::cerr << "write2d(), H5Screate_simple failed!!!" << std::endl;
         return;
     }
 
@@ -175,7 +179,7 @@ void write2d(hid_t file_id, const string& datasetname, vector<double>& data,
         H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
     if (dataset_id < 0)
     {
-        cerr << "write2d(), H5Dcreate2 failed!!!" << endl;
+        std::cerr << "write2d(), H5Dcreate2 failed!!!" << std::endl;
         return;
     }
     H5Sclose(dataspace_id);
@@ -184,28 +188,28 @@ void write2d(hid_t file_id, const string& datasetname, vector<double>& data,
         dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &data[0]);
     if (status < 0)
     {
-        cerr << "write2d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "write2d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
     status = H5Dclose(dataset_id);
     if (status < 0)
     {
-        cerr << "write2d(), H5Dclose failed!!!" << endl;
+        std::cerr << "write2d(), H5Dclose failed!!!" << std::endl;
         return;
     }
 }
 
-void write2d(hid_t file_id, const string& datasetname, vector<string>& data,
-    size_t* dims)
+void write2d(hid_t file_id, const std::string& datasetname,
+    std::vector<std::string>& data, size_t* dims)
 {
     assert(file_id >= 0);
 
-    // create type for strings of length IonData_MaxStrLength
+    // create type for std::strings of length IonData_MaxStrLength
     hid_t strtype = H5Tcopy(H5T_C_S1);
     H5Tset_size(strtype, IonData_MaxStrLength);
 
-    // cout<<"Write "<<dims[0]<<" atomic numbers..."<<endl;
+    // std::cout<<"Write "<<dims[0]<<" atomic numbers..."<<endl;
 
     if (dims[0] == 0) return;
 
@@ -215,7 +219,7 @@ void write2d(hid_t file_id, const string& datasetname, vector<string>& data,
     hid_t dataspace_id = H5Screate_simple(2, dimsm, nullptr);
     if (dataspace_id < 0)
     {
-        cerr << "write2d(), H5Screate_simple failed!!!" << endl;
+        std::cerr << "write2d(), H5Screate_simple failed!!!" << std::endl;
         return;
     }
 
@@ -224,14 +228,15 @@ void write2d(hid_t file_id, const string& datasetname, vector<string>& data,
         dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
     if (dataset_id < 0)
     {
-        cerr << "write2d(), H5Dcreate2 failed!!!" << endl;
+        std::cerr << "write2d(), H5Dcreate2 failed!!!" << std::endl;
         return;
     }
     H5Sclose(dataspace_id);
 
     // First copy the contents of the vector into a temporary container
-    vector<FixedLengthString> tc;
-    for (vector<string>::const_iterator i = data.begin(), end = data.end();
+    std::vector<FixedLengthString> tc;
+    for (std::vector<std::string>::const_iterator i   = data.begin(),
+                                                  end = data.end();
          i != end; ++i)
     {
         FixedLengthString t;
@@ -239,7 +244,7 @@ void write2d(hid_t file_id, const string& datasetname, vector<string>& data,
         tc.push_back(t);
     }
 
-    string attname("String_Length");
+    std::string attname("String_Length");
     hsize_t dimsA[1]    = { 1 };
     hid_t dataspaceA_id = H5Screate_simple(1, dimsA, nullptr);
     hid_t attribute_id = H5Acreate2(dataset_id, attname.c_str(), H5T_NATIVE_INT,
@@ -248,15 +253,15 @@ void write2d(hid_t file_id, const string& datasetname, vector<string>& data,
         = H5Awrite(attribute_id, H5T_NATIVE_USHORT, &IonData_MaxStrLength);
     if (status < 0)
     {
-        cerr << "write2d(), Attribute: " << attname << " --- H5Awrite failed!!!"
-             << endl;
+        std::cerr << "write2d(), Attribute: " << attname
+                  << " --- H5Awrite failed!!!" << std::endl;
     }
 
     status
         = H5Dwrite(dataset_id, strtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &tc[0]);
     if (status < 0)
     {
-        cerr << "write2d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "write2d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
@@ -265,25 +270,25 @@ void write2d(hid_t file_id, const string& datasetname, vector<string>& data,
     status = H5Sclose(dataspaceA_id);
     if (status < 0)
     {
-        cerr << "write2d(), H5Sclose failed!!!" << endl;
+        std::cerr << "write2d(), H5Sclose failed!!!" << std::endl;
     }
 
     status = H5Aclose(attribute_id);
     if (status < 0)
     {
-        cerr << "write2d(), H5Aclose failed!!!" << endl;
+        std::cerr << "write2d(), H5Aclose failed!!!" << std::endl;
     }
 
     status = H5Dclose(dataset_id);
     if (status < 0)
     {
-        cerr << "write2d(), H5Dclose failed!!!" << endl;
+        std::cerr << "write2d(), H5Dclose failed!!!" << std::endl;
     }
 }
 
 #ifdef MGMOL_USE_HDF5P
-void parallelWrite2d(hid_t file_id, const string& datasetname,
-    vector<int>& data, size_t* dims, MPI_Comm comm)
+void parallelWrite2d(hid_t file_id, const std::string& datasetname,
+    std::vector<int>& data, size_t* dims, MPI_Comm comm)
 {
     assert(file_id >= 0);
     assert(!data.empty());
@@ -298,15 +303,17 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     hid_t filespace = H5Screate_simple(2, dimsf, nullptr);
     if (filespace < 0)
     {
-        cerr << "parallelWrite2d(), H5Screate_simple failed for filespace!!!"
-             << endl;
+        std::cerr
+            << "parallelWrite2d(), H5Screate_simple failed for filespace!!!"
+            << std::endl;
         return;
     }
     hid_t memspace = H5Screate_simple(2, dimsm, nullptr);
     if (memspace < 0)
     {
-        cerr << "parallelWrite2d(), H5Screate_simple failed for memspace!!!"
-             << endl;
+        std::cerr
+            << "parallelWrite2d(), H5Screate_simple failed for memspace!!!"
+            << std::endl;
         return;
     }
 
@@ -317,8 +324,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT);
     if (dset_id < 0)
     {
-        cerr << "parallelWrite2d() for dataset " << datasetname
-             << ", H5Dcreate2() failed!!!" << endl;
+        std::cerr << "parallelWrite2d() for dataset " << datasetname
+                  << ", H5Dcreate2() failed!!!" << std::endl;
         return;
     }
     H5Pclose(plist_id);
@@ -333,7 +340,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         filespace, H5S_SELECT_SET, offset, stride, count, block);
     if (status < 0)
     {
-        cerr << "parallelWrite2d(), H5Sselect_hyperslab() failed!!!" << endl;
+        std::cerr << "parallelWrite2d(), H5Sselect_hyperslab() failed!!!"
+                  << std::endl;
         return;
     }
 
@@ -345,7 +353,7 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         dset_id, H5T_NATIVE_INT, memspace, filespace, plist_id, &data[0]);
     if (status < 0)
     {
-        cerr << "parallelWrite2d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "parallelWrite2d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
@@ -355,8 +363,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     H5Sclose(memspace);
 }
 
-void parallelWrite2d(hid_t file_id, const string& datasetname,
-    vector<unsigned short>& data, size_t* dims, MPI_Comm comm)
+void parallelWrite2d(hid_t file_id, const std::string& datasetname,
+    std::vector<unsigned short>& data, size_t* dims, MPI_Comm comm)
 {
     assert(file_id >= 0);
     assert(!data.empty());
@@ -371,15 +379,17 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     hid_t filespace = H5Screate_simple(2, dimsf, nullptr);
     if (filespace < 0)
     {
-        cerr << "parallelWrite2d(), H5Screate_simple failed for filespace!!!"
-             << endl;
+        std::cerr
+            << "parallelWrite2d(), H5Screate_simple failed for filespace!!!"
+            << std::endl;
         return;
     }
     hid_t memspace = H5Screate_simple(2, dimsm, nullptr);
     if (memspace < 0)
     {
-        cerr << "parallelWrite2d(), H5Screate_simple failed for memspace!!!"
-             << endl;
+        std::cerr
+            << "parallelWrite2d(), H5Screate_simple failed for memspace!!!"
+            << std::endl;
         return;
     }
 
@@ -390,8 +400,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT);
     if (dset_id < 0)
     {
-        cerr << "parallelWrite2d() for dataset " << datasetname
-             << ", H5Dcreate2() failed!!!" << endl;
+        std::cerr << "parallelWrite2d() for dataset " << datasetname
+                  << ", H5Dcreate2() failed!!!" << std::endl;
         return;
     }
     H5Pclose(plist_id);
@@ -406,7 +416,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         filespace, H5S_SELECT_SET, offset, stride, count, block);
     if (status < 0)
     {
-        cerr << "parallelWrite2d(), H5Sselect_hyperslab() failed!!!" << endl;
+        std::cerr << "parallelWrite2d(), H5Sselect_hyperslab() failed!!!"
+                  << std::endl;
         return;
     }
 
@@ -418,7 +429,7 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         dset_id, H5T_NATIVE_USHORT, memspace, filespace, plist_id, &data[0]);
     if (status < 0)
     {
-        cerr << "parallelWrite2d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "parallelWrite2d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
@@ -428,8 +439,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     H5Sclose(memspace);
 }
 
-void parallelWrite2d(hid_t file_id, const string& datasetname,
-    vector<double>& data, size_t* dims, MPI_Comm comm)
+void parallelWrite2d(hid_t file_id, const std::string& datasetname,
+    std::vector<double>& data, size_t* dims, MPI_Comm comm)
 {
     assert(file_id >= 0);
     assert(!data.empty());
@@ -444,15 +455,17 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     hid_t filespace = H5Screate_simple(2, dimsf, nullptr);
     if (filespace < 0)
     {
-        cerr << "parallelWrite2d(), H5Screate_simple failed for filespace!!!"
-             << endl;
+        std::cerr
+            << "parallelWrite2d(), H5Screate_simple failed for filespace!!!"
+            << std::endl;
         return;
     }
     hid_t memspace = H5Screate_simple(2, dimsm, nullptr);
     if (memspace < 0)
     {
-        cerr << "parallelWrite2d(), H5Screate_simple failed for memspace!!!"
-             << endl;
+        std::cerr
+            << "parallelWrite2d(), H5Screate_simple failed for memspace!!!"
+            << std::endl;
         return;
     }
 
@@ -463,8 +476,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         filespace, H5P_DEFAULT, plist_id, H5P_DEFAULT);
     if (dset_id < 0)
     {
-        cerr << "parallelWrite2d() for dataset " << datasetname
-             << ", H5Dcreate2() failed!!!" << endl;
+        std::cerr << "parallelWrite2d() for dataset " << datasetname
+                  << ", H5Dcreate2() failed!!!" << std::endl;
         return;
     }
     H5Pclose(plist_id);
@@ -479,7 +492,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         filespace, H5S_SELECT_SET, offset, stride, count, block);
     if (status < 0)
     {
-        cerr << "parallelWrite2d(), H5Sselect_hyperslab() failed!!!" << endl;
+        std::cerr << "parallelWrite2d(), H5Sselect_hyperslab() failed!!!"
+                  << std::endl;
         return;
     }
 
@@ -491,7 +505,7 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, plist_id, &data[0]);
     if (status < 0)
     {
-        cerr << "parallelWrite2d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "parallelWrite2d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
@@ -501,13 +515,13 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     H5Sclose(memspace);
 }
 
-void parallelWrite2d(hid_t file_id, const string& datasetname,
-    vector<string>& data, size_t* dims, MPI_Comm comm)
+void parallelWrite2d(hid_t file_id, const std::string& datasetname,
+    std::vector<std::string>& data, size_t* dims, MPI_Comm comm)
 {
     assert(file_id >= 0);
     assert(!data.empty());
 
-    // create type for strings of length IonData_MaxStrLength
+    // create type for std::strings of length IonData_MaxStrLength
     hid_t strtype = H5Tcopy(H5T_C_S1);
     H5Tset_size(strtype, IonData_MaxStrLength);
 
@@ -521,15 +535,17 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     hid_t filespace = H5Screate_simple(2, dimsf, nullptr);
     if (filespace < 0)
     {
-        cerr << "parallelWrite2d(), H5Screate_simple failed for filespace!!!"
-             << endl;
+        std::cerr
+            << "parallelWrite2d(), H5Screate_simple failed for filespace!!!"
+            << std::endl;
         return;
     }
     hid_t memspace = H5Screate_simple(2, dimsm, nullptr);
     if (memspace < 0)
     {
-        cerr << "parallelWrite2d(), H5Screate_simple failed for memspace!!!"
-             << endl;
+        std::cerr
+            << "parallelWrite2d(), H5Screate_simple failed for memspace!!!"
+            << std::endl;
         return;
     }
 
@@ -540,8 +556,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         H5P_DEFAULT, plist_id, H5P_DEFAULT);
     if (dset_id < 0)
     {
-        cerr << "parallelWrite2d() for dataset " << datasetname
-             << ", H5Dcreate2() failed!!!" << endl;
+        std::cerr << "parallelWrite2d() for dataset " << datasetname
+                  << ", H5Dcreate2() failed!!!" << std::endl;
         return;
     }
     H5Pclose(plist_id);
@@ -556,7 +572,8 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
         filespace, H5S_SELECT_SET, offset, stride, count, block);
     if (status < 0)
     {
-        cerr << "parallelWrite2d(), H5Sselect_hyperslab() failed!!!" << endl;
+        std::cerr << "parallelWrite2d(), H5Sselect_hyperslab() failed!!!"
+                  << std::endl;
         return;
     }
 
@@ -565,8 +582,9 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
 
     // First copy the contents of the vector into a temporary container
-    vector<FixedLengthString> tc;
-    for (vector<string>::const_iterator i = data.begin(), end = data.end();
+    std::vector<FixedLengthString> tc;
+    for (std::vector<std::string>::const_iterator i   = data.begin(),
+                                                  end = data.end();
          i != end; ++i)
     {
         FixedLengthString t;
@@ -576,7 +594,7 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
     status = H5Dwrite(dset_id, strtype, memspace, filespace, plist_id, &tc[0]);
     if (status < 0)
     {
-        cerr << "parallelWrite2d(), H5Dwrite failed!!!" << endl;
+        std::cerr << "parallelWrite2d(), H5Dwrite failed!!!" << std::endl;
         return;
     }
 
@@ -589,7 +607,7 @@ void parallelWrite2d(hid_t file_id, const string& datasetname,
 #endif
 
 void addAttribute2Dataset(
-    hid_t dset_id, const char* attname, const vector<double>& attr_data)
+    hid_t dset_id, const char* attname, const std::vector<double>& attr_data)
 {
     assert(dset_id > -1);
 
@@ -600,30 +618,30 @@ void addAttribute2Dataset(
     hid_t dataspace_id = H5Screate_simple(1, &dim, nullptr);
     if (dataspace_id < 0)
     {
-        cerr << "H5Screate failed!!!" << endl;
+        std::cerr << "H5Screate failed!!!" << std::endl;
         return;
     }
     hid_t attribute_id = H5Acreate2(dset_id, attname, H5T_NATIVE_DOUBLE,
         dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
     if (attribute_id < 0)
     {
-        cerr << "H5Acreate failed!!!" << endl;
+        std::cerr << "H5Acreate failed!!!" << std::endl;
         return;
     }
 
     herr_t status = H5Sclose(dataspace_id);
-    if (status < 0) cerr << "H5Sclose failed!!!" << endl;
+    if (status < 0) std::cerr << "H5Sclose failed!!!" << std::endl;
 
     //(*MPIdata::sout)<<"Write attribute "<<attname<<endl;
     status = H5Awrite(attribute_id, H5T_NATIVE_DOUBLE, &attr_data[0]);
-    if (status < 0) cerr << "H5Awrite failed!!!" << endl;
+    if (status < 0) std::cerr << "H5Awrite failed!!!" << std::endl;
 
     status = H5Aclose(attribute_id);
-    if (status < 0) cerr << "H5Aclose failed!!!" << endl;
+    if (status < 0) std::cerr << "H5Aclose failed!!!" << std::endl;
 }
 
 void addAttribute2Dataset(
-    hid_t dset_id, const char* attname, const vector<int>& attr_data)
+    hid_t dset_id, const char* attname, const std::vector<int>& attr_data)
 {
     assert(dset_id > -1);
 
@@ -634,21 +652,21 @@ void addAttribute2Dataset(
     hid_t dataspace_id = H5Screate_simple(1, &dim, nullptr);
     if (dataspace_id < 0)
     {
-        cerr << "H5Screate failed!!!" << endl;
+        std::cerr << "H5Screate failed!!!" << std::endl;
         return;
     }
     hid_t attribute_id = H5Acreate2(dset_id, attname, H5T_NATIVE_INT,
         dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
     if (attribute_id < 0)
     {
-        cerr << "H5Acreate failed!!!" << endl;
+        std::cerr << "H5Acreate failed!!!" << std::endl;
         return;
     }
 
     herr_t status = H5Sclose(dataspace_id);
     if (status < 0)
     {
-        cerr << "H5Sclose failed!!!" << endl;
+        std::cerr << "H5Sclose failed!!!" << std::endl;
         return;
     }
 
@@ -656,14 +674,14 @@ void addAttribute2Dataset(
     status = H5Awrite(attribute_id, H5T_NATIVE_INT, &attr_data[0]);
     if (status < 0)
     {
-        cerr << "H5Awrite failed!!!" << endl;
+        std::cerr << "H5Awrite failed!!!" << std::endl;
         return;
     }
 
     status = H5Aclose(attribute_id);
     if (status < 0)
     {
-        cerr << "H5Aclose failed!!!" << endl;
+        std::cerr << "H5Aclose failed!!!" << std::endl;
     }
 }
 
@@ -676,7 +694,7 @@ int whatisopen(hid_t fid)
 
     if (cnt <= 0) return cnt;
 
-    if (cnt > 1) cout << "HDF5 file: " << cnt << " object(s) open\n";
+    if (cnt > 1) std::cout << "HDF5 file: " << cnt << " object(s) open\n";
 
     // objs = malloc(cnt * sizeof(hid_t));
     hid_t* objs = new hid_t[cnt];
@@ -699,4 +717,80 @@ int whatisopen(hid_t fid)
 
     return howmany;
 }
+
+int write_matrix(
+    hid_t file_id, std::string& name, const double* matrix, const int dim)
+{
+    if (file_id < 0) return 0;
+
+    hsize_t dims[2] = { (hsize_t)dim, (hsize_t)dim };
+
+    // filespace identifier
+    hid_t dataspace = H5Screate_simple(2, dims, nullptr);
+
+    hid_t dset_id = H5Dcreate2(file_id, name.c_str(), H5T_NATIVE_DOUBLE,
+        dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    if (dset_id < 0)
+    {
+        MGMOL_HDF5_FAIL("H5Dcreate2 failed!!!");
+        return -1;
+    }
+
+    hid_t memspace  = dataspace;
+    hid_t filespace = dataspace;
+
+    herr_t status = H5Dwrite(
+        dset_id, H5T_NATIVE_DOUBLE, memspace, filespace, H5P_DEFAULT, matrix);
+    if (status < 0)
+    {
+        MGMOL_HDF5_FAIL("H5Dwrite failed!!!");
+        return -1;
+    }
+
+    status = H5Dclose(dset_id);
+    if (status < 0)
+    {
+        MGMOL_HDF5_FAIL("H5Dclose failed!!!");
+        return -1;
+    }
+    status = H5Sclose(dataspace);
+    if (status < 0)
+    {
+        MGMOL_HDF5_FAIL("H5Sclose failed!!!");
+        return -1;
+    }
+
+    return 0;
+}
+
+int read_matrix(hid_t file_id, std::string& name, double* matrix)
+{
+    int ierr      = 0;
+    hid_t dset_id = H5Dopen2(file_id, name.c_str(), H5P_DEFAULT);
+    if (dset_id < 0)
+    {
+        MGMOL_HDF5_FAIL("H5Dopen failed!!");
+        ierr = -1;
+    }
+    else
+    {
+        herr_t status = H5Dread(
+            dset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, matrix);
+        if (status < 0)
+        {
+            MGMOL_HDF5_FAIL("H5Dread failed!!");
+            ierr = -1;
+        }
+
+        status = H5Dclose(dset_id);
+        if (status < 0)
+        {
+            MGMOL_HDF5_FAIL("H5Dclose failed!!!");
+            ierr = -1;
+        }
+    }
+
+    return ierr;
 }
+
+} // namespace
diff --git a/src/hdf_tools.h b/src/hdf_tools.h
index e1a0ed1f..7a64bdce 100644
--- a/src/hdf_tools.h
+++ b/src/hdf_tools.h
@@ -6,8 +6,8 @@
 // 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
-
 #include "hdf5.h"
+
 #include <mpi.h>
 #include <string>
 #include <vector>
@@ -39,4 +39,8 @@ void addAttribute2Dataset(
 void addAttribute2Dataset(
     hid_t dset_id, const char* attname, const std::vector<int>& attr_data);
 int whatisopen(hid_t fid);
+
+int write_matrix(
+    hid_t file_id, std::string& name, const double* matrix, const int dim);
+int read_matrix(hid_t file_id, std::string& name, double* matrix);
 }
diff --git a/src/md.cc b/src/md.cc
index 314fd0dc..36236f85 100644
--- a/src/md.cc
+++ b/src/md.cc
@@ -226,7 +226,7 @@ void checkMaxForces(const std::vector<double>& fion,
 }
 
 template <class OrbitalsType>
-int MGmol<OrbitalsType>::dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
+int MGmol<OrbitalsType>::dumpMDrestartFile(OrbitalsType& orbitals, Ions& ions,
     Rho<OrbitalsType>& rho, const bool write_extrapolated_wf, const short count)
 {
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
@@ -245,9 +245,11 @@ int MGmol<OrbitalsType>::dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
 
     HDFrestart h5file(filename, myPEenv, gdim, ct.out_restart_file_type);
 
-    OrbitalsType previous_orbitals("ForDumping", **orbitals, false);
+    OrbitalsType previous_orbitals("ForDumping", orbitals, false);
     if (!orbitals_extrapol_->getRestartData(previous_orbitals))
-        previous_orbitals.assign(**orbitals);
+        previous_orbitals.assign(orbitals);
+
+    // write all restart info in HDF5 file
     int ierr = write_hdf5(h5file, rho.rho_, ions, previous_orbitals, lrs_);
     mmpi.allreduce(&ierr, 1, MPI_MIN);
 
@@ -259,12 +261,11 @@ int MGmol<OrbitalsType>::dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
                 << std::endl;
         return ierr;
     }
-    // write_hdf5(h5file, rho.rho_, ions, *orbitals_minus1);
-    // stepper->write_hdf5(h5file);
 
     if (write_extrapolated_wf && ct.out_restart_info > 2)
     {
-        ierr = (*orbitals)->write_func_hdf5(h5file, "ExtrapolatedFunction");
+        // write extra  info needed for seamless MD restart
+        ierr = orbitals.write(h5file, "ExtrapolatedFunction");
         mmpi.allreduce(&ierr, 1, MPI_MIN);
         if (ierr < 0)
         {
@@ -274,6 +275,10 @@ int MGmol<OrbitalsType>::dumpMDrestartFile(OrbitalsType** orbitals, Ions& ions,
                                  << std::endl;
             return ierr;
         }
+
+        // write DM associated with non-extrapolated wavefunctions
+        // (last computed solution of KS equations)
+        proj_matrices_->writeSavedDM(h5file);
     }
 
     ierr = h5file.close();
@@ -354,16 +359,11 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
 
     if (ct.restart_info > 1)
     {
-        int flag_extrapolated_data = 0;
-        if (onpe0)
-        {
+        int flag_extrapolated_data
+            = h5f_file_->checkDataExists("ExtrapolatedFunction0000");
+        if (flag_extrapolated_data == 0)
             flag_extrapolated_data
-                = h5f_file_->dset_exists("ExtrapolatedFunction0000");
-            if (flag_extrapolated_data == 0)
-                flag_extrapolated_data
-                    = h5f_file_->dset_exists("ExtrapolatedFunction0");
-        }
-        MPI_Bcast(&flag_extrapolated_data, 1, MPI_INT, 0, comm_);
+                = h5f_file_->checkDataExists("ExtrapolatedFunction0");
 
         if (ct.restart_info > 2)
         {
@@ -377,7 +377,6 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
 
                 // need to reset a few things as we just read new orbitals
                 (*orbitals)->computeGramAndInvS();
-                dm_strategy_->update(*current_orbitals_);
             }
 
             DFTsolver<OrbitalsType>::setItCountLarge();
@@ -622,6 +621,13 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
             lrs_->clearOldCenters();
         }
 
+        // save DM for possible restart write
+        // note: extrapolation is going to modify it!
+        if ((ct.out_restart_info > 2)
+            && (((md_iteration_ % ct.checkpoint) == 0)
+                   || (mdstep == ct.num_MD_steps)))
+            proj_matrices_->saveDM();
+
         preWFextrapolation();
 
         if (ct.dt > 0.
@@ -653,7 +659,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
                     {
                         dump_tm_.start();
                         ierr = dumpMDrestartFile(
-                            orbitals, ions, *rho_, extrapolated_flag, count);
+                            **orbitals, ions, *rho_, extrapolated_flag, count);
                         dump_tm_.stop();
                         if (onpe0 && ierr < 0 && count < (DUMP_MAX_NUM_TRY - 1))
                             std::cout
@@ -692,7 +698,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
         {
             dump_tm_.start();
             ierr = dumpMDrestartFile(
-                orbitals, ions, *rho_, extrapolated_flag, count);
+                **orbitals, ions, *rho_, extrapolated_flag, count);
             dump_tm_.stop();
 
             if (onpe0 && ierr < 0 && count < (DUMP_MAX_NUM_TRY - 1))
@@ -712,6 +718,7 @@ void MGmol<OrbitalsType>::md(OrbitalsType** orbitals, Ions& ions)
 template <class OrbitalsType>
 void MGmol<OrbitalsType>::loadRestartFile(const std::string filename)
 {
+    if (onpe0) std::cout << "loadRestartFile..." << std::endl;
     MGmol_MPI& mmpi(*(MGmol_MPI::instance()));
     Control& ct              = *(Control::instance());
     Mesh* mymesh             = Mesh::instance();
@@ -734,6 +741,12 @@ void MGmol<OrbitalsType>::loadRestartFile(const std::string filename)
 
         global_exit(0);
     }
+    if (!ct.fullyOccupied())
+    {
+        // overwrite DM with restart data in dataset Density_Matrix_WF
+        if (h5file.checkDataExists("Density_Matrix_WF"))
+            ierr = proj_matrices_->readWFDM(h5file);
+    }
 
     ierr = h5file.close();
     mmpi.allreduce(&ierr, 1, MPI_MIN);
diff --git a/src/restart.cc b/src/restart.cc
index 944b1d8e..2fa76ecb 100644
--- a/src/restart.cc
+++ b/src/restart.cc
@@ -163,9 +163,14 @@ int MGmol<OrbitalsType>::write_hdf5(HDFrestart& h5f_file,
     // Write wavefunctions and old centers.
     if (ct.out_restart_info > 2)
     {
-        int ierr = orbitals.write_hdf5(h5f_file);
+        int ierr = orbitals.write(h5f_file);
         if (ierr < 0) return ierr;
 
+        if (!ct.fullyOccupied())
+        {
+            ierr = proj_matrices_->writeDM(h5f_file);
+            if (ierr < 0) return ierr;
+        }
         if (ct.isLocMode()
             && ct.WFExtrapolation() == WFExtrapolationType::Reversible)
         {
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 6de34971..53793986 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -247,6 +247,8 @@ add_executable(testDensityMatrix
                ${CMAKE_SOURCE_DIR}/src/tools/MGmol_MPI.cc
                ${CMAKE_SOURCE_DIR}/src/tools/mgmol_mpi_tools.cc
                ${CMAKE_SOURCE_DIR}/src/tools/random.cc
+               ${CMAKE_SOURCE_DIR}/src/ReplicatedWorkSpace.cc
+               ${CMAKE_SOURCE_DIR}/src/hdf_tools.cc
                ${CMAKE_SOURCE_DIR}/tests/ut_magma_main.cc)
 add_executable(testEnergyAndForces
                ${CMAKE_SOURCE_DIR}/tests/EnergyAndForces/testEnergyAndForces.cc)
@@ -491,6 +493,14 @@ add_test(NAME testMD_D72
          ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/coords.in
          ${CMAKE_CURRENT_SOURCE_DIR}/MD_D72/lrs.in
          ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
+add_test(NAME testMD_MVP
+         COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/MD_MVP/test.py
+         ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
+         ${CMAKE_CURRENT_BINARY_DIR}/../src/mgmol-opt
+         ${CMAKE_CURRENT_SOURCE_DIR}/MD_MVP/quench.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/MD_MVP/md.cfg
+         ${CMAKE_CURRENT_SOURCE_DIR}/MD_MVP/li2.xyz
+         ${CMAKE_CURRENT_SOURCE_DIR}/../potentials)
 add_test(NAME testLBFGS
          COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/LBFGS/test.py
          ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 4 ${MPIEXEC_PREFLAGS}
@@ -544,7 +554,7 @@ target_include_directories(testConditionDistMatrix PRIVATE ${Boost_INCLUDE_DIRS}
 target_include_directories(testConditionDistMatrixPower PRIVATE ${Boost_INCLUDE_DIRS})
 target_include_directories(testPower PRIVATE ${Boost_INCLUDE_DIRS})
 target_include_directories(testPowerDistMatrix PRIVATE ${Boost_INCLUDE_DIRS})
-target_include_directories(testDensityMatrix PRIVATE ${Boost_INCLUDE_DIRS})
+target_include_directories(testDensityMatrix PRIVATE ${Boost_INCLUDE_DIRS} ${HDF5_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})
@@ -558,6 +568,7 @@ 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)
+target_link_libraries(testDensityMatrix PRIVATE ${HDF5_LIBRARIES})
 
 if(${MAGMA_FOUND})
   target_link_libraries(testDistVector PRIVATE ${SCALAPACK_LIBRARIES}
diff --git a/tests/MD_MVP/li2.xyz b/tests/MD_MVP/li2.xyz
new file mode 100644
index 00000000..fce3cd72
--- /dev/null
+++ b/tests/MD_MVP/li2.xyz
@@ -0,0 +1,4 @@
+2
+
+Li  1.33 0.0 0.0
+Li -1.33 0.0 0.0
diff --git a/tests/MD_MVP/md.cfg b/tests/MD_MVP/md.cfg
new file mode 100644
index 00000000..83d5aa4c
--- /dev/null
+++ b/tests/MD_MVP/md.cfg
@@ -0,0 +1,39 @@
+verbosity=0
+xcFunctional=LDA
+FDtype=4th
+[Mesh]
+nx=32
+ny=32
+nz=32
+[Domain]
+ox=-5.
+oy=-5.
+oz=-5.
+lx=10.
+ly=10.
+lz=10.
+[Potentials]
+pseudopotential=pseudo.Li_ONCVPSP_LDA
+[Run]
+type=MD
+[MD]
+num_steps=5
+dt=15.
+[Quench]
+solver=PSD
+max_steps=30
+atol=1.e-8
+[Orbitals]
+nempty=1
+initial_type=random
+temperature=300.
+[ProjectedMatrices]
+solver=exact
+[DensityMatrix]
+solver=MVP
+nb_inner_it=1
+mixing=1.
+[Restart]
+input_filename=wave.out
+input_level=3
+output_level=3
diff --git a/tests/MD_MVP/quench.cfg b/tests/MD_MVP/quench.cfg
new file mode 100644
index 00000000..135dbde7
--- /dev/null
+++ b/tests/MD_MVP/quench.cfg
@@ -0,0 +1,34 @@
+verbosity=0
+xcFunctional=LDA
+FDtype=4th
+[Mesh]
+nx=32
+ny=32
+nz=32
+[Domain]
+ox=-5.
+oy=-5.
+oz=-5.
+lx=10.
+ly=10.
+lz=10.
+[Potentials]
+pseudopotential=pseudo.Li_ONCVPSP_LDA
+[Run]
+type=QUENCH
+[Quench]
+solver=PSD
+max_steps=120
+atol=1.e-8
+[Orbitals]
+nempty=1
+initial_type=random
+temperature=300.
+[ProjectedMatrices]
+solver=exact
+[DensityMatrix]
+solver=MVP
+nb_inner_it=1
+mixing=1.
+[Restart]
+output_level=3
diff --git a/tests/MD_MVP/test.py b/tests/MD_MVP/test.py
new file mode 100755
index 00000000..489c9c8c
--- /dev/null
+++ b/tests/MD_MVP/test.py
@@ -0,0 +1,89 @@
+#!/usr/bin/env python
+import sys
+import os
+import subprocess
+import string
+import shutil
+
+print("Test MD with MVP solver...")
+
+nargs=len(sys.argv)
+
+mpicmd = sys.argv[1]+" "+sys.argv[2]+" "+sys.argv[3]
+for i in range(4,nargs-5):
+  mpicmd = mpicmd + " "+sys.argv[i]
+print("MPI run command: {}".format(mpicmd))
+
+exe = sys.argv[nargs-5]
+inp1 = sys.argv[nargs-4]
+inp2 = sys.argv[nargs-3]
+coords = sys.argv[nargs-2]
+print("coordinates file: %s"%coords)
+
+#create links to potentials files
+dst = 'pseudo.Li_ONCVPSP_LDA'
+src = sys.argv[-1] + '/' + dst
+
+if not os.path.exists(dst):
+  print("Create link to %s"%dst)
+  os.symlink(src, dst)
+
+#run quench
+command = "{} {} -c {} -i {}".format(mpicmd,exe,inp1,coords)
+print("Run command: {}".format(command))
+output1 = subprocess.check_output(command,shell=True)
+lines=output1.split(b'\n')
+
+#analyse output of quench
+for line in lines:
+  num_matches = line.count(b'%%')
+  if num_matches:
+    print(line)
+
+#run MD
+for i in range(2):
+  command = "ls -ld snapshot0* | awk '{ print $9 }' | tail -n1"
+  print(command)
+  restart_file = subprocess.check_output(command,shell=True)
+  restart_file=str(restart_file[:-1],'utf-8')
+  print(restart_file)
+
+  os.rename(restart_file, 'wave.out')
+
+  #run MGmol
+  command = "{} {} -c {} -i {}".format(mpicmd,exe,inp2,coords)
+  output2 = subprocess.check_output(command,shell=True)
+
+  #remove used restart files
+  shutil.rmtree('wave.out')
+
+  #analyse mgmol standard output
+  lines=output2.split(b'\n')
+
+  print("Check energy conservation...")
+  tol = 1.e-4
+  energy = 0.
+  count = 0
+  for line in lines:
+    if line.count(b'%%'):
+      print(line)
+    if line.count(b'Total') and line.count(b'Energy'):
+      print(line)
+      count=count+1
+      words=line.split()
+    
+      energy=eval(words[2])
+      if count==1:
+        first_energy=energy
+
+      if count>1 and abs(energy-first_energy)>tol:
+        print("ERROR Energy = {} != {}".format(energy,first_energy))
+        sys.exit(1)
+
+#remove last restart files
+command = "ls -ld snapshot0* | awk '{ print $9 }' | tail -n1"
+restart_file = subprocess.check_output(command,shell=True)
+restart_file=str(restart_file[:-1],'utf-8')
+shutil.rmtree(restart_file)
+
+sys.exit(0)