Make sure position, velocity etc. information is preserved when reinitializeing openmm context

modules/mol/mm/doc/simulation.rst

@@ -201,6 +201,10 @@ mapped back to the attached structure at any time.
   .. method:: AddPositionConstraint(index)
 
     Fixes the position of the particle with index given in the argument.
+    This requires to reinitialize the internal openmm Context 
+    (this is expensive!). Positions, velocities, forces, energies etc. 
+    are properly preserved but e.g. states of random number generators etc.
+    might be lost.
 
     :param index:       Particle to be fixed
     :type index:        :class:`int`
@@ -211,6 +215,10 @@ mapped back to the attached structure at any time.
   .. method:: AddPositionConstraints(indices)
 
     Fixes the position of the atoms with the indices given in in the argument.
+    This requires to reinitialize the internal openmm Context 
+    (this is expensive!). Positions, velocities, forces, energies etc. 
+    are properly preserved but e.g. states of random number generators etc.
+    might be lost.
 
     :param indices:     Particles to be fixed
     :type indices:      :class:`list`
@@ -221,6 +229,10 @@ mapped back to the attached structure at any time.
   .. method:: ResetPositionConstraints()
 
     Removes all position constraints.
+    This requires to reinitialize the internal openmm Context 
+    (this is expensive!). Positions, velocities, forces, energies etc. 
+    are properly preserved but e.g. states of random number generators etc.
+    might be lost.
 
   .. method:: ResetHarmonicBond(index, bond_length, force_constant)
 
@@ -347,7 +359,11 @@ mapped back to the attached structure at any time.
   .. method:: ResetDistanceConstraint(index, constraint_length)
 
     Update of the distance constraint parameters in the simulation **and**
-    in the attached :class:`Topology`
+    in the attached :class:`Topology`. 
+    This requires to reinitialize the internal openmm Context 
+    (this is expensive!). Positions, velocities, forces, energies etc. 
+    are properly preserved but e.g. states of random number generators etc.
+    might be lost. 
     
     :param index:       Distance constraint to be reset
     :param constraint_length: New constraint length in nm
@@ -445,6 +461,10 @@ mapped back to the attached structure at any time.
 
     Update of the mass in the simulation **and**
     in the attached :class:`Topology`
+    This requires to reinitialize the internal openmm Context 
+    (this is expensive!). Positions, velocities, forces, energies etc. 
+    are properly preserved but e.g. states of random number generators etc.
+    might be lost.
 
     :param index:       Mass to be reset
     :param mass:        New mass

modules/mol/mm/src/simulation.cc

@@ -694,7 +694,7 @@ void Simulation::ResetDistanceConstraint(uint index, Real constraint_length){
   int particle1, particle2;
   system_->getConstraintParameters(index,particle1,particle2,dummy);
   system_->setConstraintParameters(index,particle1,particle2,constraint_length);
-  context_->reinitialize();
+  this->ReinitializeContext();
   top_->SetDistanceConstraintParameters(index, constraint_length);
 }
 
@@ -703,7 +703,7 @@ void Simulation::AddPositionConstraint(uint index){
     throw ost::Error("Provided index exceeds number of atoms!");
   }
   system_->setParticleMass(index,0.0);
-  context_->reinitialize();
+  this->ReinitializeContext();
   top_->AddPositionConstraint(index);
 }
 
@@ -717,7 +717,7 @@ void Simulation::AddPositionConstraints(const std::vector<uint>& index){
     system_->setParticleMass(*i,0.0);
     top_->AddPositionConstraint(*i);
   }
-  context_->reinitialize();
+  this->ReinitializeContext();
 }
 
 void Simulation::ResetPositionConstraints(){
@@ -726,7 +726,7 @@ void Simulation::ResetPositionConstraints(){
     system_->setParticleMass(i,original_masses[i]);
   }
   top_->ResetPositionConstraints();
-  context_->reinitialize();
+  this->ReinitializeContext();
 }
 
 void Simulation::ResetHarmonicPositionRestraint(uint index, const geom::Vec3& ref_position, Real k,
@@ -872,7 +872,7 @@ void Simulation::ResetMass(uint index, Real mass){
     throw ost::Error("Provided index exceeds number of atoms!");
   }
   system_->setParticleMass(index,mass);
-  context_->reinitialize();
+  this->ReinitializeContext();
   top_->SetMass(index,mass);
 }
 
@@ -896,5 +896,19 @@ void Simulation::SetPeriodicBoxExtents(geom::Vec3& vec){
   context_->setPeriodicBoxVectors(ucell_a,ucell_b,ucell_c);
 }
 
+void Simulation::ReinitializeContext() {
+  // reinitializing requires to reset all those things!
+  // Be aware, state of random number generators etc might not be
+  // preserved!
+  OpenMM::State state = context_->getState(OpenMM::State::Positions |
+                                           OpenMM::State::Velocities |
+                                           OpenMM::State::Forces |
+                                           OpenMM::State::Energy |
+                                           OpenMM::State::Parameters |
+                                           OpenMM::State::ParameterDerivatives);
+  context_->reinitialize();
+  context_->setState(state); 
+}
+
 }}}
 

modules/mol/mm/src/simulation.hh

@@ -147,6 +147,8 @@ private:
 
   int TimeToNextNotification();
 
+  void ReinitializeContext();
+
   // loads plugins from directory for OpenMM BUT only once per unique path!
   static void EnsurePluginsLoaded(const String& plugin_path);