diff --git a/src/HDFrestart.cc b/src/HDFrestart.cc index 0fac518e..87982290 100644 --- a/src/HDFrestart.cc +++ b/src/HDFrestart.cc @@ -1440,8 +1440,8 @@ int HDFrestart::read_1func_hdf5(T* vv, const std::string& datasetname) } template -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 -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 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="< - 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 int readData(T* vv, hid_t memspace, hid_t dset_id, const short precision); diff --git a/src/Ion.cc b/src/Ion.cc index 33465635..90350477 100644 --- a/src/Ion.cc +++ b/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]); diff --git a/src/Ion.h b/src/Ion.h index ae870820..871205a3 100644 --- a/src/Ion.h +++ b/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 kbproj() { return kbproj_; } const std::shared_ptr 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++) diff --git a/src/Ions.cc b/src/Ions.cc index e93d82ea..09b97d2e 100644 --- a/src/Ions.cc +++ b/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 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 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 data; if (h5f_file.gatherDataX()) @@ -2535,6 +2596,28 @@ void Ions::gatherPositions( if (mype == root) positions = data; } +void Ions::gatherPreviousPositions( + std::vector& positions, const int root, const MPI_Comm comm) const +{ + std::vector 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 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& forces, const int root, const MPI_Comm comm) const { diff --git a/src/Ions.h b/src/Ions.h index ae4f6adb..b03972d1 100644 --- a/src/Ions.h +++ b/src/Ions.h @@ -147,6 +147,8 @@ class Ions const MPI_Comm comm) const; void gatherPositions(std::vector& positions, const int root, const MPI_Comm comm) const; + void gatherPreviousPositions(std::vector& positions, const int root, + const MPI_Comm comm) const; void gatherLockedNames(std::vector& 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(); }; diff --git a/src/restart.cc b/src/restart.cc index ed7f478c..532043af 100644 --- a/src/restart.cc +++ b/src/restart.cc @@ -102,6 +102,8 @@ int MGmol::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);