fixes to dcd import byte swapping

modules/io/pymod/__init__.py

@@ -109,7 +109,7 @@ def LoadPDB(filename, restrict_chains="", no_hetatms=None,
   :param join_spread_atom_records: If set, overrides the value of 
       :attr:`IOProfile.join_spread_atom_records`.
   
-  :param remote: If set to true, the method tries to load the pdb from the 
+  :param remote: If set to True, the method tries to load the pdb from the 
      remote pdb repository www.pdb.org. The filename is then interpreted as the 
      pdb id.
      
@@ -119,7 +119,7 @@ def LoadPDB(filename, restrict_chains="", no_hetatms=None,
   :param dialect: Specifies the particular dialect to use. If set, overrides 
     the value of :attr:`IOProfile.dialect`
 
-  :param seqres: Whether to read SEQRES records. If set to true, the loaded 
+  :param seqres: Whether to read SEQRES records. If set to True, the loaded 
     entity and seqres entry will be returned as a tuple.
 
   :type dialect: :class:`str`
@@ -238,7 +238,7 @@ LoadMapList=LoadImageList
 
 def LoadCHARMMTraj(crd, dcd_file=None, profile='CHARMM',
                    lazy_load=False, stride=1, 
-                   dialect=None):
+                   dialect=None, detect_swap=True,swap_bytes=False):
   """
   Load CHARMM trajectory file.
   
@@ -257,6 +257,10 @@ def LoadCHARMMTraj(crd, dcd_file=None, profile='CHARMM',
       set, overrides the value of the profile
   :param profile: The IO profile to use for loading the PDB file. See 
       :doc:`profile`.
+  :param detect_swap: if True (the default), then automatic detection of endianess
+      is attempted, otherwise the swap_bytes parameter is used
+  :param swap_bytes: is detect_swap is False, this flag determines whether bytes
+      are swapped upon loading or not
   """
   if not isinstance(crd, mol.EntityHandle):
     if dcd_file==None:
@@ -266,7 +270,7 @@ def LoadCHARMMTraj(crd, dcd_file=None, profile='CHARMM',
   else:
     if not dcd_file:
       raise ValueError("No DCD filename given")
-  return LoadCHARMMTraj_(crd, dcd_file, stride, lazy_load)
+  return LoadCHARMMTraj_(crd, dcd_file, stride, lazy_load, detect_swap, swap_bytes)
 
 def LoadMMCIF(filename, restrict_chains="", fault_tolerant=None, calpha_only=None, profile='DEFAULT', remote=False, strict_hydrogens=None, seqres=False, info=False):
   """
@@ -282,7 +286,7 @@ def LoadMMCIF(filename, restrict_chains="", fault_tolerant=None, calpha_only=Non
   :param fault_tolerant: Enable/disable fault-tolerant import. If set, overrides
      the value of :attr:`IOProfile.fault_tolerant`.
   
-  :param remote: If set to true, the method tries to load the pdb from the 
+  :param remote: If set to True, the method tries to load the pdb from the 
      remote pdb repository www.pdb.org. The filename is then interpreted as the 
      pdb id.
      
@@ -291,7 +295,7 @@ def LoadMMCIF(filename, restrict_chains="", fault_tolerant=None, calpha_only=Non
   :param strict_hydrogens: If set, overrides the value of 
      :attr:`IOProfile.strict_hydrogens`.
 
-  :param seqres: Whether to read SEQRES records. If set to true, the loaded 
+  :param seqres: Whether to read SEQRES records. If set to True, the loaded 
     entity and seqres entry will be returned as second item.
 
   :param info: Whether to return an info container with the other output.

modules/io/pymod/wrap_io.cc

@@ -112,7 +112,9 @@ BOOST_PYTHON_MODULE(_ost_io)
 
   def("LoadCRD", &LoadCRD);
   def("LoadCHARMMTraj_", &LoadCHARMMTraj, (arg("ent"), arg("trj_filename"), 
-      arg("stride")=1, arg("lazy_load")=false));
+                                           arg("stride")=1, arg("lazy_load")=false,
+                                           arg("detect_swap")=true,arg("swap_bytes")=false))
+;
   def("LoadMAE", &LoadMAE);
 
   export_pdb_io();

modules/io/src/mol/dcd_io.cc

@@ -67,31 +67,41 @@ bool less_index(const mol::AtomHandle& a1, const mol::AtomHandle& a2)
 }
 
 bool read_dcd_header(std::istream& istream, DCDHeader& header, bool& swap_flag, 
-                     bool& ucell_flag, bool& gap_flag)
+                     bool& ucell_flag, bool& gap_flag, bool detect_swap, bool byte_swap)
 {
   if (!istream) {
     return false;
   }
   char dummy[4];  
   gap_flag=true;
-  swap_flag=false;
+  swap_flag = detect_swap ? false : byte_swap;
+
   ucell_flag=false;
   if(gap_flag) istream.read(dummy,sizeof(dummy));
   istream.read(header.hdrr,sizeof(char)*4);
   if(header.hdrr[0]!='C' || header.hdrr[1]!='O' || header.hdrr[2]!='R' || header.hdrr[3]!='D') {
     throw IOException("LoadCHARMMTraj: missing CORD magic in header");
   }
+
   istream.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) {
-      std::ostringstream msg;
-      msg << "LoadCHARMMTraj: nonsense atom count (" << header.icntrl[1] << ") in header";
-      throw IOException(msg.str());
-    } else {
-      LOG_VERBOSE("LoadCHARMMTraj: byte-swapping");
-      swap_flag=true;
+
+  if(detect_swap) {
+    if(header.icntrl[1]<0 || header.icntrl[1]>1e7) {
+      // nonsense atom count, try byte swapping
+      swap_int(header.icntrl,20);
+      if(header.icntrl[1]<0 || header.icntrl[1]>1e7) {
+        // still weird? swap back and keep fingers crossed
+        LOG_WARNING("LoadCHARMMTraj: byte swap detection failed, trying to continue");
+        swap_int(header.icntrl,20);
+        swap_flag=false;
+      } else {
+        LOG_VERBOSE("LoadCHARMMTraj: byte-swapping (auto-detected)");
+        swap_flag=true;
+      }
+    }
+  } else {
+    if(byte_swap) {
+      LOG_VERBOSE("LoadCHARMMTraj: byte-swapping (manually set)");
     }
   }
 
@@ -241,7 +251,9 @@ bool read_frame(std::istream& istream, const DCDHeader& header,
   
 mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom list is already sorted!
                                const String& trj_fn,
-                               unsigned int stride)
+                               unsigned int stride,
+                               bool detect_swap,
+                               bool byte_swap)
 {
   std::ifstream istream(trj_fn.c_str(), std::ios::binary);
   if(!istream) {
@@ -253,7 +265,7 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
 
   DCDHeader header; 
   bool swap_flag=false, ucell_flag=false, gap_flag=false;
-  read_dcd_header(istream, header, swap_flag, ucell_flag, gap_flag);
+  read_dcd_header(istream, header, swap_flag, ucell_flag, gap_flag, detect_swap, byte_swap);
   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 
@@ -296,9 +308,10 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist, // this atom li
 class  DCDCoordSource : public mol::CoordSource {
 public:
   DCDCoordSource(const mol::AtomHandleList& atoms, const String& filename, 
-                 uint stride): 
+                 uint stride, bool detect_swap, bool byte_swap): 
     mol::CoordSource(atoms), filename_(filename), 
-    stream_(filename.c_str(), std::ios::binary), loaded_(false), stride_(stride)
+    stream_(filename.c_str(), std::ios::binary), loaded_(false), stride_(stride),
+    detect_swap_(detect_swap), byte_swap_(byte_swap)
   {
     frame_count_=0;
     this->SetMutable(false);
@@ -334,6 +347,7 @@ private:
   uint                 frame_count_;
   uint                 curr_frame_;
   uint                 stride_;
+  bool                 detect_swap_, byte_swap_;
   size_t               frame_start_;
   mol::CoordFramePtr   frame_;
 };
@@ -342,7 +356,7 @@ private:
 void DCDCoordSource::FetchFrame(uint frame)
 {
   if (!loaded_) {
-    read_dcd_header(stream_, header_, swap_flag_, ucell_flag_, gap_flag_);
+    read_dcd_header(stream_, header_, swap_flag_, ucell_flag_, gap_flag_, detect_swap_,byte_swap_);
     frame_start_=stream_.tellg();
     loaded_=true;
     frame_count_=header_.num/stride_;
@@ -365,17 +379,20 @@ typedef boost::shared_ptr<DCDCoordSource> DCDCoordSourcePtr;
 
 mol::CoordGroupHandle LoadCHARMMTraj(const mol::EntityHandle& ent,
                                      const String& trj_fn,
-                                     unsigned int stride, bool lazy_load)
+                                     unsigned int stride,
+                                     bool lazy_load,
+                                     bool detect_swap,
+                                     bool byte_swap)
 {
   mol::AtomHandleList alist(ent.GetAtomList());
   std::sort(alist.begin(),alist.end(),less_index);
   if (lazy_load) {
     LOG_VERBOSE("LoadCHARMMTraj: importing with lazy_load=true");
-    DCDCoordSourcePtr source(new DCDCoordSource(alist, trj_fn, stride));
+    DCDCoordSourcePtr source(new DCDCoordSource(alist, trj_fn, stride, detect_swap, byte_swap));
     return mol::CoordGroupHandle(source);
   }
   LOG_VERBOSE("LoadCHARMMTraj: importing with lazy_load=false");  
-  return load_dcd(alist, trj_fn, stride);
+  return load_dcd(alist, trj_fn, stride, detect_swap, byte_swap);
 }
 
 namespace {

modules/io/src/mol/dcd_io.hh

@@ -36,9 +36,11 @@ namespace ost { namespace io {
     atom layout of the entity must match the trajectory file
 */
 mol::CoordGroupHandle DLLEXPORT_OST_IO LoadCHARMMTraj(const mol::EntityHandle& ent,
-                                                       const String& trj_filename,
-                                                       unsigned int stride=1,
-                                                       bool lazy_load=false);
+                                                      const String& trj_filename,
+                                                      unsigned int stride=1,
+                                                      bool lazy_load=false,
+                                                      bool detect_swap=true,
+                                                      bool byte_swap=false);
 
 
 /*! \brief export coord group as PDB file and DCD trajectory