Add support for reading Desmond DTR (Anton) trajectories (#796)

* DRR - Add the DTR plugin from VMD

* DRR - Comment out vmd-specific includes

* DRR - Rename since likely will want dms support as well

* DRR - DMS plugin from vmd

* DRR - Update README

* DRR - Plugin is really from VMD

* DRR - Add more files from VMD plugin framework to get dtr plugin to comple properly without too much modification

* DRR - Restore original includes except for the byte swap routines; use the cpptraj versions for those

* DRR - Update README

* DRR - Initial version of TrajectoryIO class for DTR

* DRR - Update dependencies

* DRR - Function to determine file or directory

* DRR - Try adding automatic ID

* DRR - Add DTR to TrajectoryFile

* DRR - Start adding memory allocation, deallocation. Use forward declare

* DRR - Start trajectory setup

* DRR - Try to return number of frames

* DRR - Use version of FindDepend from clusterrevamp that can handle different directory levels

* DRR - Follow cpptraj naming convention

* DRR - Update dependencies

* DRR - Fix init check, add some debug info

* DRR - Make it easier to set box info from 6 floats. Finish coord and box read.

* DRR - Refactor reading of mol2 files for topologies since e.g. sometimes VMD writes mol2s without type names

* DRR - Hide debug info

* DRR - Add Time

* DRR - Put the VMD plugin interface behind ifdefs

* DRR - Revision bump for addition of dtr

* DRR - Recognize file name clickme.dtr as desmond DTR for compat. with vmd behavior

* DRR - Add DTR to read traj array

* DRR - Add DTR to trajectory table

* DRR - Actually read the first frame to get box info instead of guessing

* DRR - Add a sanity check

* DRR - Add note about VMD dtr plugin to README

* DRR - Try to fix the windows build by changing ifdef

* DRR - Trouble compiling DTR plugin on windows, so disable for now.

* DRR - Change from NO_DTR to ENABLE_DTR; disabled unless configure is used on non-windows platform.

README.md

@@ -220,3 +220,5 @@ External libraries bundled with CPPTRAJ
 * CPPTRAJ uses the [GROMACS TNG](https://github.com/gromacs/tng) library for reading TNG files. See `sec/tng/README` for details.
 
 * The reciprocal part of the PME calculation is handled by the [helPME](https://github.com/andysim/helpme) library by Andy Simmonett.
+
+* Support for reading DTR trajectories uses the VMD DTR plugin.

configure

@@ -781,6 +781,12 @@ SetupFinalFlags() {
     CPPTRAJ_LIB="$CPPTRAJ_LIB -lpthread"
   fi
   LDFLAGS="$ldf $LDFLAGS"
+  # Platform-specific directives
+  if [ "$PLATFORM" = 'windows' ] ; then
+    echo "Warning: DTR trajectory not supported on windows."
+  else
+    DIRECTIVES="$DIRECTIVES -DENABLE_DTR"
+  fi
 }
 
 #-------------------------------------------------------------------------------

doc/cpptraj.lyx

@@ -12931,7 +12931,7 @@ Cpptraj currently understands the following trajectory file formats:
 \begin_layout Standard
 \align center
 \begin_inset Tabular
-<lyxtabular version="3" rows="21" columns="4">
+<lyxtabular version="3" rows="22" columns="4">
 <features tabularvalignment="middle">
 <column alignment="center" valignment="top" width="25text%">
 <column alignment="center" valignment="top" width="25text%">
@@ -13710,6 +13710,44 @@ xyz
 
 \end_layout
 
+\end_inset
+</cell>
+</row>
+<row>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Desmond DTR (Anton)
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+dtr
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+.dtr
+\end_layout
+
+\end_inset
+</cell>
+<cell alignment="center" valignment="top" topline="true" leftline="true" rightline="true" usebox="none">
+\begin_inset Text
+
+\begin_layout Plain Layout
+Read Only
+\end_layout
+
 \end_inset
 </cell>
 </row>

src/Box.cpp

@@ -122,6 +122,16 @@ void Box::SetBox(Matrix_3x3 const& ucell) {
   SetBoxType();
 }
 
+void Box::SetBox(float A, float B, float C, float alpha, float beta, float gamma)
+{
+  box_[0] = A;
+  box_[1] = B;
+  box_[2] = C;
+  box_[3] = alpha;
+  box_[4] = beta;
+  box_[5] = gamma;
+  SetBoxType();
+}
 // Box::SetTruncOct()
 /** Set as truncated octahedron with all lengths set to whatever X is. */
 void Box::SetTruncOct() {

src/Box.h

@@ -23,6 +23,7 @@ class Box {
     void SetBox(const double*);
     void SetBox(const float*);
     void SetBox(Matrix_3x3 const&);
+    void SetBox(float,float,float,float,float,float);
     void SetTruncOct();
     void SetNoBox();
     void SetMissingInfo(const Box&);

src/FileType.cpp

@@ -0,0 +1,18 @@
+#include "FileType.h"
+#include <sys/stat.h> // stat
+
+Cpptraj::File::Type Cpptraj::File::IdType(std::string const& filenameIn) {
+  if (filenameIn.empty()) return T_UNKNOWN;
+
+  struct stat frame_stat;
+  if (stat(filenameIn.c_str(), &frame_stat) == -1) {
+    return T_UNKNOWN;
+  }
+
+  if ( frame_stat.st_mode & S_IFDIR ) {
+    return T_DIR;
+  } else if ( frame_stat.st_mode & S_IFREG ) {
+    return T_FILE;
+  }
+  return T_UNKNOWN;
+}

src/FileType.h

@@ -0,0 +1,11 @@
+#ifndef INC_FILETYPE_H
+#define INC_FILETYPE_H
+#include <string>
+namespace Cpptraj {
+namespace File {
+  enum Type { T_FILE = 0, T_DIR, T_UNKNOWN };
+  /// \return File type
+  Type IdType(std::string const&);
+}
+}
+#endif

src/FindDepend.cpp

@@ -52,10 +52,21 @@ bool IgnoreHeader(const char* headername) {
 void GetDependencies(string const& filename) {
   char buffer[BUFFERSIZE+1];
   char headername[BUFFERSIZE+1];
+  // Determine path
+  string baseName, dirPrefix;
+  size_t found = filename.find_last_of("/");
+  if (found == std::string::npos) {
+    baseName = filename;
+    //dirPrefix_.clear();
+  } else {
+    baseName = filename.substr(found+1);
+    dirPrefix = filename.substr(0, found+1);
+  }
+  //printf("DEBUG: Dir='%s' Base='%s'\n", dirPrefix.c_str(), baseName.c_str());
   // Determine type
   FileType type;
   string ext;
-  size_t found = filename.find_last_of(".");
+  found = filename.find_last_of(".");
   if (found != string::npos)
     ext = filename.substr(found);
 
@@ -107,7 +118,7 @@ void GetDependencies(string const& filename) {
           size_t pos = strlen(headername);
           if (headername[pos-1]=='"') headername[pos-1]='\0';
           if (!IgnoreHeader(headername))
-            depends.insert( string(headername) );
+            depends.insert( dirPrefix + string(headername) );
         } //else
           //printf("\tSkipping system header line: %s", buffer);
       }

src/Frame.h

@@ -106,8 +106,10 @@ class Frame {
     int CrdIdx()                      const { return crdidx_;        }
     /// Set box alpha, beta, and gamma
     inline void SetBoxAngles(const double*);
-    /// Set box
+    /// Set box from another box
     void SetBox( Box const& b ) { box_ = b; }
+    /// Modify box in place
+    Box& SetBox() { return box_; }
     /// Set temperature
     void SetTemperature(double tIn) { T_ = tIn;     }
     /// Set step

src/Mol2File.cpp

@@ -1,6 +1,9 @@
 #include <cstring> //strlen
 #include <cstdio> //sscanf
+#include <cstdlib> // atof
 #include "Mol2File.h"
+#include "Atom.h"
+#include "Residue.h"
 #include "CpptrajStdio.h" // To print debug info
 #include "StringRoutines.h" // RemoveTrailingWhitespace
 
@@ -129,27 +132,75 @@ int Mol2File::Mol2Bond(int& at1, int& at2) {
 // Mol2File::Mol2XYZ()
 int Mol2File::Mol2XYZ(double *X) {
   if ( Gets(linebuffer_, BUF_SIZE) != 0 ) return 1;
-  sscanf(linebuffer_,"%*i %*s %lf %lf %lf",X, X+1, X+2);
+  int nread = sscanf(linebuffer_,"%*i %*s %lf %lf %lf", X, X+1, X+2);
+  if (nread != 3) {
+    mprintf("Warning: When reading Mol2 coordinates, expected 3, got %i\n", nread);
+    return 1;
+  }
   return 0;
 }
 
 // Mol2File::Mol2Atom()
-Atom Mol2File::Mol2Atom() {
-  char mol2name[10], mol2type[10];
-  double mol2q;
+int Mol2File::Mol2Atom(Atom& atom, Residue& res, double* XYZ) {
+  if ( Gets(linebuffer_, BUF_SIZE) != 0 ) return 1;
+
+  static const unsigned int m2w = 32; // Max width for any mol2 field
+  char field1[m2w];
+  char field2[m2w];
+  char field3[m2w];
+  char field4[m2w];
+  char field5[m2w];
+  char field6[m2w];
+  char field7[m2w];
+  char field8[m2w];
+
+  //         1         2 3 4 5          6         7           8
   // atom_id atom_name x y z atom_type [subst_id [subst_name [charge [status_bit]]]]
-  sscanf(linebuffer_, "%*i %s %*f %*f %*f %s %*i %*s %lf", mol2name, mol2type, &mol2q);
-  NameType m2name( mol2name );
-  return Atom( m2name, mol2type, mol2q );
-}
+  int nread = sscanf(linebuffer_,"%*s %s %s %s %s %s %s %s %s",
+                     field1, field2, field3, field4,
+                     field5, field6, field7, field8);
+
+  const char* mol2name = 0;
+  const char* mol2x = 0;
+  const char* mol2y = 0;
+  const char* mol2z = 0;
+  const char* mol2type = 0;
+  const char* mol2resnum = 0;
+  const char* mol2resname = 0;
+  const char* mol2charge = 0;
 
-// Mol2File::Mol2Residue()
-Residue Mol2File::Mol2Residue() {
-  char resname[10];
-  int current_res;
-  sscanf(linebuffer_,"%*i %*s %*f %*f %*f %*s %i %s", &current_res, resname);
-  NameType rname( resname );
-  return Residue(rname, current_res, ' ', ' ');
+  if (nread == 8) {
+    mol2name = field1;
+    mol2x = field2;
+    mol2y = field3;
+    mol2z = field4;
+    mol2type = field5;
+    mol2resnum = field6;
+    mol2resname = field7;
+    mol2charge = field8;
+  } else if (nread == 7) {
+    // This is a hack for VMD files that get written with no atom type.
+    mol2name = field1;
+    mol2x = field2;
+    mol2y = field3;
+    mol2z = field4;
+    mol2type = field1;
+    mol2resnum = field5;
+    mol2resname = field6;
+    mol2charge = field7;
+  } else {
+    mprinterr("Error: Malformed mol2 line: %s\n", linebuffer_);
+    return 1;
+  }
+
+  atom = Atom( NameType(mol2name), NameType(mol2type), atof(mol2charge) );
+  res = Residue( NameType(mol2resname), atoi(mol2resnum), ' ', ' ' );
+
+  XYZ[0] = atof( mol2x );
+  XYZ[1] = atof( mol2y );
+  XYZ[2] = atof( mol2z );
+
+  return 0;
 }
 
 void Mol2File::WriteMol2Atom(int atnum, Atom const& atomIn,

src/Mol2File.h

@@ -1,9 +1,10 @@
 #ifndef INC_MOL2FILE_H
 #define INC_MOL2FILE_H
 #include "CpptrajFile.h"
-#include "Atom.h"
-#include "Residue.h"
+#include "NameType.h"
 #include <map>
+class Atom;
+class Residue;
 /// Used to access mol2 files.
 class Mol2File : private CpptrajFile {
   public: 
@@ -25,10 +26,9 @@ class Mol2File : private CpptrajFile {
     int Mol2Bond(int &, int &);
     /// Read in the next Mol2 ATOM line. Get the X Y and Z coords.
     int Mol2XYZ(double *);
-    /// Convert current line to Atom 
-    Atom Mol2Atom();
-    /// Convert current line to Residue
-    Residue Mol2Residue();
+    /// Read in next Mol2 ATOM line. Convert to Atom/Residue/Coords 
+    int Mol2Atom(Atom&, Residue&, double*);
+
     /// Write mol2 atom line: at#, atom, res#, res, coords
     void WriteMol2Atom(int, Atom const&, int, const char*, const double*);
     /// Write mol2 bond line; bond#, atom1, atom2

src/Parm_Mol2.cpp

@@ -23,8 +23,10 @@ int Parm_Mol2::ReadParm(FileName const& fname, Topology &parmOut) {
   double XYZ[3];
   Frame Coords;
   for (int atom=0; atom < infile.Mol2Natoms(); atom++) {
-    if ( infile.Mol2XYZ(XYZ) ) return 1;
-    parmOut.AddTopAtom( infile.Mol2Atom(), infile.Mol2Residue() );
+    Atom mol2atom;
+    Residue mol2res;
+    if ( infile.Mol2Atom(mol2atom, mol2res, XYZ) ) return 1;
+    parmOut.AddTopAtom( mol2atom, mol2res );
     Coords.AddXYZ( XYZ );
   }