Fix segmentation faults when reading in trajectories without forces/velocities after reading in trajectories with forces/velocities (#1188)

* Add tanh

* Add hyperbolic sin/cos

* Add a test that segfaults

* Put into separate functions

* Add test for changing velocity information

* Check if the current coordinate info (specifically velocity/force info) differs from the coordinate info output trajectories were set up with. If set up for velocities/forces and no velocity/force info in the current frame, set all zeroes in the current frame.

* Move to Test_VelFrc

* Revert general tests

* Add test with forces

* Use OutputTrajCommon instead of trajio

* Add ensemble test

* No need for separate CoordInfo function, already stored in OutputTrajCommon

* Add function to check ensembles as well

* If initial traj had vel/frc info but this one does not, make sure the velocity/force info is zeroed out.

* Add test trajectory files

* Protect tests in parallel

* Need to modify FrameArray, not FramePtrArray

* Enable the ensemble test

* Save index/value maps for each ensemble so we can tell if the index value mapping is the same between different trajectories.

* Do a better FP comparison

* When replica indices do not map to same values as in the previous ensemble force sorting by value

* Fix frames in parallel ensemble

* Update ensemble test files with good ones from pmemd (with correct replica indices)

* Used the wrong restarts

* Ensemble test

* Protect in parallel

* Handle case where outgoing Frame already has forces and/or velocities allocated

* Record when ensemble changes

* If the input ensemble changes see if the outgoing frame needs to be modified

* Code to gather ensemble values in parallel

* Hide debug info

* Version 7.10.0. Minor version bump for fixing handling remd ensembles with changing vel/frc info. Also added 3 math funcitons.

* Add bytemp to ensemble help

* Add bytemp option to ensemble help

* Need to use NcTest to compare NetCDF files, not DoTest...

* More protections in parallel

* In serial N_THREADS is not set

---------

Co-authored-by: Daniel R. Roe <daniel.roe@nih.gov>

Author: drroe

doc/CpptrajManual.pdf

@@ -1,3 +1,3 @@
 f6f8cb1a79951d80a9d2656fd9c30f55  CpptrajDevelopmentGuide.lyx
-30a19d148a0011fdb02ecc0d2f26e81f  cpptraj.lyx
+8ba86aca8b37cdcb5b01ef05141e5dce  cpptraj.lyx
 5d9b5b5ed47a3ded57b6464df99b3585  CpptrajManual.lyx

doc/DocumentChecksums.txt

@@ -20551,6 +20551,10 @@ ensemble <file0> {[<start>] [<stop> | last] [offset]} | lastframe
 	       bycrdidx | 
 \end_layout
 
+\begin_layout LyX-Code
+           bytemp |
+\end_layout
+
 \begin_layout LyX-Code
            remlog <remlogfile> [nstlim <nstlim> ntwx <ntwx>]}]
 \begin_inset Separator latexpar
@@ -20628,6 +20632,10 @@ last] Frame to stop reading ensemble at;
  This is preferred over sorting via 'remlog'.
 \end_layout
 
+\begin_layout Description
+[bytemp] If specified the ensemble will be sorted by temperature values instead of replica indices.
+\end_layout
+
 \begin_layout Description
 [remlog
 \begin_inset space ~

doc/cpptraj.lyx

@@ -682,6 +682,8 @@ int CpptrajState::RunEnsemble() {
                                      EnsembleParm[0].Nframes() );
       lastPindex = currentTop->Pindex();
     }
+    // Check output trajectories against current coord info
+    ensembleOut_.CheckEnsembleOutCoordInfo(currentCoordInfo, FrameEnsemble);
 #   ifdef TIMER
     setup_time.Stop();
 #   endif
@@ -895,6 +897,8 @@ int CpptrajState::RunParaEnsemble() {
     return 1;
   // Allocate FrameEnsemble here in case preload is needed.
   FrameEnsemble.SetupFrames( NAV.FirstParm()->Atoms(), NAV.EnsCoordInfo() );
+  // Check output trajectories against current coord info
+  ensembleOut_.CheckEnsembleOutCoordInfo(NAV.EnsCoordInfo(), FrameEnsemble);
   // Figure out if any frames need to be preloaded on ranks
   int preload_err = 0;
   if (!TrajComm.Master()) {
@@ -931,6 +935,13 @@ int CpptrajState::RunParaEnsemble() {
       break;
     }
     if (!NAV.CurrentEns()->BadEnsemble()) {
+      // Has the ensemble changed?
+      if (set != my_start && NAV.EnsembleHasChanged()) {
+        if (debug_ > 0)
+          rprintf("DEBUG: Ensemble has changed at set %i\n", set);
+        // Check output trajectories against current input ensemble
+        ensembleOut_.CheckEnsembleOutCoordInfo(NAV.CurrentEns()->EnsembleCoordInfo(), FrameEnsemble);
+      }
       // Since Frame can be modified by actions, save original and use currentFrame
       ActionFrame currentFrame( SortedFrames[0], set );
       if ( currentFrame.Frm().CheckCoordsInvalid() )
@@ -1337,6 +1348,8 @@ int CpptrajState::RunNormal() {
                                  currentSetup.Nframes() );
       lastPindex = currentSetup.Top().Pindex();
     }
+    // Check output trajectories against current coord info
+    trajoutList_.CheckTrajoutCoordInfo(currentSetup.CoordInfo(), TrajFrame);
 #   ifdef TIMER
     setup_time.Stop();
 #   endif

src/CpptrajState.cpp

@@ -10,6 +10,7 @@ Frame DataSet_Coords::AllocateFrame() const {
 void DataSet_Coords::CommonInfo() const {
   if (cInfo_.HasBox()) mprintf(" Box Coords,");
   if (cInfo_.HasVel()) mprintf(" Velocities,");
+  if (cInfo_.HasForce()) mprintf(" Forces,");
   mprintf(" %i atoms", top_.Natom());
 }
 

src/DataSet_Coords.cpp

@@ -111,8 +111,18 @@ void DataSet_Coords_TRJ::GetFrame(int idx, Frame& fIn) {
       Traj_ = trajinList_[ IDX_.CurrentTrajNum() ];
       // NOTE: Currently enforcing all traj have same # atoms, no need to check topology.
       // See if frame needs (re-)allocation.
-      if (previousTraj == 0 || previousTraj->TrajCoordInfo() != Traj_->TrajCoordInfo())
+      if (previousTraj == 0 || previousTraj->TrajCoordInfo() != Traj_->TrajCoordInfo()) {
         readFrame_.SetupFrameV( Traj_->Traj().Parm()->Atoms(), Traj_->TrajCoordInfo() );
+        // If previously had vel/frc info but now do not, zero it out.
+        if (CoordsInfo().HasVel() && !Traj_->TrajCoordInfo().HasVel())
+          readFrame_.AddVelocities(Frame::Darray(readFrame_.size(), 0.0));
+        if (fIn.HasVelocity() && !Traj_->TrajCoordInfo().HasVel())
+          fIn.ZeroVelocities();
+        if (CoordsInfo().HasForce() && !Traj_->TrajCoordInfo().HasForce())
+          readFrame_.AddForces(Frame::Darray(readFrame_.size(), 0.0));
+        if (fIn.HasForce() && !Traj_->TrajCoordInfo().HasForce())
+          fIn.ZeroForces();
+      }
       // Open traj.
       if (Traj_->BeginTraj()) {
         mprinterr("Error: Could not open trajectory %i '%s'\n", IDX_.CurrentTrajNum(),

src/DataSet_Coords_TRJ.cpp

@@ -1,5 +1,6 @@
 #include "EnsembleIn.h"
 #include "CpptrajStdio.h"
+#include "Constants.h" // FNE
 
 EnsembleIn::EnsembleIn() :
   targetType_(ReplicaInfo::NONE), badEnsemble_(0), debug_(0)
@@ -35,6 +36,27 @@ int EnsembleIn::GatherIndices(int* iAddress, std::vector<RemdIdxType>& allIndice
   return 0;
 }
 
+int EnsembleIn::GatherValues(Frame const& frameIn, ReplicaDimArray const& replicaDims,
+                             std::vector<Darray>& allValues,
+                             Parallel::Comm const& commIn)
+{
+  int Ndims = replicaDims.Ndims();
+  std::vector<double> all_values(allValues.size() * Ndims, 0);
+  std::vector<double> dtemp(Ndims);
+  for (int idim = 0; idim != replicaDims.Ndims(); idim++)
+    dtemp[idim] = frameIn.RemdValue(idim, replicaDims);
+  if (commIn.AllGather( &dtemp[0], Ndims, MPI_DOUBLE, &all_values[0])) {
+    rprinterr("Error: Gathering replica values.\n");
+    return 1;
+  }
+  std::vector<double>::const_iterator idx_it = all_values.begin();
+  for (std::vector<Darray>::iterator it = allValues.begin();
+                                     it != allValues.end();
+                                   ++it, idx_it += Ndims)
+    it->assign(idx_it, idx_it + Ndims);
+  return 0;
+}
+
 #ifdef TIMER
 double EnsembleIn::total_mpi_allgather_ = 0.0;
 double EnsembleIn::total_mpi_sendrecv_ = 0.0;
@@ -93,3 +115,92 @@ int EnsembleIn::SetIndicesMap(std::vector<RemdIdxType> const& allIndices) {
   }
   return 0;
 }
+
+int EnsembleIn::SetIdxValMap(std::vector<RemdIdxType> const& allIndices,
+                             std::vector<Darray> const& allValues)
+{
+  if (allIndices.size() != allValues.size()) {
+    rprinterr("Internal Error: EnsembleIn::SetIdxValMap(): allIndices size %zu != allValues size %zu\n",
+              allIndices.size(), allValues.size());
+    return 1;
+  }
+
+  if (debug_ > 0) {
+    mprintf("DEBUG: index/value map:\n");
+    for (unsigned int member = 0; member != allIndices.size(); member++) {
+      mprintf("\t");
+      for (Frame::RemdIdxType::const_iterator idx = allIndices[member].begin();
+                                              idx != allIndices[member].end(); ++idx)
+        mprintf(" %6i", *idx);
+      mprintf(" : ");
+      for (Darray::const_iterator dval = allValues[member].begin();
+                                  dval != allValues[member].end(); ++dval)
+        mprintf(" %8.3f", *dval);
+      mprintf("\n");
+    }
+  }
+
+  typedef std::pair<RemdIdxType, Darray> ivPair;
+  for (unsigned int idx = 0; idx != allIndices.size(); idx++)
+    RemdIdxValMap_.insert( ivPair(allIndices[idx], allValues[idx]) );
+  return 0;
+}
+
+static inline void indexprint(Frame::RemdIdxType const& IDXin) {
+  for (Frame::RemdIdxType::const_iterator idx = IDXin.begin();
+                                          idx != IDXin.end(); ++idx)
+    mprintf(" %6i", *idx);
+}
+
+static inline void valprint(std::vector<double> const& VALin) {
+  for (std::vector<double>::const_iterator dval = VALin.begin();
+                                           dval != VALin.end(); ++dval)
+    mprintf(" %8.3f", *dval);
+}
+
+/// Floating point not equals.
+static inline bool FNE(double v1, double v2) {
+  double delta = v1 - v2;
+  if (delta < 0.0) delta = -delta;
+  return (delta > Constants::SMALL);
+}
+
+// Compare Darrays
+static inline bool DNE(std::vector<double> const& lhs,
+                       std::vector<double> const& rhs)
+{
+  if (lhs.size() != rhs.size()) return true;
+  for (unsigned int idx = 0; idx < lhs.size(); idx++)
+    if (FNE(lhs[idx], rhs[idx])) return true;
+  return false;
+}
+
+int EnsembleIn::CheckIdxValMap(RemdIdxValMapType const& mapIn) const {
+  if (mapIn.size() != RemdIdxValMap_.size()) {
+    mprintf("Warning: REMD index/value map size %zu does not match previous map size %zu\n",
+            RemdIdxValMap_.size(), mapIn.size());
+    return 1;
+  }
+  int match = 0;
+  for (RemdIdxValMapType::const_iterator it = RemdIdxValMap_.begin();
+                                         it != RemdIdxValMap_.end(); ++it)
+  {
+    RemdIdxValMapType::const_iterator jt = mapIn.find( it->first );
+    if (jt == mapIn.end()) {
+      mprintf("Warning: Index");
+      indexprint( it->first );
+      mprintf(" for this ensemble not found in previous ensemble.\n");
+      match = 1;
+    } else if (DNE(jt->second, it->second)) {
+      mprintf("Warning: Index");
+      indexprint( it->first );
+      mprintf(" values ");
+      valprint( it->second );
+      mprintf(" do not match previous map values ");
+      valprint( jt->second );
+      mprintf("\n");
+      match = 1;
+    }
+  }
+  return match;
+}

src/EnsembleIn.cpp

@@ -16,7 +16,12 @@ class FrameArray;
 class CoordinateInfo;
 /// Read in an array of frames at a time.
 class EnsembleIn {
+  protected:
+    typedef Frame::RemdIdxType RemdIdxType;
+    typedef std::vector<double> Darray;
   public:
+    typedef std::map<RemdIdxType, Darray> RemdIdxValMapType;
+
     EnsembleIn();
     virtual ~EnsembleIn() {}
     virtual int SetupEnsembleRead(FileName const&, ArgList&, Topology*) = 0;
@@ -36,19 +41,28 @@ class EnsembleIn {
     InputTrajCommon const& Traj() const { return traj_;        }
     /// \return true if there was a problem reading the ensemble.
     int BadEnsemble()             const { return badEnsemble_; }
+    /// \return 1 if this index/value map does not match the given map
+    int CheckIdxValMap(RemdIdxValMapType const&) const;
+    /// \return index/value map
+    RemdIdxValMapType const& RemdIdxValMap() const { return RemdIdxValMap_; }
+    /// \return The sorting mode
+    ReplicaInfo::TargetType TargetMode() const { return targetType_; }
+    /// Set the debug level
     void SetDebug(int d)                { debug_ = d;          }
     /// Write current replica mapping to STDOUT.
     void PrintReplicaInfo() const;
   protected:
-    typedef Frame::RemdIdxType RemdIdxType;
 
     int SetTemperatureMap(std::vector<double> const&);
     int SetIndicesMap(std::vector<RemdIdxType> const&);
+    int SetIdxValMap(std::vector<RemdIdxType> const&, std::vector<Darray> const&);
     InputTrajCommon& SetTraj() { return traj_; }
     /// For converting temperature to replica index
     ReplicaInfo::Map<double> TemperatureMap_;
     /// For converting indices to replica index
     ReplicaInfo::Map<RemdIdxType> IndicesMap_;
+
+    RemdIdxValMapType RemdIdxValMap_;
     ReplicaInfo::TargetType targetType_; ///< Hold type of REMD frame being searched for.
     int badEnsemble_;         ///< Set to 1 if problem reading ensemble, 0 otherwise.
     int debug_;
@@ -60,6 +74,8 @@ class EnsembleIn {
     Parallel::Comm const& EnsembleComm() const { return Parallel::EnsembleComm(); }
     static int GatherTemperatures(double*, std::vector<double>&, Parallel::Comm const&);
     static int GatherIndices(int*, std::vector<RemdIdxType>&, int, Parallel::Comm const&);
+    static int GatherValues(Frame const&, ReplicaDimArray const&,
+                            std::vector<Darray>&, Parallel::Comm const&);
 #   ifdef TIMER
     Timer mpi_allgather_timer_;
     Timer mpi_sendrecv_timer_;
@@ -68,8 +84,8 @@ class EnsembleIn {
 #   endif
 #   endif
   private:
-    static void PrintReplicaTmap(ReplicaInfo::Map<double> const&);
-    static void PrintReplicaImap(ReplicaInfo::Map<RemdIdxType> const&);
+    //static void PrintReplicaTmap(ReplicaInfo::Map<double> const&);
+    //static void PrintReplicaImap(ReplicaInfo::Map<RemdIdxType> const&);
 
     InputTrajCommon traj_;
 #   ifdef MPI

src/EnsembleIn.h

@@ -23,6 +23,7 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn,
   double remlog_ntwx      = argIn.getKeyDouble("ntwx",   1.0);
   bool no_sort = argIn.hasKey("nosort");
   bool by_crdidx = argIn.hasKey("bycrdidx");
+  bool by_temp = argIn.hasKey("bytemp");
   if (no_sort && by_crdidx) {
     mprinterr("Error: Should only specify either 'nosort' or 'bycrdidx', not both.\n");
     return 1;
@@ -35,6 +36,12 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn,
     mprinterr("Error: Should only specify either 'bycrdidx' or 'remlog', not both.\n");
     return 1;
   }
+  if (by_temp) {
+    if (by_crdidx || no_sort) {
+      mprinterr("Error: 'bycrdidx'/'nosort' should not be specified with 'bytemp'\n");
+      return 1;
+    }
+  }
   // CRDIDXARG: Parse out 'crdidx <indices list>' now so it is not processed
   //            by SetupTrajIO.
   ArgList crdidxarg;
@@ -127,8 +134,14 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn,
       }
     } else {
       // If dimensions are present index by replica indices, otherwise index
-      // by temperature. 
-      if (cInfo_.ReplicaDimensions().Ndims() > 0)
+      // by temperature.
+      if (by_temp) {
+        targetType_ = ReplicaInfo::TEMP;
+        if (cInfo_.ReplicaDimensions().Ndims() > 1) {
+          mprinterr("Error: 'bytemp' not compatible with multiple dimensions.\n");
+          return 1;
+        }
+      } else if (cInfo_.ReplicaDimensions().Ndims() > 0)
         targetType_ = ReplicaInfo::INDICES;
       else
         targetType_ = ReplicaInfo::TEMP;
@@ -147,10 +160,14 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn,
     frameIn.SetupFrameV( Traj().Parm()->Atoms(), cInfo_ );
     std::vector<double> allTemps;
     std::vector<RemdIdxType> allIndices;
+    std::vector<Darray> allValues;
     if (targetType_ == ReplicaInfo::TEMP)
       allTemps.assign(REMDtraj_.size(), -1.0);
-    else if (targetType_ == ReplicaInfo::INDICES)
+    else if (targetType_ == ReplicaInfo::INDICES) {
       allIndices.resize( REMDtraj_.size() );
+      if (cInfo_.UseRemdValues())
+        allValues.resize( REMDtraj_.size(), Darray(cInfo_.ReplicaDimensions().Ndims(), -1.0) );
+    }
 #   ifdef MPI
     int err = 0;
     if (Parallel::TrajComm().Master()) {
@@ -174,6 +191,10 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn,
           if (GatherIndices(frameIn.iAddress(), allIndices, cInfo_.ReplicaDimensions().Ndims(),
                             EnsembleComm()))
             err = 4;
+          if (cInfo_.UseRemdValues()) {
+            if (GatherValues(frameIn, cInfo_.ReplicaDimensions(), allValues, EnsembleComm()))
+              err = 5;
+          }
         }
       }
     }
@@ -196,6 +217,18 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn,
             Parallel::TrajComm().MasterBcast( &((*it)[0]), it->size(), MPI_INT );
           }
         }
+        if (cInfo_.UseRemdValues()) {
+          for (std::vector<Darray>::iterator it = allValues.begin();
+                                             it != allValues.end(); ++it)
+          {
+            if (Parallel::TrajComm().Master())
+              Parallel::TrajComm().MasterBcast( &((*it)[0]), it->size(), MPI_DOUBLE );
+            else {
+              it->resize( cInfo_.ReplicaDimensions().Ndims() );
+              Parallel::TrajComm().MasterBcast( &((*it)[0]), it->size(), MPI_DOUBLE );
+            }
+          }
+        }
       }
     }
 #   else /* MPI */
@@ -205,14 +238,24 @@ int EnsembleIn_Multi::SetupEnsembleRead(FileName const& tnameIn, ArgList& argIn,
       REMDtraj_[member]->closeTraj();
       if (targetType_ == ReplicaInfo::TEMP)
         allTemps[member] = frameIn.Temperature();
-      else if (targetType_ == ReplicaInfo::INDICES)
+      else if (targetType_ == ReplicaInfo::INDICES) {
         allIndices[member] = frameIn.RemdIndices();
+        if (cInfo_.UseRemdValues()) {
+          for (int idim = 0; idim != cInfo_.ReplicaDimensions().Ndims(); idim++) {
+            double dval = frameIn.RemdValue(idim, cInfo_.ReplicaDimensions());
+            allValues[member][idim] = dval;
+          }
+        } // END if UseRemdValues
+      }
     }
 #   endif /* MPI */
     if (targetType_ == ReplicaInfo::TEMP) {
       if (SetTemperatureMap(allTemps)) return 1;
     } else if (targetType_ == ReplicaInfo::INDICES) {
       if (SetIndicesMap(allIndices)) return 1;
+      if (cInfo_.UseRemdValues()) {
+        if (SetIdxValMap(allIndices, allValues)) return 1;
+      }
     }
   }  // Otherwise NONE, no sorting