sort before superposing

modules/io/src/mol/dcd_io.cc

@@ -144,14 +144,16 @@ size_t calc_frame_size(bool skip_flag, bool gap_flag, size_t num_atoms)
 bool read_frame(std::istream& istream, const DCDHeader& header, 
                 size_t frame_size, bool skip_flag, bool gap_flag, 
                 bool swap_flag, std::vector<float>& xlist,
-                std::vector<geom::Vec3>& frame)
+                std::vector<geom::Vec3>& frame,
+                uint frame_num)
 {
   char dummy[4];
   if(skip_flag) istream.seekg(14*4,std::ios_base::cur);
   // read each frame
   if(!istream) {
     /* premature EOF */
-    LOG_ERROR("LoadCHARMMTraj: premature end of file, frames read");
+    LOG_ERROR("LoadCHARMMTraj: premature end of file while trying to read frame "
+              << frame_num);
     return false;
   }
   // x coord
@@ -217,7 +219,7 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
   int i=0;
   for(;i<header.num;i+=stride) {
     if (!read_frame(istream, header, frame_size, skip_flag, gap_flag, 
-                    swap_flag, xlist, clist)) {
+                    swap_flag, xlist, clist, i)) {
       break;
     }
     cg.AddFrame(clist);
@@ -296,7 +298,7 @@ void DCDCoordSource::FetchFrame(uint frame)
   stream_.seekg(pos,std::ios_base::beg);
   std::vector<float> xlist(header_.t_atom_count);
   if (!read_frame(stream_, header_, frame_size, skip_flag_, gap_flag_, 
-                  swap_flag_, xlist, *frame_.get())) {
+                  swap_flag_, xlist, *frame_.get(), frame)) {
   }  
 }
 

modules/mol/alg/src/superpose_frames.cc

@@ -26,6 +26,13 @@
 
 namespace ost { namespace mol { namespace alg {
 
+namespace {
+bool less_index(const mol::AtomHandle& a1, const mol::AtomHandle& a2)
+{
+  return a1.GetIndex()<a2.GetIndex();
+}
+
+}
 
 typedef Eigen::Matrix<Real, 3, 3> EMat3;
 typedef Eigen::Matrix<Real, 1, 3> ECVec3;
@@ -76,7 +83,7 @@ inline EMatX3 col_sub(const EMatX3& m, const ECVec3& s)
 }
 
   
-  void AddSuperposedFrame(CoordGroupHandle& superposed, EMatX3& ref_mat,EMatX3& ref_centered,ECVec3& ref_center,CoordFramePtr frame,std::vector<unsigned long>& indices)
+void AddSuperposedFrame(CoordGroupHandle& superposed, EMatX3& ref_mat,EMatX3& ref_centered,ECVec3& ref_center,CoordFramePtr frame,std::vector<unsigned long>& indices)
 {
   // This function superposes and then adds a CoordFrame (frame) to a CoordGroup (superposed).
   // ref_mat, ref_centered and ref_center contain respectively the positions, centered positions and
@@ -144,7 +151,9 @@ CoordGroupHandle SuperposeFrames(CoordGroupHandle& cg, EntityView& sel,
       indices.push_back(i->GetIndex());
     }
   }
-  CoordGroupHandle superposed=CreateCoordGroup(cg.GetEntity().GetAtomList());
+  mol::AtomHandleList alist(cg.GetEntity().GetAtomList());
+  std::sort(alist.begin(), alist.end(),less_index);
+  CoordGroupHandle superposed=CreateCoordGroup(alist);
   EMatX3 ref_mat=EMatX3::Zero(indices.size(), 3);
   EMatX3 ref_centered=EMatX3::Zero(indices.size(), 3);
   ECVec3 ref_center;
@@ -208,7 +217,9 @@ CoordGroupHandle SuperposeFrames(CoordGroupHandle& cg, EntityView& sel,
   if (int(indices.size())!=ref_view.GetAtomCount()){
     throw std::runtime_error("atom counts of the two views are not equal");
   }
-  CoordGroupHandle superposed=CreateCoordGroup(cg.GetEntity().GetAtomList());
+  mol::AtomHandleList alist(cg.GetEntity().GetAtomList());
+  std::sort(alist.begin(), alist.end(),less_index);
+  CoordGroupHandle superposed=CreateCoordGroup(alist);
   EMatX3 ref_mat=EMatX3::Zero(indices.size(), 3);
   EMatX3 ref_centered=EMatX3::Zero(indices.size(), 3);
   ECVec3 ref_center;