Skip to content
Snippets Groups Projects
Commit 72b3280f authored by Niklaus Johner's avatar Niklaus Johner
Browse files

Saving periodic cell information to CoordFrames when loading a trajectory.

Specifically: modified CoordGroupHandle::load_dcd and io::read_frame
parent 661b7cf3
No related branches found
No related tags found
No related merge requests found
...@@ -141,14 +141,22 @@ size_t calc_frame_size(bool skip_flag, bool gap_flag, size_t num_atoms) ...@@ -141,14 +141,22 @@ size_t calc_frame_size(bool skip_flag, bool gap_flag, size_t num_atoms)
return frame_size; return frame_size;
} }
bool read_frame(std::istream& istream, const DCDHeader& header, bool read_frame(std::istream& istream, const DCDHeader& header,
size_t frame_size, bool skip_flag, bool gap_flag, size_t frame_size, bool skip_flag, bool gap_flag,
bool swap_flag, std::vector<float>& xlist, bool swap_flag, std::vector<float>& xlist,
std::vector<geom::Vec3>& frame, std::vector<geom::Vec3>& frame,
uint frame_num) uint frame_num,geom::Vec3& cell_size,geom::Vec3& cell_angles)
{ {
char dummy[4]; char dummy[4];
if(skip_flag) istream.seekg(14*4,std::ios_base::cur); //if(skip_flag) istream.seekg(14*4,std::ios_base::cur);
if(skip_flag){
istream.read(dummy,sizeof(dummy));
istream.read(reinterpret_cast<char*>(&cell_size[0]),sizeof(float)*3);
istream.read(reinterpret_cast<char*>(&cell_angles[0]),sizeof(float)*3);
istream.read(dummy,sizeof(dummy));
}
// read each frame // read each frame
if(!istream) { if(!istream) {
/* premature EOF */ /* premature EOF */
...@@ -201,7 +209,17 @@ bool read_frame(std::istream& istream, const DCDHeader& header, ...@@ -201,7 +209,17 @@ bool read_frame(std::istream& istream, const DCDHeader& header,
return true; return true;
} }
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,uint frame_num)
{
geom::Vec3 cell_size=geom::Vec3(),cell_angles=geom::Vec3();
return read_frame(istream,header, frame_size,skip_flag, gap_flag,
swap_flag, xlist,frame,frame_num, cell_size, cell_angles);
}
mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom list is already sorted! mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom list is already sorted!
const String& trj_fn, const String& trj_fn,
unsigned int stride) unsigned int stride)
...@@ -227,20 +245,23 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li ...@@ -227,20 +245,23 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
mol::CoordGroupHandle cg=CreateCoordGroup(alist); mol::CoordGroupHandle cg=CreateCoordGroup(alist);
std::vector<geom::Vec3> clist(header.t_atom_count); std::vector<geom::Vec3> clist(header.t_atom_count);
std::vector<float> xlist(header.t_atom_count); std::vector<float> xlist(header.t_atom_count);
geom::Vec3 cell_size, cell_angles;
size_t frame_size=calc_frame_size(skip_flag, gap_flag, xlist.size()); size_t frame_size=calc_frame_size(skip_flag, gap_flag, xlist.size());
int i=0; int i=0;
for(;i<header.num;i+=stride) { for(;i<header.num;i+=stride) {
std::cout << i << " " << header.num << std::endl; std::cout << i << " " << header.num << std::endl;
if (!read_frame(istream, header, frame_size, skip_flag, gap_flag, if (!read_frame(istream, header, frame_size, skip_flag, gap_flag,
swap_flag, xlist, clist, i)) { swap_flag, xlist, clist, i,cell_size,cell_angles)) {
break; break;
} }
cg.AddFrame(clist); if(skip_flag) {
cg.AddFrame(clist,cell_size,cell_angles);
}
else cg.AddFrame(clist);
// skip frames (defined by stride) // skip frames (defined by stride)
if(stride>1) istream.seekg(frame_size*(stride-1),std::ios_base::cur); if(stride>1) istream.seekg(frame_size*(stride-1),std::ios_base::cur);
} }
istream.get(); istream.get();
if(!istream.eof()) { if(!istream.eof()) {
LOG_VERBOSE("LoadCHARMMTraj: unexpected trailing file data, bytes read: " LOG_VERBOSE("LoadCHARMMTraj: unexpected trailing file data, bytes read: "
...@@ -279,6 +300,7 @@ public: ...@@ -279,6 +300,7 @@ public:
} }
virtual void AddFrame(const std::vector<geom::Vec3>& coords) {} virtual void AddFrame(const std::vector<geom::Vec3>& coords) {}
virtual void AddFrame(const std::vector<geom::Vec3>& coords,const geom::Vec3& box_size,const geom::Vec3& box_angles) {}
virtual void InsertFrame(int pos, const std::vector<geom::Vec3>& coords) {} virtual void InsertFrame(int pos, const std::vector<geom::Vec3>& coords) {}
private: private:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment