re-organized DCD io

git-svn-id: https://dng.biozentrum.unibas.ch/svn/openstructure/trunk@2035 5a81b35b-ba03-0410-adc8-b2c5c5119f08

modules/io/pymod/wrap_io.cc

@@ -30,7 +30,7 @@ using namespace boost::python;
 #include <ost/io/mol/entity_io_crd_handler.hh>
 #include <ost/io/mol/entity_io_sdf_handler.hh>
 #include <ost/io/mol/pdb_reader.hh>
-#include <ost/io/mol/save_dcd.hh>
+#include <ost/io/mol/dcd_io.hh>
 using namespace ost;
 using namespace ost::io;
 
@@ -65,10 +65,14 @@ BOOST_PYTHON_FUNCTION_OVERLOADS(save_entity_handle_ov,
 BOOST_PYTHON_FUNCTION_OVERLOADS(save_entity_view_ov,
                                 save_ent_view, 2, 3)
 
-BOOST_PYTHON_FUNCTION_OVERLOADS(load_charmm_trj_ov,
-                                LoadCHARMMTraj, 2, 3)
 BOOST_PYTHON_FUNCTION_OVERLOADS(save_charmm_trj_ov,
                                 SaveCHARMMTraj, 3, 4)
+
+mol::CoordGroupHandle load_dcd1(const String& c, const String& t) {return LoadCHARMMTraj(c,t);}
+mol::CoordGroupHandle load_dcd2(const String& c, const String& t, unsigned int s) {return LoadCHARMMTraj(c,t,s);}
+mol::CoordGroupHandle load_dcd3(const mol::EntityHandle& e, const String& t) {return LoadCHARMMTraj(e,t);}
+mol::CoordGroupHandle load_dcd4(const mol::EntityHandle& e, const String& t, unsigned int s) {return LoadCHARMMTraj(e,t,s);}
+
 }
 
 void export_pdb_io();
@@ -111,7 +115,10 @@ BOOST_PYTHON_MODULE(_io)
   def("LoadSDF", &LoadSDF);
 
   def("LoadCRD", &LoadCRD);
-  def("LoadCHARMMTraj",LoadCHARMMTraj,load_charmm_trj_ov());
+  def("LoadCHARMMTraj",load_dcd1);
+  def("LoadCHARMMTraj",load_dcd2);
+  def("LoadCHARMMTraj",load_dcd3);
+  def("LoadCHARMMTraj",load_dcd4);
   def("SaveCHARMMTraj",SaveCHARMMTraj,save_charmm_trj_ov());
   export_pdb_io();
 #if OST_IMG_ENABLED  

modules/io/src/mol/CMakeLists.txt

@@ -8,14 +8,14 @@ save_entity.cc
 load_entity.cc
 surface_io_msms_handler.cc
 load_surface.cc
-save_dcd.cc
+dcd_io.cc
 star_parser.cc
 PARENT_SCOPE
 )
 
 set(OST_IO_MOL_HEADERS
 star_parser.hh
-save_dcd.hh
+dcd_io.hh
 entity_io_crd_handler.hh	
 pdb_io.hh
 entity_io_handler.hh		

modules/io/src/mol/save_dcd.cc → modules/io/src/mol/dcd_io.cc

@@ -17,16 +17,211 @@
 // 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 //------------------------------------------------------------------------------
 
+/*
+  Authors: Ansgar Philippsen, Marco Biasini
+*/
+
 #include <fstream>
+#include <boost/iostreams/filter/gzip.hpp>
+#include <boost/filesystem/convenience.hpp>
 
-#include <ost/io/mol/pdb_writer.hh>
-#include <ost/io/mol/save_dcd.hh>
 #include <ost/stdint.hh>
-/*
-  Author: Marco Biasini
- */
+#include <ost/log.hh>
+#include <ost/profile.hh>
+#include <ost/conop/conop.hh>
+#include <ost/mol/xcs_editor.hh>
+
+#include "dcd_io.hh"
+#include "pdb_reader.hh"
+#include "pdb_writer.hh"
+#include "entity_io_crd_handler.hh"
+#include <ost/io/io_exception.hh>
+#include <ost/io/swap_util.hh>
+
 namespace ost { namespace io {
+
+/* 
+   icntrl[1]: number of coordinate sets in file
+   icntrl[2]: number of previous dynamic steps
+   icntrl[3]: frequency for saving coordinates
+   icntlr[4]: number of steps for creation run
+*/
+
+namespace {
+
+struct DCDHeader {
+  char hdrr[4];
+  int icntrl[20];
+  int ntitle;
+  char title[1024];
+  int num, istep, freq,nstep;
+  int t_atom_count,f_atom_count, atom_count;
+};
+
+mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist,
+			       const String& trj_fn,
+			       unsigned int stride)
+{
+  Profile profile_load("LoadCHARMMTraj");
+
+  bool gap_flag = true;
+
+  boost::filesystem::path trj_f(trj_fn);
+  boost::filesystem::ifstream ff(trj_f);
   
+  DCDHeader header;
+  char dummy[4];
+  bool swap_flag=false;
+
+  LOGN_MESSAGE("importing trajectory data");
+
+  if(gap_flag) ff.read(dummy,sizeof(dummy));
+  ff.read(header.hdrr,sizeof(char)*4);
+  ff.read(reinterpret_cast<char*>(header.icntrl),sizeof(int)*20);
+
+  if(header.icntrl[1]<0 || header.icntrl[1]>1e6) {
+    // nonsense atom count, try swapping
+    swap_int(header.icntrl,20);
+    if(header.icntrl[1]<0 || header.icntrl[1]>1e6) {
+      throw(IOException("LoadCHARMMTraj: nonsense atom count in header"));
+    } else {
+      LOGN_MESSAGE("LoadCHARMMTraj: byte-swapping");
+      swap_flag=true;
+    }
+  }
+
+  bool skip_flag=false;
+
+  if(header.icntrl[19]!=0) { // CHARMM format
+    skip_flag=(header.icntrl[10]!=0);
+    if(skip_flag) {
+      LOGN_VERBOSE("LoadCHARMMTraj: using CHARMM format with per-frame header");
+    } else {
+      LOGN_VERBOSE("LoadCHARMMTraj: using CHARMM format");
+    }
+  } else {
+    // XPLOR format
+    LOGN_VERBOSE("LoadCHARMMTraj: using XPLOR format");
+  }
+
+  if(gap_flag) ff.read(dummy,sizeof(dummy));
+  ff.read(reinterpret_cast<char*>(&header.ntitle),sizeof(int));
+  if(swap_flag) swap_int(&header.ntitle,1);
+  if(gap_flag) ff.read(dummy,sizeof(dummy));
+  ff.read(header.title,sizeof(char)*header.ntitle);
+  header.title[header.ntitle]='\0';
+  if(gap_flag) ff.read(dummy,sizeof(dummy));
+  ff.read(reinterpret_cast<char*>(&header.t_atom_count),sizeof(int));
+  if(swap_flag) swap_int(&header.t_atom_count,1);
+  if(gap_flag) ff.read(dummy,sizeof(dummy));
+  header.num=header.icntrl[0];
+  header.istep=header.icntrl[1];
+  header.freq=header.icntrl[2];
+  header.nstep=header.icntrl[3];
+  header.f_atom_count=header.icntrl[8];
+  header.atom_count=header.t_atom_count-header.f_atom_count;
+
+  LOGN_MESSAGE("LoadCHARMMTraj: " << header.num << " trajectories with " 
+               << header.atom_count << " atoms (" << header.f_atom_count 
+               << " fixed) each");
+
+  if(alist.size() != static_cast<size_t>(header.t_atom_count)) {
+    LOGN_ERROR("LoadCHARMMTraj: atom count missmatch: " << alist.size() 
+               << " in coordinate file, " << header.t_atom_count 
+               << " in each traj frame");
+    throw(IOException("invalid trajectory"));
+  }
+
+  mol::CoordGroupHandle cg=CreateCoordGroup(alist);
+  std::vector<geom::Vec3> clist(header.t_atom_count);
+  std::vector<float> xlist(header.t_atom_count);
+
+  for(int i=0;i<header.num;++i) {
+    if(skip_flag) ff.seekg(14*4,std::ios_base::cur);
+    // read each frame
+    if(!ff) {
+      /* premature EOF */
+      LOGN_ERROR("LoadCHARMMTraj: premature end of file, " << i 
+                 << " frames read");
+      break;
+    }
+    // x coord
+    if(gap_flag) ff.read(dummy,sizeof(dummy));
+    ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
+    if(swap_flag) swap_float(&xlist[0],xlist.size());
+    if(gap_flag) ff.read(dummy,sizeof(dummy));
+    for(uint j=0;j<clist.size();++j) {
+      clist[j].SetX(xlist[j]);
+    }
+    // y coord
+    if(gap_flag) ff.read(dummy,sizeof(dummy));
+    ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
+    if(swap_flag) swap_float(&xlist[0],xlist.size());
+    if(gap_flag) ff.read(dummy,sizeof(dummy));
+    for(uint j=0;j<clist.size();++j) {
+      clist[j].SetY(xlist[j]);
+    }
+    // z coord
+    if(gap_flag) ff.read(dummy,sizeof(dummy));
+    ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
+    if(swap_flag) swap_float(&xlist[0],xlist.size());
+    if(gap_flag) ff.read(dummy,sizeof(dummy));
+    for(uint j=0;j<clist.size();++j) {
+      clist[j].SetZ(xlist[j]);
+    }
+    if(i%stride) {
+      cg.AddFrame(clist);
+    }
+  }
+
+  ff.get();
+  if(!ff.eof()) {
+    LOGN_VERBOSE("LoadCHARMMTraj: unexpected trailing file data, bytes read: " 
+                 << ff.tellg());
+  }
+  return cg;
+}
+
+} // anon ns
+
+mol::CoordGroupHandle LoadCHARMMTraj(const String& crd_fn,
+                                     const String& trj_fn,
+                                     unsigned int stride)
+{
+  boost::filesystem::path crd_f(crd_fn);
+
+  // first the coordinate reference file
+  conop::BuilderP builder = conop::Conopology::Instance().GetBuilder();  
+  mol::EntityHandle ent=mol::CreateEntity();
+  mol::XCSEditor editor=ent.RequestXCSEditor(mol::BUFFERED_EDIT);
+  std::vector<mol::AtomHandle> alist;
+  if(boost::filesystem::extension(crd_f)==".pdb") {
+    PDBReader reader(crd_f);
+    reader.SetFlags(PDB::SEQUENTIAL_ATOM_IMPORT);
+    LOGN_MESSAGE("importing coordinate data");
+    reader.Import(ent);
+    alist = reader.GetSequentialAtoms();
+  } else if (boost::filesystem::extension(crd_f)==".crd") {
+    CRDReader reader(crd_f);
+    LOGN_MESSAGE("importing coordinate data");
+    reader.Import(ent);
+    alist = reader.GetSequentialAtoms();
+  } else {
+    throw(IOException("unsupported coordinate file format"));
+  }
+  conop::Conopology::Instance().ConnectAll(builder,ent);
+
+  return load_dcd(alist,trj_fn,stride);
+}
+
+
+mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const mol::EntityHandle& e,
+						      const String& trj_fn,
+						      unsigned int stride)
+{
+  return load_dcd(e.GetAtomList(),trj_fn,stride);
+}
+
 
 namespace {
 
@@ -82,8 +277,10 @@ void SaveCHARMMTraj(const mol::CoordGroupHandle& coord_group,
                     unsigned int stepsize)
 {  
   if(stepsize==0) stepsize=1;
-  PDBWriter writer(pdb_filename);
-  writer.Write(coord_group.GetAtomList());
+  if(!pdb_filename.empty()) {
+    PDBWriter writer(pdb_filename);
+    writer.Write(coord_group.GetAtomList());
+  }
   std::ofstream out(dcd_filename.c_str(), std::ios::binary);
   write_dcd_hdr(out, coord_group,stepsize);
   int frame_count=coord_group.GetFrameCount();
@@ -115,4 +312,4 @@ void SaveCHARMMTraj(const mol::CoordGroupHandle& coord_group,
   }
 }
 
-}}
+}} // ns

modules/io/src/mol/save_dcd.hh → modules/io/src/mol/dcd_io.hh

@@ -16,23 +16,47 @@
 // along with this library; if not, write to the Free Software Foundation, Inc.,
 // 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 //------------------------------------------------------------------------------
-#ifndef OST_SAVE_DCD_HH
-#define OST_SAVE_DCD_HH
+#ifndef OST_IO_ENTITY_DCD_IO_HH
+#define OST_IO_ENTITY_DCD_IO_HH
 
 /*
-  Author: Marco Biasini
+  CHARMM trajectory IO
+  Authors: Ansgar Philippsen, Marco Biasini
  */
- 
+
 #include <ost/io/module_config.hh>
 #include <ost/mol/coord_group.hh>
 
 namespace ost { namespace io {
-  
-/// \brief export coord group as PDB file and DCD trajectory
+
+/*! \brief import a CHARMM trajectory in dcd format
+    requires the coordinate and the trajectory file; the format
+    of the coordinate file will be automatically deduced from the extension
+    the optional stride parameter will cause only every nth frame to be loaded
+*/
+mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const String& coord,
+						      const String& trj,
+						      unsigned int stride=1);
+
+/*! \brief import a CHARMM trajectory in dcd format with an existing entity
+    requires the existing entity and the trajectory file - obviously the
+    atom layout of the entity must match the trajectory file
+*/
+mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const mol::EntityHandle& e,
+						      const String& trj,
+						      unsigned int stride=1);
+
+
+/*! \brief export coord group as PDB file and DCD trajectory
+    if the pdb filename is an empty string, it won't be exported
+    the optional stride parameter will cause every nth frame to be exported
+ */
 void DLLEXPORT_OST_IO SaveCHARMMTraj(const mol::CoordGroupHandle& coord_group, 
                                      const String& pdb_filename, 
                                      const String& dcd_filename,
-                                     unsigned int stepsize=1);
-}}
+                                     unsigned int stride=1);
+
+
+}} // ns
 
 #endif

modules/io/src/mol/entity_io_crd_handler.cc

@@ -35,10 +35,10 @@
 #include <ost/profile.hh>
 
 #include <ost/io/io_exception.hh>
-#include "entity_io_crd_handler.hh"
-#include "pdb_reader.hh"
 #include <ost/io/swap_util.hh>
 
+#include "entity_io_crd_handler.hh"
+
 namespace ost { namespace io {
 
 using boost::format;
@@ -219,173 +219,6 @@ mol::EntityHandle LoadCRD(const String& file_name)
   return ent;
 }
 
-namespace {
-
-struct TrjHeader {
-  char hdrr[4];
-  int icntrl[20];
-  int ntitle;
-  char title[1024];
-  int num, istep, freq,nstep;
-  int t_atom_count,f_atom_count, atom_count;
-};
-
-}
-
-/* 
-   icntrl[1]: number of coordinate sets in file
-   icntrl[2]: number of previous dynamic steps
-   icntrl[3]: frequency for saving coordinates
-   icntlr[4]: number of steps for creation run
-*/
-
-mol::CoordGroupHandle LoadCHARMMTraj(const String& crd_fn,
-                                     const String& trj_fn,
-                                     int flags)
-{
-  Profile profile_load("LoadCHARMMTraj");
-
-  bool gap_flag = true;
-
-  boost::filesystem::path crd_f(crd_fn);
-  boost::filesystem::path trj_f(trj_fn);
-
-  
-
-  // first the coordinate reference file
-  conop::BuilderP builder = conop::Conopology::Instance().GetBuilder();  
-  mol::EntityHandle ent=mol::CreateEntity();
-  mol::XCSEditor editor=ent.RequestXCSEditor(mol::BUFFERED_EDIT);
-  std::vector<mol::AtomHandle> alist;
-  if(boost::filesystem::extension(crd_f)==".pdb") {
-    PDBReader reader(crd_f);
-    reader.SetFlags(PDB::SEQUENTIAL_ATOM_IMPORT);
-    LOGN_MESSAGE("importing coordinate data");
-    reader.Import(ent);
-    alist = reader.GetSequentialAtoms();
-  } else if (boost::filesystem::extension(crd_f)==".crd") {
-    CRDReader reader(crd_f);
-    LOGN_MESSAGE("importing coordinate data");
-    reader.Import(ent);
-    alist = reader.GetSequentialAtoms();
-  } else {
-    throw(IOException("unsupported coordinate file format"));
-  }
-  conop::Conopology::Instance().ConnectAll(builder,ent);
-
-  // now the trajectory
-  boost::filesystem::ifstream ff(trj_f);
-  
-  TrjHeader header;
-  char dummy[4];
-  bool swap_flag=false;
-
-  LOGN_MESSAGE("importing trajectory data");
-
-  if(gap_flag) ff.read(dummy,sizeof(dummy));
-  ff.read(header.hdrr,sizeof(char)*4);
-  ff.read(reinterpret_cast<char*>(header.icntrl),sizeof(int)*20);
-
-  if(header.icntrl[1]<0 || header.icntrl[1]>1e6) {
-    // nonsense atom count, try swapping
-    swap_int(header.icntrl,20);
-    if(header.icntrl[1]<0 || header.icntrl[1]>1e6) {
-      throw(IOException("LoadCHARMMTraj: nonsense atom count in header"));
-    } else {
-      LOGN_MESSAGE("LoadCHARMMTraj: byte-swapping");
-      swap_flag=true;
-    }
-  }
-
-  bool skip_flag=false;
-
-  if(header.icntrl[19]!=0) { // CHARMM format
-    skip_flag=(header.icntrl[10]!=0);
-    if(skip_flag) {
-      LOGN_VERBOSE("LoadCHARMMTraj: using CHARMM format with per-frame header");
-    } else {
-      LOGN_VERBOSE("LoadCHARMMTraj: using CHARMM format");
-    }
-  } else {
-    // XPLOR format
-    LOGN_VERBOSE("LoadCHARMMTraj: using XPLOR format");
-  }
-
-  if(gap_flag) ff.read(dummy,sizeof(dummy));
-  ff.read(reinterpret_cast<char*>(&header.ntitle),sizeof(int));
-  if(swap_flag) swap_int(&header.ntitle,1);
-  if(gap_flag) ff.read(dummy,sizeof(dummy));
-  ff.read(header.title,sizeof(char)*header.ntitle);
-  header.title[header.ntitle]='\0';
-  if(gap_flag) ff.read(dummy,sizeof(dummy));
-  ff.read(reinterpret_cast<char*>(&header.t_atom_count),sizeof(int));
-  if(swap_flag) swap_int(&header.t_atom_count,1);
-  if(gap_flag) ff.read(dummy,sizeof(dummy));
-  header.num=header.icntrl[0];
-  header.istep=header.icntrl[1];
-  header.freq=header.icntrl[2];
-  header.nstep=header.icntrl[3];
-  header.f_atom_count=header.icntrl[8];
-  header.atom_count=header.t_atom_count-header.f_atom_count;
-
-  LOGN_MESSAGE("LoadCHARMMTraj: " << header.num << " trajectories with " 
-               << header.atom_count << " atoms (" << header.f_atom_count 
-               << " fixed) each");
-
-  if(alist.size() != static_cast<size_t>(header.t_atom_count)) {
-    LOGN_ERROR("LoadCHARMMTraj: atom count missmatch: " << alist.size() 
-               << " in coordinate file, " << header.t_atom_count 
-               << " in each traj frame");
-    throw(IOException("invalid trajectory"));
-  }
-
-  mol::CoordGroupHandle cg=CreateCoordGroup(alist);
-  std::vector<geom::Vec3> clist(header.t_atom_count);
-  std::vector<float> xlist(header.t_atom_count);
-
-  for(int i=0;i<header.num;++i) {
-    if(skip_flag) ff.seekg(14*4,std::ios_base::cur);
-    // read each frame
-    if(!ff) {
-      /* premature EOF */
-      LOGN_ERROR("LoadCHARMMTraj: premature end of file, " << i 
-                 << " frames read");
-      break;
-    }
-    // x coord
-    if(gap_flag) ff.read(dummy,sizeof(dummy));
-    ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
-    if(swap_flag) swap_float(&xlist[0],xlist.size());
-    if(gap_flag) ff.read(dummy,sizeof(dummy));
-    for(uint j=0;j<clist.size();++j) {
-      clist[j].SetX(xlist[j]);
-    }
-    // y coord
-    if(gap_flag) ff.read(dummy,sizeof(dummy));
-    ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
-    if(swap_flag) swap_float(&xlist[0],xlist.size());
-    if(gap_flag) ff.read(dummy,sizeof(dummy));
-    for(uint j=0;j<clist.size();++j) {
-      clist[j].SetY(xlist[j]);
-    }
-    // z coord
-    if(gap_flag) ff.read(dummy,sizeof(dummy));
-    ff.read(reinterpret_cast<char*>(&xlist[0]),sizeof(float)*xlist.size());
-    if(swap_flag) swap_float(&xlist[0],xlist.size());
-    if(gap_flag) ff.read(dummy,sizeof(dummy));
-    for(uint j=0;j<clist.size();++j) {
-      clist[j].SetZ(xlist[j]);
-    }
-    cg.AddFrame(clist);
-  }
-
-  ff.get();
-  if(!ff.eof()) {
-    LOGN_VERBOSE("LoadCHARMMTraj: unexpected trailing file data, bytes read: " 
-                 << ff.tellg());
-  }
-  return cg;
-}
 
 void EntityIOCRDHandler::Import(mol::EntityHandle& ent, 
                                 std::istream& stream)

modules/io/src/mol/entity_io_crd_handler.hh

@@ -25,8 +25,6 @@
  */
 #include <ost/io/mol/entity_io_handler.hh>
 
-#include <ost/mol/coord_group.hh>
-
 #include <boost/iostreams/filtering_stream.hpp>
 #include <boost/filesystem/fstream.hpp>
 
@@ -80,14 +78,6 @@ typedef EntityIOHandlerFactory<EntityIOCRDHandler> EntityIOCRDHandlerFactory;
 
 mol::EntityHandle DLLEXPORT_OST_IO LoadCRD(const String& file_name);
 
-/*
-  the flag can tweak the behavior of the importer; right now, the zero bit
-  can be used to switch old format compatibility on
-*/
-mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const String& crd,
-                                                 const String& trj,
-                                                 int flags=0);
-
 }} // ns
 
 #endif