Add functionalities for extra info in restart file

src/HDFrestart.cc

@@ -1440,8 +1440,8 @@ int HDFrestart::read_1func_hdf5(T* vv, const std::string& datasetname)
 }
 
 template <class T>
-int HDFrestart::write_1func_hdf5(
-    T* vv, const std::string& datasetname, double* ll, double* cell_origin)
+int HDFrestart::write_1func_hdf5(const T* const vv,
+    const std::string& datasetname, double* ll, double* cell_origin)
 {
     assert(ll != nullptr);
     assert(cell_origin != nullptr);
@@ -1652,7 +1652,7 @@ int HDFrestart::readData(
 }
 
 template <class T>
-int HDFrestart::writeData(T* data, hid_t space_id, hid_t memspace,
+int HDFrestart::writeData(const T* const data, hid_t space_id, hid_t memspace,
     hid_t dset_id, const short precision)
 {
     if (precision == 1)
@@ -2017,7 +2017,8 @@ int HDFrestart::readAtomicData(
 {
     Control& ct = *(Control::instance());
     if (onpe0 && ct.verbose > 0)
-        (*MPIdata::sout) << "HDFrestart::readAtomicNames()..." << std::endl;
+        (*MPIdata::sout) << "HDFrestart::readAtomicData(), dataset = "
+                         << datasetname << std::endl;
 
     std::vector<char> buffer;
     short name_length = 7; // default, value used before February 2016
@@ -2095,10 +2096,9 @@ int HDFrestart::readAtomicData(
     {
         std::string t(&buffer[i], name_length);
         assert(t.size() > 0);
-        // cout<<"name="<<t<<endl;
 
         stripLeadingAndTrailingBlanks(t);
-        // cout<<"stripped name="<<t<<endl;
+        // std::cout<<"stripped name="<<t<<std::endl;
 
         assert(t.size() > 0);
         data.push_back(t);
@@ -2232,13 +2232,13 @@ template int HDFrestart::read_1func_hdf5(float*, const std::string&);
 template int HDFrestart::read_1func_hdf5(double*, const std::string&);
 
 template int HDFrestart::write_1func_hdf5(
-    double*, const std::string&, double* ll, double* origin);
+    const double* const, const std::string&, double* ll, double* origin);
 
 template int HDFrestart::readData(
     double*, hid_t memspace, hid_t dset_id, const short precision);
 template int HDFrestart::readData(
     float*, hid_t memspace, hid_t dset_id, const short precision);
-template int HDFrestart::writeData(double* vv, hid_t filespace, hid_t memspace,
-    hid_t dset_id, const short precision);
-template int HDFrestart::writeData(float* vv, hid_t filespace, hid_t memspace,
-    hid_t dset_id, const short precision);
+template int HDFrestart::writeData(const double* const vv, hid_t filespace,
+    hid_t memspace, hid_t dset_id, const short precision);
+template int HDFrestart::writeData(const float* const vv, hid_t filespace,
+    hid_t memspace, hid_t dset_id, const short precision);

src/HDFrestart.h

@@ -246,16 +246,16 @@ class HDFrestart
     int read_1func_hdf5(T*, const std::string&);
 
     template <class T>
-    int write_1func_hdf5(
-        T*, const std::string&, double* ll = nullptr, double* origin = nullptr);
+    int write_1func_hdf5(const T* const, const std::string&,
+        double* ll = nullptr, double* origin = nullptr);
 
     int read_att(const hid_t dset_id, const std::string& attname,
         std::vector<double>& attr_data);
 
     // write data in file with precision "precision"
     template <class T>
-    int writeData(T* vv, hid_t filespace, hid_t memspace, hid_t dset_id,
-        const short precision);
+    int writeData(const T* const vv, hid_t filespace, hid_t memspace,
+        hid_t dset_id, const short precision);
 
     template <class T>
     int readData(T* vv, hid_t memspace, hid_t dset_id, const short precision);

src/Ion.cc

@@ -294,12 +294,18 @@ void Ion::getIonData(IonData& idata) const
     }
 }
 
+void Ion::resetPositionsToPrevious()
+{
+    for (int pos = 0; pos < 3; pos++)
+        position_[pos] = old_position_[pos];
+}
+
 void Ion::setFromIonData(const IonData& data)
 {
     // random state
     setRandomState(data.rand_state[0], data.rand_state[1], data.rand_state[2]);
     // previous position
-    setOldPosition(
+    setPreviousPosition(
         data.old_position[0], data.old_position[1], data.old_position[2]);
     // velocity
     setVelocity(data.velocity[0], data.velocity[1], data.velocity[2]);

src/Ion.h

@@ -74,12 +74,6 @@ class Ion
 
         position_[i] += shift;
     }
-    void setOldPosition(const double x, const double y, const double z)
-    {
-        old_position_[0] = x;
-        old_position_[1] = y;
-        old_position_[2] = z;
-    }
 
 public:
     Ion(const Species& species, const std::string& name, const double crds[3],
@@ -96,6 +90,13 @@ class Ion
     void init(const double crds[3], const double velocity[3], const bool lock);
     void setup();
 
+    void setPreviousPosition(const double x, const double y, const double z)
+    {
+        old_position_[0] = x;
+        old_position_[1] = y;
+        old_position_[2] = z;
+    }
+
     std::shared_ptr<KBprojector> kbproj() { return kbproj_; }
     const std::shared_ptr<KBprojector> kbproj() const { return kbproj_; }
 
@@ -189,6 +190,9 @@ class Ion
 
         kbproj_->clear();
     }
+
+    void resetPositionsToPrevious();
+
     void shiftPositionXLBOMDTest(Vector3D shift)
     {
         for (short dir = 0; dir < 3; dir++)

src/Ions.cc

@@ -853,6 +853,42 @@ void Ions::writePositions(HDFrestart& h5f_file)
     }
 }
 
+void Ions::writePreviousPositions(HDFrestart& h5f_file)
+{
+    Control& ct(*(Control::instance()));
+
+    if (onpe0 && ct.verbose > 1)
+    {
+        (*MPIdata::sout) << "Ions::writePositions" << std::endl;
+    }
+
+    std::vector<double> data;
+    if (h5f_file.gatherDataX())
+    {
+        Mesh* mymesh             = Mesh::instance();
+        const pb::PEenv& myPEenv = mymesh->peenv();
+        MPI_Comm comm            = myPEenv.comm_x();
+
+        gatherPreviousPositions(data, 0, comm);
+    }
+    else
+    {
+        for (auto& ion : local_ions_)
+        {
+            data.push_back(ion->getPreviousPosition(0));
+            data.push_back(ion->getPreviousPosition(1));
+            data.push_back(ion->getPreviousPosition(2));
+        }
+    }
+
+    hid_t file_id = h5f_file.file_id();
+    if (file_id >= 0)
+    {
+        std::string datasetname("/Ionic_previous_positions");
+        writeData2d(h5f_file, datasetname, data, 3, 1.e32);
+    }
+}
+
 void Ions::initFromRestartFile(HDFrestart& h5_file)
 {
     assert(list_ions_.empty());
@@ -995,14 +1031,39 @@ void Ions::readRestartPositions(HDFrestart& h5_file)
     }
 }
 
+void Ions::readRestartPreviousPositions(HDFrestart& h5_file)
+{
+    Control& ct = *(Control::instance());
+    if (onpe0 && ct.verbose > 0)
+        (*MPIdata::sout) << "Read ionic positions from hdf5 file" << std::endl;
+
+    std::vector<double> data;
+    std::string datasetname("/Ionic_previous_positions");
+    h5_file.readAtomicData(datasetname, data);
+    assert(data.size() == 3 * local_ions_.size());
+
+    int i = 0;
+    for (auto& ion : local_ions_)
+    {
+        ion->setPreviousPosition(data[3 * i], data[3 * i + 1], data[3 * i + 2]);
+        i++;
+    }
+}
+
+void Ions::resetPositionsToPrevious()
+{
+    for (auto& ion : local_ions_)
+    {
+        ion->resetPositionsToPrevious();
+    }
+}
+
 void Ions::writeVelocities(HDFrestart& h5f_file)
 {
     Control& ct(*(Control::instance()));
 
     if (onpe0 && ct.verbose > 1)
-    {
         (*MPIdata::sout) << "Ions::writeVelocities" << std::endl;
-    }
 
     std::vector<double> data;
     if (h5f_file.gatherDataX())
@@ -2535,6 +2596,28 @@ void Ions::gatherPositions(
     if (mype == root) positions = data;
 }
 
+void Ions::gatherPreviousPositions(
+    std::vector<double>& positions, const int root, const MPI_Comm comm) const
+{
+    std::vector<double> local_positions;
+
+    for (auto& ion : local_ions_)
+    {
+        local_positions.push_back(ion->getPreviousPosition(0));
+        local_positions.push_back(ion->getPreviousPosition(1));
+        local_positions.push_back(ion->getPreviousPosition(2));
+    }
+
+    // gather data to PE root
+    std::vector<double> data;
+    mgmol_tools::gatherV(local_positions, data, root, comm);
+
+    int mype = 0;
+    MPI_Comm_rank(comm, &mype);
+    positions.clear();
+    if (mype == root) positions = data;
+}
+
 void Ions::gatherForces(
     std::vector<double>& forces, const int root, const MPI_Comm comm) const
 {

src/Ions.h

@@ -147,6 +147,8 @@ class Ions
         const MPI_Comm comm) const;
     void gatherPositions(std::vector<double>& positions, const int root,
         const MPI_Comm comm) const;
+    void gatherPreviousPositions(std::vector<double>& positions, const int root,
+        const MPI_Comm comm) const;
     void gatherLockedNames(std::vector<std::string>& names, const int root,
         const MPI_Comm comm) const;
     void gatherIndexes(
@@ -198,6 +200,7 @@ class Ions
             ion++;
         }
     }
+    void resetPositionsToPrevious();
     void removeMassCenterMotion();
 
     bool hasNLprojectors()
@@ -237,6 +240,7 @@ class Ions
     double kinetic_E(void) const;
 
     void writePositions(HDFrestart& h5f_file);
+    void writePreviousPositions(HDFrestart& h5f_file);
     void writeVelocities(HDFrestart& h5f_file);
     void writeRandomStates(HDFrestart& h5f_file);
     void writeForces(HDFrestart& h5f_file);
@@ -350,6 +354,7 @@ class Ions
     void addIonToList(const Species& sp, const std::string& name,
         const double crds[3], const double velocity[3], const bool lock);
 
+    void readRestartPreviousPositions(HDFrestart& h5_file);
     // void checkUnicityLocalIons();
 };
 

src/restart.cc

@@ -102,6 +102,8 @@ int MGmol<OrbitalsType>::write_hdf5(HDFrestart& h5f_file,
         ions.writeAtomicNLprojIDs(h5f_file);
         ions.writePositions(h5f_file);
         if (ct.LangevinThermostat()) ions.writeRandomStates(h5f_file);
+        if (ct.AtomsDynamic() == AtomsDynamicType::MD)
+            ions.writePreviousPositions(h5f_file);
         ions.writeVelocities(h5f_file);
         ions.writeForces(h5f_file);