From a978edb0371835921a4ad3a52412e6d366c93030 Mon Sep 17 00:00:00 2001
From: Marco Biasini <marco.biasini@unibas.ch>
Date: Tue, 15 Nov 2011 11:05:55 +0100
Subject: [PATCH] sort before superposing
---
modules/io/src/mol/dcd_io.cc | 10 ++++++----
modules/mol/alg/src/superpose_frames.cc | 17 ++++++++++++++---
2 files changed, 20 insertions(+), 7 deletions(-)
diff --git a/modules/io/src/mol/dcd_io.cc b/modules/io/src/mol/dcd_io.cc
index e5f48a93d..cfc001cc3 100644
--- a/modules/io/src/mol/dcd_io.cc
+++ b/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)) {
}
}
diff --git a/modules/mol/alg/src/superpose_frames.cc b/modules/mol/alg/src/superpose_frames.cc
index 465f776e7..85514d6d7 100644
--- a/modules/mol/alg/src/superpose_frames.cc
+++ b/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;
--
GitLab