fixes to dcd io and coord frame

fixed buggy dcd unit cell import
fixed missing dcd unit cell export
added unit cell for dcd io
some code cleanup for dcd export
some code cleanup for CoordFrame and Vec3List

modules/geom/src/composite3_op.cc

@@ -178,76 +178,5 @@ bool IsInSphere(const Sphere& s, const Vec3& v){
   return Length(s.GetOrigin()-v)<=s.GetRadius();
 }
 
-Line3 Vec3List::FitCylinder(const Vec3 initial_direction, const Vec3 center){
-  //This function fits a cylinder to the positions in Vec3List
-  //It takes as argument an initial guess for the direction and the geometric
-  //center of the atoms. The center is not changed during optimisation as the
-  //best fitting cylinder can be shown to have its axis pass through the geometric center
-  Line3 axis=Line3(center,center+initial_direction), axis_old;
-  Real radius,res_sum_old,res_sum,delta_0=0.01,prec=0.0000001,err,norm,delta;
-  unsigned long n_step=1000, n_res=this->size();
-  Vec3 v,gradient;
-  
-  radius=0.0;
-  delta=delta_0;
-  for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-    radius+=geom::Distance(axis,(*i));
-  }
-  radius/=Real(n_res);
-  res_sum=0.0;
-  for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-    res_sum+=pow(Distance(axis,(*i))-radius,2.);
-  }
-  unsigned long k=0;
-  err=2.0*prec;
-  while (err>prec and k<n_step) {
-    res_sum_old=res_sum;
-    axis_old=axis;
-    radius=0.0;
-    if (k>50) {
-      delta=delta_0*pow((50./k),2.0);
-    }
-    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-      radius+=Distance(axis,(*i));
-    }
-    radius/=Real(n_res);
-    for (int j=0; j!=3; ++j){
-      res_sum=0.0;
-      v=Vec3(0.0,0.0,0.0);
-      v[j]=delta;
-      axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()+v);
-      radius=0.0;
-      for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-        radius+=Distance(axis,(*i));
-      }
-      radius/=Real(n_res);
-      for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-        res_sum+=pow(Distance(axis,(*i))-radius,2.);
-      }
-      gradient[j]=(res_sum-res_sum_old)/delta;
-    }
-    norm=Dot(gradient,gradient);
-    if (norm>1.) {
-      gradient=Normalize(gradient);
-    }
-    axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()-delta*gradient);
-    radius=0.0;
-    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-      radius+=Distance(axis,(*i));
-    }
-    radius/=Real(n_res);
-    res_sum=0.0;
-    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
-      res_sum+=pow(Distance(axis,(*i))-radius,2.);
-    }
-    err=fabs((res_sum-res_sum_old)/float(n_res));
-    k++;
-  }
-  if (err>prec) {
-    std::cout<<"axis fitting did not converge"<<std::endl;
-  }
-  return axis;
-}
-
 } // ns
 

modules/geom/src/vec3.cc

@@ -19,8 +19,12 @@
 #include <algorithm>
 
 #include <Eigen/SVD>
-#include "vec3.hh"
 
+#include "vec3.hh"
+// TODO: these are for the (misplaced) Vec3List algorithm functions
+#include "vecmat3_op.hh"
+#include "composite3.hh"
+#include "composite3_op.hh"
 
 namespace geom {
 
@@ -68,18 +72,86 @@ Vec3 Vec3List::GetCenter() const
   return center/=this->size();
 }
 
-Line3 Vec3List::GetODRLine()
+Line3 Vec3List::GetODRLine() const
 {
   Vec3 center=this->GetCenter();
   Vec3 direction=this->GetPrincipalAxes().GetRow(2);
   return Line3(center,center+direction);
 }
 
-Plane Vec3List::GetODRPlane()
+Plane Vec3List::GetODRPlane() const
 {
   Vec3 origin=this->GetCenter();
   Vec3 normal=this->GetPrincipalAxes().GetRow(0);
   return Plane(origin,normal);
 }
   
+Line3 Vec3List::FitCylinder(const Vec3& initial_direction, const Vec3& center) const
+{
+  Line3 axis=Line3(center,center+initial_direction), axis_old;
+  Real radius,res_sum_old,res_sum,delta_0=0.01,prec=0.0000001,err,norm,delta;
+  unsigned long n_step=1000, n_res=this->size();
+  Vec3 v,gradient;
+  
+  radius=0.0;
+  delta=delta_0;
+  for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+    radius+=geom::Distance(axis,(*i));
+  }
+  radius/=Real(n_res);
+  res_sum=0.0;
+  for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+    res_sum+=pow(Distance(axis,(*i))-radius,2.);
+  }
+  unsigned long k=0;
+  err=2.0*prec;
+  while (err>prec and k<n_step) {
+    res_sum_old=res_sum;
+    axis_old=axis;
+    radius=0.0;
+    if (k>50) {
+      delta=delta_0*pow((50./k),2.0);
+    }
+    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+      radius+=Distance(axis,(*i));
+    }
+    radius/=Real(n_res);
+    for (int j=0; j!=3; ++j){
+      res_sum=0.0;
+      v=Vec3(0.0,0.0,0.0);
+      v[j]=delta;
+      axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()+v);
+      radius=0.0;
+      for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+        radius+=Distance(axis,(*i));
+      }
+      radius/=Real(n_res);
+      for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+        res_sum+=pow(Distance(axis,(*i))-radius,2.);
+      }
+      gradient[j]=(res_sum-res_sum_old)/delta;
+    }
+    norm=Dot(gradient,gradient);
+    if (norm>1.) {
+      gradient=Normalize(gradient);
+    }
+    axis=Line3(axis_old.GetOrigin(),axis_old.GetOrigin()+axis_old.GetDirection()-delta*gradient);
+    radius=0.0;
+    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+      radius+=Distance(axis,(*i));
+    }
+    radius/=Real(n_res);
+    res_sum=0.0;
+    for (Vec3List::const_iterator i=this->begin(),e=this->end(); i!=e; ++i) {
+      res_sum+=pow(Distance(axis,(*i))-radius,2.);
+    }
+    err=fabs((res_sum-res_sum_old)/float(n_res));
+    k++;
+  }
+  if (err>prec) {
+    std::cout<<"axis fitting did not converge"<<std::endl;
+  }
+  return axis;
+}
+
 }

modules/geom/src/vec3.hh

@@ -189,15 +189,13 @@ inline std::ostream& operator<<(std::ostream& os, const Vec3& v)
   os << "[" << v.x << ", " << v.y << ", " << v.z << "]";
   return os;
 }
-}
-
-#include <ost/geom/vec2.hh>
-#include <ost/geom/vec4.hh>
-#include <ost/geom/mat3.hh>
-#include <ost/geom/composite3.hh>
+} // ns geom
 
 namespace geom {
 
+  // TODO: move to separate file
+  class Mat3;
+
 class DLLEXPORT_OST_GEOM Vec3List : public std::vector<Vec3> {
 public:
   typedef std::vector<Vec3> base_type;
@@ -213,30 +211,40 @@ public:
     base_type::operator=(rhs);
     return *this;
   }
+
+  // TODO: move some or all of these to stand-alone functions
   Mat3 GetInertia() const;
-  
   Vec3 GetCenter() const;
-  
   Mat3 GetPrincipalAxes() const;
-  Line3 GetODRLine();
-  Plane GetODRPlane();
-  Line3 FitCylinder(const Vec3 initial_direction, const Vec3 center);
-
+  Line3 GetODRLine() const;
+  Plane GetODRPlane() const;
+
+  //This function fits a cylinder to the positions in Vec3List
+  //It takes as argument an initial guess for the direction and the geometric
+  //center of the atoms. The center is not changed during optimisation as the
+  //best fitting cylinder can be shown to have its axis pass through the geometric center
+  Line3 FitCylinder(const Vec3& initial_direction, const Vec3& center) const;
 };
+} // ns geom
 
 
-inline Vec3::Vec3(const Vec2& v): x(v.x), y(v.y), z(0.0) { }
+#include <ost/geom/vec2.hh>
+#include <ost/geom/vec4.hh>
+#include <ost/geom/mat3.hh>
+#include <ost/geom/composite3.hh>
 
-inline Vec3::Vec3(const Vec4& v): x(v.x), y(v.y), z(v.z) 
-{ 
-  if (std::fabs(v.w)<1e-10) {
-    throw DivideByZeroException();
-  }
-  x/=v.w;
-  y/=v.w;
-  z/=v.w;
-}
+namespace geom {
+  inline Vec3::Vec3(const Vec2& v): x(v.x), y(v.y), z(0.0) { }
   
+  inline Vec3::Vec3(const Vec4& v): x(v.x), y(v.y), z(v.z) 
+  { 
+    if (std::fabs(v.w)<1e-10) {
+      throw DivideByZeroException();
+    }
+    x/=v.w;
+    y/=v.w;
+    z/=v.w;
+  }
 } // namespace geom
 
 

modules/io/src/mol/dcd_io.cc

@@ -67,7 +67,7 @@ bool less_index(const mol::AtomHandle& a1, const mol::AtomHandle& a2)
 }
 
 bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag, 
-                     bool& skip_flag, bool& gap_flag)
+                     bool& ucell_flag, bool& gap_flag)
 {
   if (!istream) {
     return false;
@@ -75,7 +75,7 @@ bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag,
   char dummy[4];  
   gap_flag=true;
   swap_flag=false;
-  skip_flag=false;
+  ucell_flag=false;
   if(gap_flag) istream.read(dummy,sizeof(dummy));
   istream.read(header.hdrr,sizeof(char)*4);
   istream.read(reinterpret_cast<char*>(header.icntrl),sizeof(int)*20);
@@ -91,8 +91,8 @@ bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag,
   }
 
   if(header.icntrl[19]!=0) { // CHARMM format
-    skip_flag=(header.icntrl[10]!=0);
-    if(skip_flag) {
+    ucell_flag=(header.icntrl[10]!=0);
+    if(ucell_flag) {
       LOG_VERBOSE("LoadCHARMMTraj: using CHARMM format with per-frame header");
     } else {
       LOG_VERBOSE("LoadCHARMMTraj: using CHARMM format");
@@ -128,10 +128,10 @@ bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag,
 }
 
 
-size_t calc_frame_size(bool skip_flag, bool gap_flag, size_t num_atoms)
+size_t calc_frame_size(bool ucell_flag, bool gap_flag, size_t num_atoms)
 {
   size_t frame_size=0;
-  if (skip_flag) {
+  if (ucell_flag) {
     frame_size+=14*sizeof(int);
   }
   if (gap_flag) {
@@ -143,33 +143,12 @@ 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, 
+                size_t frame_size, bool ucell_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)
 {
   char dummy[4];
-  //if(skip_flag) istream.seekg(14*4,std::ios_base::cur);
-  if(skip_flag){
-    istream.read(dummy,sizeof(dummy));
-    double x,y,z,a,b,c;
-    istream.read(reinterpret_cast<char*>(&x),sizeof(double));
-    istream.read(reinterpret_cast<char*>(&a),sizeof(double));
-    istream.read(reinterpret_cast<char*>(&y),sizeof(double));
-    istream.read(reinterpret_cast<char*>(&b),sizeof(double));
-    istream.read(reinterpret_cast<char*>(&c),sizeof(double));
-    istream.read(reinterpret_cast<char*>(&z),sizeof(double));
-    istream.read(dummy,sizeof(dummy));
-    cell_size[0]=x;
-    cell_size[1]=y;
-    cell_size[2]=z;
-    cell_angles[0]=acos(a);
-    cell_angles[1]=acos(b);
-    cell_angles[2]=acos(c);
-    if(a!=0.0||b!=0.0||c!=0.0){
-       LOG_ERROR("LoadCHARMMTraj: periodic cell not parallelepipedic, cell angles might be wrong, handle carefully")
-    }
-  }
 
   // read each frame
   if(!istream) {
@@ -178,6 +157,21 @@ bool read_frame(std::istream& istream, const DCDHeader& header,
               << frame_num << "Nothing left to be read");
     return false;
   }
+
+  if(ucell_flag){
+    istream.read(dummy,sizeof(dummy));
+    double tmp[6];
+    istream.read(reinterpret_cast<char*>(tmp),sizeof(double)*6);
+    // a,alpha,b,beta,gamma,c (don't ask)
+    cell_size[0]=static_cast<Real>(tmp[0]);
+    cell_size[1]=static_cast<Real>(tmp[2]);
+    cell_size[2]=static_cast<Real>(tmp[5]);
+    cell_angles[0]=static_cast<Real>(acos(tmp[1]));
+    cell_angles[1]=static_cast<Real>(acos(tmp[3]));
+    cell_angles[2]=static_cast<Real>(acos(tmp[4]));
+    istream.read(dummy,sizeof(dummy));
+  }
+
   // x coord
   if(gap_flag) istream.read(dummy,sizeof(dummy));
   istream.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
@@ -224,12 +218,12 @@ 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 ucell_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, 
+  return read_frame(istream,header, frame_size,ucell_flag, gap_flag, 
                     swap_flag, xlist,frame,frame_num, cell_size, cell_angles);
 }
   
@@ -247,8 +241,8 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
   Profile profile_load("LoadCHARMMTraj");
 
   DCDHeader header; 
-  bool swap_flag=false, skip_flag=false, gap_flag=false;
-  read_dcd_header(istream, header, swap_flag, skip_flag, gap_flag);
+  bool swap_flag=false, ucell_flag=false, gap_flag=false;
+  read_dcd_header(istream, header, swap_flag, ucell_flag, gap_flag);
   if(alist.size() != static_cast<size_t>(header.t_atom_count)) {
     LOG_ERROR("LoadCHARMMTraj: atom count missmatch: " << alist.size() 
                << " in coordinate file, " << header.t_atom_count 
@@ -260,18 +254,18 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
   std::vector<geom::Vec3> clist(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(ucell_flag, gap_flag, xlist.size());
   int i=0;
   for(;i<header.num;i+=stride) {
-    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, ucell_flag, gap_flag, 
                     swap_flag, xlist, clist, i,cell_size,cell_angles)) {
       break;
     }
-    if(skip_flag) {
+    if(ucell_flag) {
       cg.AddFrame(clist,cell_size,cell_angles);
+    } else {
+      cg.AddFrame(clist);
     }
-    else cg.AddFrame(clist);
 
     // skip frames (defined by stride)
     if(stride>1) istream.seekg(frame_size*(stride-1),std::ios_base::cur);
@@ -321,7 +315,7 @@ private:
   void FetchFrame(uint frame);
   String               filename_;
   DCDHeader            header_;
-  bool                 skip_flag_;
+  bool                 ucell_flag_;
   bool                 swap_flag_;
   bool                 gap_flag_;
   std::ifstream        stream_;
@@ -337,17 +331,17 @@ private:
 void DCDCoordSource::FetchFrame(uint frame)
 {
   if (!loaded_) {
-    read_dcd_header(stream_, header_, swap_flag_, skip_flag_, gap_flag_);
+    read_dcd_header(stream_, header_, swap_flag_, ucell_flag_, gap_flag_);
     frame_start_=stream_.tellg();
     loaded_=true;
     frame_count_=header_.num/stride_;
   }
-  size_t frame_size=calc_frame_size(skip_flag_, gap_flag_, 
+  size_t frame_size=calc_frame_size(ucell_flag_, gap_flag_, 
                                     header_.t_atom_count);  
   size_t pos=frame_start_+frame_size*frame*stride_;
   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_, 
+  if (!read_frame(stream_, header_, frame_size, ucell_flag_, gap_flag_, 
                   swap_flag_, xlist, *frame_.get(), frame)) {
   }  
 }
@@ -379,48 +373,65 @@ void write_dcd_hdr(std::ofstream& out,
                    const mol::CoordGroupHandle& coord_group, 
                    unsigned int stepsize)
 {
+  // size of first header block in bytes
   int32_t magic_number=84;
-  char crd[]={'C', 'O', 'R', 'D'};
   out.write(reinterpret_cast<char*>(&magic_number), 4);
+
+  // magic string
+  char crd[]={'C', 'O', 'R', 'D'};
   out.write(crd, 4);
+
+  // icntrl[0], NSET, number of frames
   int32_t num_frames=coord_group.GetFrameCount()/stepsize;
   out.write(reinterpret_cast<char*>(&num_frames), 4);
+
+  // icntrl[1], ISTART, starting timestep
   int32_t zero=0;
-  // write zero for istart
   out.write(reinterpret_cast<char*>(&zero), 4);
+
+  // icntrl[2], NSAVC, timesteps between DCD saves
   int32_t one=1;
-  // number of timesteps between dcd saves.
   out.write(reinterpret_cast<char*>(&one), 4);
-  // write spacer of 5 blank integers
-  for (int i=0; i<5; ++i) {
+
+  // icntrl[3] to icntrl[7], unused
+  for (int i=3; i<=7; ++i) {
     out.write(reinterpret_cast<char*>(&zero), 4);
   }
-  // write number of fixed atoms
 
+  // icntrl[8], NAMNF, number of fixed atoms
   out.write(reinterpret_cast<char*>(&zero), 4);
+
+  // icntrl[9], DELTA, timestep as float for CHARMM format
   float delta=1.0;
   out.write(reinterpret_cast<char*>(&delta), 4);
-  // write spacer of 10 blank integers
-  for (int i=0; i <10; ++i) {
+
+  // icntrl[10], CHARMM format: ucell per frame
+  out.write(reinterpret_cast<char*>(&one), 4);
+
+  // icntrl[11] to icntrl[18], unused
+  for (int i=11; i<=18; ++i) {
     out.write(reinterpret_cast<char*>(&zero), 4);
   }
+
+  // icntrl[19], charmm version
+  int32_t charmm_version=24;
+  out.write(reinterpret_cast<char*>(&charmm_version), 4);
+  // bracket first header block
   out.write(reinterpret_cast<char*>(&magic_number), 4);
-  // we don't write any titles for now. This means that the block has only to 
-  // accomodate one int, harr, harr, harr.
-  int32_t title_block_size=4;
-  out.write(reinterpret_cast<char*>(&title_block_size), 4);
-  out.write(reinterpret_cast<char*>(&zero), 4);
+
+  //  no titles in title block
   int32_t four=4;
-  // again block size for titles?
   out.write(reinterpret_cast<char*>(&four), 4);
-  // has to be 4
+  out.write(reinterpret_cast<char*>(&zero), 4);
+  out.write(reinterpret_cast<char*>(&four), 4);
+
+  // atom count block
   out.write(reinterpret_cast<char*>(&four), 4);
   int32_t atom_count=coord_group.GetAtomCount();  
   out.write(reinterpret_cast<char*>(&atom_count), 4);
   out.write(reinterpret_cast<char*>(&four), 4);
 }
-
-}
+} // anon ns
 
 void SaveCHARMMTraj(const mol::CoordGroupHandle& coord_group, 
                     const String& pdb_filename, const String& dcd_filename,
@@ -442,10 +453,26 @@ void SaveCHARMMTraj(const mol::CoordGroupHandle& coord_group,
   int32_t out_n=atom_count*4;
   for (int i=0; i<frame_count; i+=stepsize) {
     mol::CoordFramePtr frame=coord_group.GetFrame(i);
+
+    // ucell
+    int32_t bsize=48; // ucell block size, 6 doubles
+    geom::Vec3 csize=frame->GetCellSize();
+    geom::Vec3 cangles=frame->GetCellAngles();
+    // a,alpha,b,beta,gamma,c (don't ask)
+    double ucell[]={csize[0],
+                    cos(cangles[0]),
+                    csize[1],
+                    cos(cangles[1]),
+                    cos(cangles[2]),
+                    csize[2]};
+    out.write(reinterpret_cast<char*>(&bsize),4);
+    out.write(reinterpret_cast<char*>(ucell),bsize);
+    out.write(reinterpret_cast<char*>(&bsize),4);
+
+
     int k=0;
     for (mol::CoordFrame::iterator j=frame->begin(), 
          e=frame->end(); j!=e; ++j, ++k) {
-      //geom::Vec3 v=*j;
       x[k]=float((*j)[0]);
       y[k]=float((*j)[1]);
       z[k]=float((*j)[2]);

modules/io/tests/CMakeLists.txt

@@ -4,6 +4,7 @@ set(OST_IO_UNIT_TESTS
   test_clustal.cc
   test_io_pdb.cc
   test_io_crd.cc
+  test_io_dcd.cc
   test_io_sdf.cc
   test_pir.cc
   test_iomanager.cc

modules/io/tests/test_io_dcd.cc

+//------------------------------------------------------------------------------
+// This file is part of the OpenStructure project <www.openstructure.org>
+//
+// Copyright (C) 2008-2011 by the OpenStructure authors
+//
+// This library is free software; you can redistribute it and/or modify it under
+// the terms of the GNU Lesser General Public License as published by the Free
+// Software Foundation; either version 3.0 of the License, or (at your option)
+// any later version.
+// This library is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public License for more
+// details.
+//
+// You should have received a copy of the GNU Lesser General Public License
+// along with this library; if not, write to the Free Software Foundation, Inc.,
+// 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+//------------------------------------------------------------------------------
+
+#include <ost/io/mol/dcd_io.hh>
+#include <ost/mol/entity_handle.hh>
+#include <ost/mol/residue_handle.hh>
+#include <ost/mol/chain_handle.hh>
+#include <ost/mol/atom_handle.hh>
+#include <ost/mol/xcs_editor.hh>
+#include <ost/mol/coord_group.hh>
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+using namespace ost;
+using namespace ost::io;
+
+#include <boost/random.hpp>
+namespace {
+  boost::mt19937 RandomGenerator(time(NULL));
+  boost::uniform_01<boost::mt19937> UniformRandom(RandomGenerator);
+}
+
+BOOST_AUTO_TEST_SUITE( io )
+
+BOOST_AUTO_TEST_CASE(test_io_dcd_charmm_frames)
+{
+  mol::EntityHandle eh=mol::CreateEntity();
+  mol::XCSEditor ed=eh.EditXCS();
+  mol::ChainHandle chain=ed.InsertChain("A");
+  mol::ResidueHandle res=ed.AppendResidue(chain,mol::ResidueKey("UNK"));
+
+  static unsigned int natoms=13;
+
+  // create atoms and frame coords at the same time
+  mol::AtomHandleList atoms(natoms);
+  geom::Vec3List atom_pos(natoms);
+  std::ostringstream aname;
+  geom::Vec3 apos1,apos2;
+  for(size_t i=0;i<natoms;++i) {
+    for(size_t k=0;k<3;++k) {
+      apos1[k]=UniformRandom();
+      apos2[k]=UniformRandom();
+    }
+    aname.str("");
+    aname << "X" << i;
+    mol::AtomHandle atom=ed.InsertAtom(res,aname.str(),apos1);
+    atom_pos[i]=apos2;
+    atoms[i]=atom;
+  }
+  geom::Vec3 cell_size(UniformRandom(),UniformRandom(),UniformRandom());
+  geom::Vec3 cell_angles(M_PI*UniformRandom(),M_PI*UniformRandom(),M_PI*UniformRandom());
+
+  mol::CoordGroupHandle cg=mol::CreateCoordGroup(atoms);
+  cg.AddFrame(atom_pos,cell_size,cell_angles);
+
+  SaveCHARMMTraj(cg,"test_io_dcd_out.pdb","test_io_dcd_out.dcd");
+  mol::CoordGroupHandle cg2=LoadCHARMMTraj(eh,"test_io_dcd_out.dcd");
+  BOOST_CHECK_EQUAL(cg2.GetAtomCount(),natoms);
+  BOOST_CHECK_EQUAL(cg2.GetFrameCount(),1);
+
+  mol::CoordFramePtr cf2 = cg2.GetFrame(0);
+  BOOST_CHECK(geom::Distance(cf2->GetCellSize(),cell_size)<1e-6);
+  BOOST_CHECK(geom::Distance(cf2->GetCellAngles(),cell_angles)<1e-6);
+
+  geom::Vec3List atom_pos2=cg2.GetFramePositions(0);
+
+  for(size_t i=0;i<natoms;++i) {
+    BOOST_CHECK(geom::Distance(atom_pos[i],atom_pos2[i])<1e-6);
+  }
+}
+
+BOOST_AUTO_TEST_SUITE_END()

modules/mol/base/pymod/export_coord_frame.cc

@@ -16,6 +16,11 @@
 // along with this library; if not, write to the Free Software Foundation, Inc.,
 // 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 //------------------------------------------------------------------------------
+
+/*
+  Authors: Marco Biasini, Niklaus Johner, Ansgar Philippsen
+*/
+
 #include <boost/python.hpp>
 using namespace boost::python;
 
@@ -26,28 +31,44 @@ using namespace boost::python;
 using namespace ost;
 using namespace ost::mol;
 
-geom::Vec3 (CoordFrame::*GetCellSize)() = &CoordFrame::GetCellSize;
-geom::Vec3 (CoordFrame::*GetCellAngles)() = &CoordFrame::GetCellAngles;
-geom::Vec3 (CoordFrame::*get_atom_pos)(const AtomHandle& atom) = &CoordFrame::GetAtomPos;
-Real (CoordFrame::*get_dist_atom)(const AtomHandle& a1, const AtomHandle& a2) = &CoordFrame::GetDistanceBetwAtoms;
-Real (CoordFrame::*get_angle)(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3) = &CoordFrame::GetAngle;
-Real (CoordFrame::*get_dihedral)(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4) = &CoordFrame::GetDihedralAngle;
-geom::Vec3 (CoordFrame::*get_cm)(const mol::EntityView& sele) = &CoordFrame::GetCenterOfMassPos;
-Real (CoordFrame::*get_dist_cm)(const mol::EntityView& sele1, const mol::EntityView& sele2) = &CoordFrame::GetDistanceBetwCenterOfMass;
-Real (CoordFrame::*get_min_dist_cm_v)(const mol::EntityView& sele1, const mol::EntityView& sele2) = &CoordFrame::GetMinDistBetwCenterOfMassAndView;
-Real (CoordFrame::*get_rmsd)(const mol::EntityView& Reference_View, const mol::EntityView& sele_View) = &CoordFrame::GetRMSD;
-Real (CoordFrame::*get_min_dist)(const mol::EntityView& view1, const mol::EntityView& view2) = &CoordFrame::GetMinDistance;
-Real (CoordFrame::*get_alpha)(const mol::EntityView& segment) = &CoordFrame::GetAlphaHelixContent;
-geom::Line3 (CoordFrame::*get_odr_line)(const mol::EntityView& view1) = &CoordFrame::GetODRLine;
-geom::Line3 (CoordFrame::*get_odr_line2)() = &geom::Vec3List::GetODRLine;
-geom::Plane (CoordFrame::*get_odr_plane)(const mol::EntityView& view1) = &CoordFrame::GetODRPlane;
-geom::Line3 (CoordFrame::*fit_cylinder)(const mol::EntityView& view1) = &CoordFrame::FitCylinder;
+geom::Vec3 (CoordFrame::*get_atom_pos)(const AtomHandle&) const = &CoordFrame::GetAtomPos;
+Real (CoordFrame::*get_dist_atom)(const AtomHandle&, const AtomHandle&) const = &CoordFrame::GetDistanceBetwAtoms;
+Real (CoordFrame::*get_angle)(const AtomHandle&, const AtomHandle&, const AtomHandle&) const = &CoordFrame::GetAngle;
+Real (CoordFrame::*get_dihedral)(const AtomHandle&, const AtomHandle&, const AtomHandle&, const AtomHandle&) const = &CoordFrame::GetDihedralAngle;
+geom::Vec3 (CoordFrame::*get_cm)(const mol::EntityView&) const = &CoordFrame::GetCenterOfMassPos;
+Real (CoordFrame::*get_dist_cm)(const mol::EntityView&, const mol::EntityView&) const = &CoordFrame::GetDistanceBetwCenterOfMass;
+Real (CoordFrame::*get_min_dist_cm_v)(const mol::EntityView&, const mol::EntityView&) const = &CoordFrame::GetMinDistBetwCenterOfMassAndView;
+Real (CoordFrame::*get_rmsd)(const mol::EntityView&, const mol::EntityView&) const = &CoordFrame::GetRMSD;
+Real (CoordFrame::*get_min_dist)(const mol::EntityView&, const mol::EntityView&) const = &CoordFrame::GetMinDistance;
+Real (CoordFrame::*get_alpha)(const mol::EntityView&) const = &CoordFrame::GetAlphaHelixContent;
+geom::Line3 (CoordFrame::*get_odr_line)(const mol::EntityView&) const = &CoordFrame::GetODRLine;
+geom::Plane (CoordFrame::*get_odr_plane)(const mol::EntityView&) const = &CoordFrame::GetODRPlane;
+geom::Line3 (CoordFrame::*fit_cylinder)(const mol::EntityView&) const = &CoordFrame::FitCylinder;
+// TODO: move to geom
+geom::Line3 (CoordFrame::*get_odr_line2)() const = &geom::Vec3List::GetODRLine;
+
+CoordFrame create_coord_frame1(const geom::Vec3List& atom_pos)
+{
+  return CreateCoordFrame(atom_pos);
+}
+
+CoordFrame create_coord_frame2(const geom::Vec3List& atom_pos, 
+                               const geom::Vec3& cell_size,
+                               const geom::Vec3& cell_angles)
+{
+  return CreateCoordFrame(atom_pos,cell_size,cell_angles);
+}
 
 void export_CoordFrame()
 {
-  class_<CoordFrame>("CoordFrame",no_init)
-    .def("GetCellSize",GetCellSize)
-    .def("GetCellAngles",GetCellAngles)
+  // TODO: add ctors or factory functions that take python sequences
+  class_<CoordFrame>("CoordFrame",init<size_t, optional<geom::Vec3> >())
+    .def("SetCellSize",&CoordFrame::SetCellSize)
+    .def("GetCellSize",&CoordFrame::GetCellSize)
+    .add_property("cell_size",&CoordFrame::GetCellSize,&CoordFrame::SetCellSize)
+    .def("SetCellAngles",&CoordFrame::SetCellAngles)
+    .def("GetCellAngles",&CoordFrame::GetCellAngles)
+    .add_property("cell_angles",&CoordFrame::GetCellAngles,&CoordFrame::SetCellAngles)
     .def("GetAtomPos", get_atom_pos)
     .def("GetDistanceBetwAtoms", get_dist_atom)
     .def("GetAngle", get_angle)
@@ -63,5 +84,6 @@ void export_CoordFrame()
     .def("GetAlphaHelixContent",get_alpha)
     .def("FitCylinder",fit_cylinder)
   ;
-  def("CreateCoordFrame",CreateCoordFrame);
+  def("CreateCoordFrame",create_coord_frame1);
+  def("CreateCoordFrame",create_coord_frame2);
 }

modules/mol/base/src/coord_frame.hh

@@ -20,7 +20,7 @@
 #define OST_MOL_COORD_FRAME_HH
 
 /*
-  Author: Marco Biasini and Niklaus Johner
+  Authors: Marco Biasini, Niklaus Johner, Ansgar Philippsen
  */
 #include <boost/shared_ptr.hpp>
 #include <ost/geom/geom.hh>
@@ -30,10 +30,22 @@
 
 namespace ost { namespace mol {
 
+/*!
+  A coordinate frame in a trajectory, containing atom positions and
+  per-frame unit cell size and angles.
+*/
+
+/*
+  TODO: 
+  - move algorithmic code to separate functions in mol/alg
+  - clean up mix of views and atom indices methods/functions
+  - use existing UnitCell class
+*/
+
 class DLLEXPORT_OST_MOL CoordFrame : public geom::Vec3List{
 private:
-  geom::Vec3 periodic_cell_sizes_;
-  geom::Vec3 periodic_cell_angles_;
+  geom::Vec3 ucell_size_;
+  geom::Vec3 ucell_angles_;
 public:
   typedef geom::Vec3List base_type;
   
@@ -42,62 +54,194 @@ public:
   CoordFrame(base_type::iterator b, base_type::iterator e): base_type(b, e) { }
   CoordFrame(const base_type& rhs) : base_type(rhs) { }
   CoordFrame(const std::vector<geom::Vec3>& rhs) : base_type(rhs) { }
-  CoordFrame(const std::vector<geom::Vec3>& rhs,const geom::Vec3 box_size, const geom::Vec3 box_angles) : base_type(rhs) {
-    this->periodic_cell_sizes_=box_size;
-    this->periodic_cell_angles_=box_angles;
+  CoordFrame(const std::vector<geom::Vec3>& rhs,
+             const geom::Vec3 box_size, 
+             const geom::Vec3 box_angles) : base_type(rhs) {
+    ucell_size_=box_size;
+    ucell_angles_=box_angles;
   }
-  
-  geom::Vec3 GetCellSize() {
-    return this->periodic_cell_sizes_;
+
+  void SetCellSize(const geom::Vec3& s) {
+    ucell_size_=s;
   }
-  geom::Vec3 GetCellAngles() {
-    return this->periodic_cell_angles_;
+
+  geom::Vec3 GetCellSize() const {
+    return this->ucell_size_;
   }
-  geom::Vec3 GetAtomPos(const AtomHandle& atom);
-  geom::Vec3 GetAtomPos(int atom_index);
-  Real GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2);
-  Real GetDistanceBetwAtoms(int atom1_index, int atom2_index);
-  Real GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3);
-  Real GetAngle(int atom1_index, int atom2_index, int atom3_index);  
-  Real GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4);
-  Real GetDihedralAngle(int a1_index, int a2_index, int a3_index, int a4_index); 
-  geom::Vec3 GetCenterOfMassPos(const mol::EntityView& sele);
-  geom::Vec3 GetCenterOfMassPos(std::vector<unsigned long>& indices, std::vector<Real>& masses);
-  Real GetDistanceBetwCenterOfMass(const mol::EntityView& sele1, const mol::EntityView& sele2);
+
+  void SetCellAngles(const geom::Vec3& a) {
+    ucell_angles_=a;
+  }
+
+  geom::Vec3 GetCellAngles() const {
+    return this->ucell_angles_;
+  }
+
+  geom::Vec3 GetAtomPos(const AtomHandle& atom) const;
+  geom::Vec3 GetAtomPos(int atom_index) const;
+  Real GetDistanceBetwAtoms(const AtomHandle& a1, const AtomHandle& a2) const;
+  Real GetDistanceBetwAtoms(int atom1_index, int atom2_index) const;
+
+  /*!
+    Returns the angle between the three atoms in rad
+
+    a2 is considered the central atom, i.e the result is the angle between 
+    the vectors (a1.pos-a2.pos) and (a3.pos-a2.pos)
+  */
+  Real GetAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3) const;
+
+  /*!
+    Returns the angle between the three atom indices in rad
+
+    atom2 is considered the central atom, i.e the result is the angle between 
+    the vectors (atom1.pos-atom2.pos) and (atom3.pos-atom2.pos)
+  */
+  Real GetAngle(int atom1_index, int atom2_index, int atom3_index) const;  
+
+  //! Returns the Dihedral angle between the four atoms a1,a2,a3,a4
+  Real GetDihedralAngle(const AtomHandle& a1, const AtomHandle& a2, const AtomHandle& a3, const AtomHandle& a4) const;
+
+  //! Returns the Dihedral angle between the four atom indices a1,a2,a3,a4
+  Real GetDihedralAngle(int a1_index, int a2_index, int a3_index, int a4_index) const; 
+
+  //! Returns the position of the centor of mass of the atoms in the EntityView
+  geom::Vec3 GetCenterOfMassPos(const mol::EntityView& sele) const;
+
+  //! Calculates the center of mass given a list of indices and provided masses
+  /*!
+    the masses vector contains only the masses for the atoms indexed by indices, not for
+    all atoms in the frame
+  */
+  geom::Vec3 GetCenterOfMassPos(std::vector<unsigned long>& indices, 
+                                std::vector<Real>& masses) const;
+
+  //! Returns the distance between the centers of mass of the two EntityViews
+  Real GetDistanceBetwCenterOfMass(const mol::EntityView& sele1, const mol::EntityView& sele2) const;
+
+  /*! Returns the distance betweem two custom center of masses
+    
+    Calculates the centers of mass of the two groups of atoms from given by index lists and
+    custom masses. See comment regarding masses vector for GetCenterOfMassPos
+  */
   Real GetDistanceBetwCenterOfMass(std::vector<unsigned long>& indices1, std::vector<Real>& masses1, 
-                                   std::vector<unsigned long>& indices2, std::vector<Real>& masses2);
-  Real GetRMSD(const std::vector<geom::Vec3>& ref_pos, const std::vector<unsigned long>& indices_sele);
-  Real GetRMSD(const mol::EntityView& reference_view, const mol::EntityView& sele_view);
-  Real GetMinDistance(std::vector<unsigned long>& index_list1, std::vector<unsigned long>& index_list2);
-  Real GetMinDistance(const mol::EntityView& view1, const mol::EntityView& view2);
-  Real GetMinDistBetwCenterOfMassAndView(std::vector<unsigned long>& indices_cm, std::vector<Real>& masses_cm,
-                                         std::vector<unsigned long>& indices_atoms);
-  Real GetMinDistBetwCenterOfMassAndView(const mol::EntityView& view_cm, const mol::EntityView& view_atoms);
-  geom::Line3 GetODRLine(std::vector<unsigned long>& indices_ca);
-  geom::Plane GetODRPlane(std::vector<unsigned long>& indices_ca);
-  geom::Line3 GetODRLine(const mol::EntityView& view1);
-  geom::Plane GetODRPlane(const mol::EntityView& view1);
-  geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca);
-  geom::Line3 FitCylinder(const mol::EntityView& view1);
-  Real GetAlphaHelixContent(std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
-                             std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n);
-  Real GetAlphaHelixContent(const mol::EntityView& segment);
-};
+                                   std::vector<unsigned long>& indices2, std::vector<Real>& masses2) const;
+
+  /*!
+    Returns the RMSD between the positions of the atoms whose indices are given in indices_sele 
+    and the positions given in ref_pos
+  */
+  Real GetRMSD(const std::vector<geom::Vec3>& ref_pos, 
+               const std::vector<unsigned long>& indices_sele) const;
   
+  /*!
+    Returns the rmsd between two EntityViews.
+
+    The reference positions are taken directly from the reference_view,
+    whereas they are taken from this CoordFrame for the sele_view
+  */
+  Real GetRMSD(const mol::EntityView& reference_view,
+               const mol::EntityView& sele_view) const;
+
+  //! Returns the minimal distance between two groups of atoms 
+  /*!
+    with indices given by index_list1 and index_list2
+  */
+  Real GetMinDistance(std::vector<unsigned long>& index_list1, 
+                      std::vector<unsigned long>& index_list2) const;
+
+  //! Returns the minimal distance between the atoms of two views (view1 and view2)
+  Real GetMinDistance(const mol::EntityView& view1, const mol::EntityView& view2) const;
+
+  /*!
+    Returns the minimal distance between the center of mass of two groups of atoms
+    
+    atoms and their masses are defined by indices and masses for those atoms
+  */
+  Real GetMinDistBetwCenterOfMassAndView(std::vector<unsigned long>& indices_cm,
+                                         std::vector<Real>& masses_cm,
+                                         std::vector<unsigned long>& indices_atoms) const;
+  //! Returns the minimal distance between the center of mass of two views
+  Real GetMinDistBetwCenterOfMassAndView(const mol::EntityView& view_cm,
+                                         const mol::EntityView& view_atoms) const;
+
+  //! Returns the best fit line to atoms with indices in indices_ca
+  geom::Line3 GetODRLine(std::vector<unsigned long>& indices_ca) const;
+
+  //! Returns the best fit line to atoms in the EntityView view1
+  geom::Plane GetODRPlane(std::vector<unsigned long>& indices_ca) const;
+
+  //! Returns the normal to the best fit plane to atoms with indices in indices_ca
+  geom::Line3 GetODRLine(const mol::EntityView& view1) const;
+  geom::Plane GetODRPlane(const mol::EntityView& view1) const;
+
+  /*!
+    Returns a line which is the axis of a fitted cylinder to the atoms with indices given in indices_ca
+    It is assumed that we fit an alpha-helix and that the CA atoms are oredered sequentially
+    This is used for the initial guess of the helix axis
+  */
+  geom::Line3 FitCylinder(std::vector<unsigned long>& indices_ca) const;
+
+  //! see FitCylinder(std::vector<unsigned long>&)
+  geom::Line3 FitCylinder(const mol::EntityView& view1) const;
+
+  /*!
+    Returns the percentage of residues in the EntityView (segment) that are in an alpha-helix
+    Each residue is assigned as being in an alpha-helix or not, based on and backbone torsion angles
+    and the hydrogen bonding (O-N) (this term is used only if the segment contains at least 8 residues
+    The entity view should be a single segment with no gaps and no missing N,CA,C,O backbone atoms
+    the first and last residues will not be asessed for helicity
+  */
+  Real GetAlphaHelixContent(std::vector<unsigned long>& indices_ca, 
+                            std::vector<unsigned long>& indices_c,
+                            std::vector<unsigned long>& indices_o,
+                            std::vector<unsigned long>& indices_n) const;
+
+  //! see above
+  Real GetAlphaHelixContent(const mol::EntityView& segment) const;
+};
+
+  typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
+  typedef std::vector<CoordFramePtr> CoordFrameList;
+
+  // factory method
+  // create a frame from a Vec3List containing the positions of the atoms, optionally with unit cell
+  // TODO: why is this necessary? the ctor works fine 
+  DLLEXPORT_OST_MOL CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos, 
+                                                const geom::Vec3& cell_size=geom::Vec3(),
+                                                const geom::Vec3& cell_angles=geom::Vec3());
+
+  // these should really be either in the entity view interface or in mol/alg
+
+  /*!
+    This function returns a vector containing the atom indices of the atoms in an EntityView;
+    it is used to accelerate the extraction of information from a trajectory
+  */
   void GetIndices(const EntityView& sele, std::vector<unsigned long>& indices);
+
+  /*!
+    This function returns a vector containing the atom masses of the atoms in an EntityView;
+    it is used together with GetIndices to accelerate the extraction of RMSD from a trajectory
+  */
   void GetMasses(const EntityView& sele, std::vector<Real>& masses);
-  void GetIndicesAndMasses(const EntityView& sele, std::vector<unsigned long>& indices,std::vector<Real>& masses);
+
+  //! conveniece for GetIndices and GetMasses in one call
+  void GetIndicesAndMasses(const EntityView& sele,
+                           std::vector<unsigned long>& indices,
+                           std::vector<Real>& masses);
+
+  //! Writes the positions of all atoms in the EntityView into the provided vec3 list
   void GetPositions(const EntityView& sele, std::vector<geom::Vec3>& ref_pos);
+
+  //! Writes the indices of all atoms in the EntityView into the provided list
   void GetCaIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca);
-  void GetCaCONIndices(const EntityView& segment, std::vector<unsigned long>& indices_ca, std::vector<unsigned long>& indices_c,
-                      std::vector<unsigned long>& indices_o, std::vector<unsigned long>& indices_n);
 
-typedef boost::shared_ptr<CoordFrame> CoordFramePtr;
-typedef std::vector<CoordFramePtr> CoordFrameList;
+  //! Writes the backbone indices of all residues in the EntityView into the provided list
+  void GetCaCONIndices(const EntityView& segment, 
+                       std::vector<unsigned long>& indices_ca, 
+                       std::vector<unsigned long>& indices_c,
+                       std::vector<unsigned long>& indices_o, 
+                       std::vector<unsigned long>& indices_n);
 
-// factory method
-// create a frame froma Vec3List containing the positions of the atoms
-  DLLEXPORT_OST_MOL CoordFrame CreateCoordFrame(const geom::Vec3List& atom_pos);
 
 }}