diff --git a/README.md b/README.md index e416295807..b890bd2487 100644 --- a/README.md +++ b/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. diff --git a/configure b/configure index ae5c6765c1..13bd79e3a2 100755 --- a/configure +++ b/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 } #------------------------------------------------------------------------------- diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 10d4db2630..635a2ee234 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -12931,7 +12931,7 @@ Cpptraj currently understands the following trajectory file formats: \begin_layout Standard \align center \begin_inset Tabular - + @@ -13710,6 +13710,44 @@ xyz \end_layout +\end_inset + + + + +\begin_inset Text + +\begin_layout Plain Layout +Desmond DTR (Anton) +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +dtr +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +.dtr +\end_layout + +\end_inset + + +\begin_inset Text + +\begin_layout Plain Layout +Read Only +\end_layout + \end_inset diff --git a/src/Box.cpp b/src/Box.cpp index 52deef182a..167498d949 100644 --- a/src/Box.cpp +++ b/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() { diff --git a/src/Box.h b/src/Box.h index abe351a51f..b5db28cbcc 100644 --- a/src/Box.h +++ b/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&); diff --git a/src/FileType.cpp b/src/FileType.cpp new file mode 100644 index 0000000000..46182ad96e --- /dev/null +++ b/src/FileType.cpp @@ -0,0 +1,18 @@ +#include "FileType.h" +#include // 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; +} diff --git a/src/FileType.h b/src/FileType.h new file mode 100644 index 0000000000..856932c543 --- /dev/null +++ b/src/FileType.h @@ -0,0 +1,11 @@ +#ifndef INC_FILETYPE_H +#define INC_FILETYPE_H +#include +namespace Cpptraj { +namespace File { + enum Type { T_FILE = 0, T_DIR, T_UNKNOWN }; + /// \return File type + Type IdType(std::string const&); +} +} +#endif diff --git a/src/FindDepend.cpp b/src/FindDepend.cpp index 2b6a4a51eb..5d5fd9ee53 100644 --- a/src/FindDepend.cpp +++ b/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); } diff --git a/src/Frame.h b/src/Frame.h index 05c71604ee..50ebd7d675 100644 --- a/src/Frame.h +++ b/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 diff --git a/src/Mol2File.cpp b/src/Mol2File.cpp index 85445a8a10..3ae3a6b625 100644 --- a/src/Mol2File.cpp +++ b/src/Mol2File.cpp @@ -1,6 +1,9 @@ #include //strlen #include //sscanf +#include // 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", ¤t_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, diff --git a/src/Mol2File.h b/src/Mol2File.h index 2a7274de98..0599952cee 100644 --- a/src/Mol2File.h +++ b/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 +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 diff --git a/src/Parm_Mol2.cpp b/src/Parm_Mol2.cpp index 0202a5d1ba..21a541062b 100644 --- a/src/Parm_Mol2.cpp +++ b/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 ); } diff --git a/src/Traj_DTR.cpp b/src/Traj_DTR.cpp new file mode 100644 index 0000000000..7a50040be1 --- /dev/null +++ b/src/Traj_DTR.cpp @@ -0,0 +1,253 @@ +#ifdef ENABLE_DTR +#include "Traj_DTR.h" +#include "CpptrajStdio.h" +#include "CpptrajFile.h" +#include "Topology.h" +#include "vmdplugin/dtrplugin.hxx" + +using namespace desres::molfile; +/// CONSTRUCTOR +Traj_DTR::Traj_DTR() : + DTR_(0), + fbuffer_(0), + bufsize_(0) +{} + +/// DESTRUCTOR +Traj_DTR::~Traj_DTR() { + if (DTR_ != 0) delete DTR_; + if (fbuffer_ != 0) delete[] fbuffer_; +} + +/** Identify trajectory format. File should be setup for READ */ +bool Traj_DTR::ID_TrajFormat(CpptrajFile& fileIn) { + // For backwards compatibility with VMD + if (fileIn.Filename().Base() == "clickme.dtr") return true; + if (fileIn.OpenFile()) return false; + unsigned char buffer[4]; + if (fileIn.Read(buffer, 4) != 4) return false; + fileIn.CloseFile(); + if (buffer[0] != 'D' || + buffer[1] != 'E' || + buffer[2] != 'S' || + buffer[3] != 'M') + return false; + return true; +} + +/** Print trajectory info to stdout. */ +void Traj_DTR::Info() { + mprintf("is a Desmond DTR trajectory"); +} + +/** Close file. */ +void Traj_DTR::closeTraj() { + +} + +// ----------------------------------------------------------------------------- +/** Open trajectory for reading. */ +int Traj_DTR::openTrajin() { + + return 0; +} + +/** Read help */ +void Traj_DTR::ReadHelp() { + +} + +/** Process read arguments. */ +int Traj_DTR::processReadArgs(ArgList& argIn) { + + return 0; +} + +/** Set up trajectory for reading. + * \return Number of frames in trajectory. + */ +int Traj_DTR::setupTrajin(FileName const& fname, Topology* trajParm) +{ + if (DTR_ != 0) delete DTR_; + DTR_ = 0; + if (fbuffer_ != 0) delete[] fbuffer_; + fbuffer_ = 0; + + std::string initName; + // check fot .stk file + if (StkReader::recognizes(fname.full())) { + DTR_ = new StkReader; + initName = fname.Full(); + } else { + DTR_ = new DtrReader; + // DTR should be initialized with the dir name + initName = fname.DirPrefix(); + } + + if (debug_ > 0) mprintf("DEBUG: initName= %s\n", initName.c_str()); + + if (!DTR_->init( initName.c_str() )) { + mprinterr("Error: DTR init failed.\n"); + delete DTR_; + DTR_ = 0; + return TRAJIN_ERR; + } + + // Check number of atoms + if ( (int)DTR_->natoms != trajParm->Natom() ) { + mprinterr("Error: # of atoms in DTR (%u) != # atoms in associated topology (%i)\n", + DTR_->natoms, trajParm->Natom()); + return TRAJIN_ERR; + } + + ssize_t nframes = DTR_->size(); + + if (debug_ > 0) mprintf("DEBUG: %zd frames.\n", nframes); + if (nframes < 1) { + mprinterr("Error: No frames detected in DTR trajectory.\n"); + return TRAJIN_ERR; + } + + bool has_vel = DTR_->has_velocities(); + // Allocate float buffer + bufsize_ = 3 * (size_t)trajParm->Natom(); + if (has_vel) + bufsize_ *= 2; + fbuffer_ = new float[ bufsize_ ]; + + // Set Coordinate info. + // Read the first frame to get the box info. A bit wasteful but gets + // the job done for now. + // NOTE: DTR seems to always have box? Always orthogonal? + Box tmpBox; + molfile_timestep_t Tstep; + Tstep.coords = fbuffer_; + if (has_vel) + Tstep.velocities = fbuffer_ + bufsize_; + int ret = DTR_->frame(0, &Tstep); + // -1 is EOF + // 0 is success + if ( ret != 0 ) { + mprinterr("Error: Could not read first frame of DTR during setup.\n"); + return 1; + } + + tmpBox.SetBox( Tstep.A, Tstep.B, Tstep.C, + Tstep.alpha, Tstep.beta, Tstep.gamma ); + + SetCoordInfo( CoordinateInfo(tmpBox, has_vel, false, true) ); + + return nframes; +} + +static inline void FloatToDouble(double* darray, const float* farray, size_t asize) +{ + for (size_t idx = 0; idx != asize; idx++) + darray[idx] = (double)farray[idx]; +} + +/** Read specified trajectory frame. */ +int Traj_DTR::readFrame(int set, Frame& frameIn) { + molfile_timestep_t Tstep; + Tstep.coords = fbuffer_; + if (CoordInfo().HasVel()) + Tstep.velocities = fbuffer_ + bufsize_; + + int ret = DTR_->frame(set, &Tstep); + // -1 is EOF + // 0 is success + if ( ret != 0 ) return 1; + // Convert float to double + FloatToDouble(frameIn.xAddress(), Tstep.coords, bufsize_); + if (CoordInfo().HasVel()) + FloatToDouble(frameIn.vAddress(), Tstep.velocities, bufsize_); + + // Set box + frameIn.SetBox().SetBox( Tstep.A, Tstep.B, Tstep.C, + Tstep.alpha, Tstep.beta, Tstep.gamma ); + + // Time + frameIn.SetTime( Tstep.physical_time ); + + return 0; +} + +/** Read velocities from specified frame. */ +int Traj_DTR::readVelocity(int set, Frame& frameIn) { + + return 0; +} + +/** Read forces from specified frame. */ +int Traj_DTR::readForce(int set, Frame& frameIn) { + + return 0; +} + +// ----------------------------------------------------------------------------- +/** Write help. */ +void Traj_DTR::WriteHelp() { + +} + +/** Process write arguments. */ +int Traj_DTR::processWriteArgs(ArgList& argIn, DataSetList const& DSLin) { + + return 0; +} + +/** Set up trajectory for write. */ +int Traj_DTR::setupTrajout(FileName const& fname, Topology* trajParm, + CoordinateInfo const& cInfoIn, + int NframesToWrite, bool append) +{ + + return 1; +} + +/** Write specified trajectory frame. */ +int Traj_DTR::writeFrame(int set, Frame const& frameOut) { + + return 0; +} + +// ============================================================================= +#ifdef MPI +/** Open trajectory for reading in parallel. */ +int Traj_DTR::parallelOpenTrajin(Parallel::Comm const& commIn) { + return 1; +} + +/** Open trajectory for writing in parallel. */ +int Traj_DTR::parallelOpenTrajout(Parallel::Comm const& commIn) { + return 1; +} + +/** Set up trajectory for write in parallel. */ +int Traj_DTR::parallelSetupTrajout(FileName const& fname, Topology* trajParm, + CoordinateInfo const& cInfoIn, + int NframesToWrite, bool append, + Parallel::Comm const& commIn) +{ + + return 1; +} + +/** Read frame in parallel. */ +int Traj_DTR::parallelReadFrame(int set, Frame& frameIn) { + + return 1; +} + +/** Write frame in parallel. */ +int Traj_DTR::parallelWriteFrame(int set, Frame const& frameOut) { + + return 1; +} + +/** Close trajectory in parallel. */ +void Traj_DTR::parallelCloseTraj() { + +} +#endif +#endif /* ENABLE_DTR */ diff --git a/src/Traj_DTR.h b/src/Traj_DTR.h new file mode 100644 index 0000000000..1f91eeb85a --- /dev/null +++ b/src/Traj_DTR.h @@ -0,0 +1,46 @@ +#ifndef INC_TRAJ_DTR_H +#define INC_TRAJ_DTR_H +#ifdef ENABLE_DTR +#include "TrajectoryIO.h" +namespace desres { namespace molfile { class FrameSetReader; } } +/// Used to read Desmond DTR trajectory +class Traj_DTR : public TrajectoryIO { + public: + Traj_DTR(); + ~Traj_DTR(); + static BaseIOtype* Alloc() { return (BaseIOtype*)new Traj_DTR(); } + static void WriteHelp(); + static void ReadHelp(); + private: + // ----- Inherited functions ----------------- + bool ID_TrajFormat(CpptrajFile&); + int setupTrajin(FileName const&, Topology*); + int setupTrajout(FileName const&, Topology*, CoordinateInfo const&,int, bool); + int openTrajin(); + void closeTraj(); + int readFrame(int,Frame&); + int writeFrame(int,Frame const&); + void Info(); + int readVelocity(int, Frame&); + int readForce(int, Frame&); + int processWriteArgs(ArgList&, DataSetList const&); + int processReadArgs(ArgList&); + // ------------------------------------------- +# ifdef MPI + // ----- Parallel functions ------------------ + int parallelOpenTrajin(Parallel::Comm const&); + int parallelOpenTrajout(Parallel::Comm const&); + int parallelSetupTrajout(FileName const&, Topology*, CoordinateInfo const&, + int, bool, Parallel::Comm const&); + int parallelReadFrame(int, Frame&); + int parallelWriteFrame(int, Frame const&); + void parallelCloseTraj(); + // ------------------------------------------- +# endif + + desres::molfile::FrameSetReader* DTR_; + float* fbuffer_; ///< For reading in coords and velocities + size_t bufsize_; ///< Size of buffer_ +}; +#endif /* ENABLE_DTR */ +#endif diff --git a/src/TrajectoryFile.cpp b/src/TrajectoryFile.cpp index bdeace067e..bc56fb82d1 100644 --- a/src/TrajectoryFile.cpp +++ b/src/TrajectoryFile.cpp @@ -25,6 +25,7 @@ #include "Traj_XYZ.h" #include "Traj_GmxTng.h" #include "Traj_GmxDump.h" +#include "Traj_DTR.h" // ----- STATIC VARS / ROUTINES ------------------------------------------------ // NOTE: Must be in same order as TrajFormatType @@ -71,6 +72,11 @@ const FileTypes::AllocToken TrajectoryFile::TF_AllocArray[] = { { "XYZ", 0, Traj_XYZ::WriteHelp, Traj_XYZ::Alloc }, { "LMOD conflib", 0, 0, Traj_Conflib::Alloc }, { "Gromacs dump", 0, Traj_GmxDump::WriteHelp, Traj_GmxDump::Alloc }, +# ifdef ENABLE_DTR + { "Desmond DTR", 0, 0, Traj_DTR::Alloc }, +# else + { "Desmond DTR", 0, 0, 0 }, +# endif { "Unknown trajectory", 0, 0, 0 } }; @@ -104,6 +110,7 @@ const FileTypes::KeyToken TrajectoryFile::TF_KeyArray[] = { { SQM, "sqm", ".sqm" }, { SDF, "sdf", ".sdf" }, { XYZ, "xyz", ".xyz" }, + { DTR, "dtr", ".dtr" }, { UNKNOWN_TRAJ, 0, 0 } }; diff --git a/src/TrajectoryFile.h b/src/TrajectoryFile.h index 5d5638e95b..c4fcd4ce65 100644 --- a/src/TrajectoryFile.h +++ b/src/TrajectoryFile.h @@ -28,7 +28,7 @@ class TrajectoryFile { CIF, CHARMMDCD, GMXTRX, GMXXTC, GMXTNG, BINPOS, AMBERRESTART, GRO, TINKER, CHARMMCOR, CHARMMREST, AMBERTRAJ, SQM, SDF, XYZ, - CONFLIB, GMXDUMP, + CONFLIB, GMXDUMP, DTR, UNKNOWN_TRAJ }; diff --git a/src/Version.h b/src/Version.h index 3c004fe453..0cccb140fc 100644 --- a/src/Version.h +++ b/src/Version.h @@ -12,7 +12,7 @@ * Whenever a number that precedes is incremented, all subsequent * numbers should be reset to 0. */ -#define CPPTRAJ_INTERNAL_VERSION "V4.25.3" +#define CPPTRAJ_INTERNAL_VERSION "V4.25.4" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/src/cpptrajdepend b/src/cpptrajdepend index d88419aa13..8d108ec087 100644 --- a/src/cpptrajdepend +++ b/src/cpptrajdepend @@ -287,6 +287,7 @@ FileIO_Mpi.o : FileIO_Mpi.cpp FileIO.h FileIO_Mpi.h Parallel.h FileIO_MpiShared.o : FileIO_MpiShared.cpp FileIO.h FileIO_Mpi.h FileIO_MpiShared.h Parallel.h FileIO_Std.o : FileIO_Std.cpp FileIO.h FileIO_Std.h FileName.o : FileName.cpp CpptrajStdio.h FileName.h StringRoutines.h +FileType.o : FileType.cpp FileType.h FileTypes.o : FileTypes.cpp ArgList.h BaseIOtype.h CpptrajStdio.h FileTypes.h File_TempName.o : File_TempName.cpp CpptrajStdio.h FileName.h File_TempName.h StringRoutines.h ForLoop.o : ForLoop.cpp AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h FileIO.h FileName.h ForLoop.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Vec3.h @@ -364,6 +365,7 @@ Traj_CharmmCor.o : Traj_CharmmCor.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h At Traj_CharmmDcd.o : Traj_CharmmDcd.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h ByteRoutines.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Traj_CharmmDcd.h TrajectoryIO.h TypeNameHolder.h Vec3.h Traj_CharmmRestart.o : Traj_CharmmRestart.cpp Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedFrame.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Topology.h Traj_CharmmRestart.h TrajectoryIO.h TypeNameHolder.h Vec3.h Traj_Conflib.o : Traj_Conflib.cpp Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Traj_Conflib.h TrajectoryIO.h TypeNameHolder.h Vec3.h +Traj_DTR.o : Traj_DTR.cpp Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Traj_DTR.h TrajectoryIO.h TypeNameHolder.h Vec3.h vmdplugin/dtrplugin.hxx Traj_GmxDump.o : Traj_GmxDump.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Traj_GmxDump.h TrajectoryIO.h TypeNameHolder.h Vec3.h Traj_GmxTng.o : Traj_GmxTng.cpp Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Traj_GmxTng.h TrajectoryIO.h TypeNameHolder.h Vec3.h Traj_GmxTrX.o : Traj_GmxTrX.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h ByteRoutines.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h Topology.h Traj_GmxTrX.h TrajectoryIO.h TypeNameHolder.h Vec3.h @@ -376,7 +378,7 @@ Traj_SDF.o : Traj_SDF.cpp Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Traj_SQM.o : Traj_SQM.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h Topology.h Traj_SQM.h TrajectoryIO.h TypeNameHolder.h Vec3.h Traj_Tinker.o : Traj_Tinker.cpp Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h SymbolExporting.h TinkerFile.h Topology.h Traj_Tinker.h TrajectoryIO.h TypeNameHolder.h Vec3.h Traj_XYZ.o : Traj_XYZ.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h BufferedLine.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Topology.h Traj_XYZ.h TrajectoryIO.h TypeNameHolder.h Vec3.h -TrajectoryFile.o : TrajectoryFile.cpp Atom.h AtomMask.h BaseIOtype.h Box.h BufferedFrame.h BufferedLine.h CIFfile.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Mol2File.h Molecule.h NameType.h NetcdfFile.h PDBfile.h Parallel.h ReplicaDimArray.h Residue.h SDFfile.h SymbolExporting.h TextFormat.h TinkerFile.h Traj_AmberCoord.h Traj_AmberNetcdf.h Traj_AmberRestart.h Traj_AmberRestartNC.h Traj_Binpos.h Traj_CIF.h Traj_CharmmCor.h Traj_CharmmDcd.h Traj_CharmmRestart.h Traj_Conflib.h Traj_GmxDump.h Traj_GmxTng.h Traj_GmxTrX.h Traj_GmxXtc.h Traj_Gro.h Traj_Mol2File.h Traj_NcEnsemble.h Traj_PDBfile.h Traj_SDF.h Traj_SQM.h Traj_Tinker.h Traj_XYZ.h TrajectoryFile.h TrajectoryIO.h Vec3.h +TrajectoryFile.o : TrajectoryFile.cpp Atom.h AtomMask.h BaseIOtype.h Box.h BufferedFrame.h BufferedLine.h CIFfile.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h MaskToken.h Matrix_3x3.h Mol2File.h Molecule.h NameType.h NetcdfFile.h PDBfile.h Parallel.h ReplicaDimArray.h Residue.h SDFfile.h SymbolExporting.h TextFormat.h TinkerFile.h Traj_AmberCoord.h Traj_AmberNetcdf.h Traj_AmberRestart.h Traj_AmberRestartNC.h Traj_Binpos.h Traj_CIF.h Traj_CharmmCor.h Traj_CharmmDcd.h Traj_CharmmRestart.h Traj_Conflib.h Traj_DTR.h Traj_GmxDump.h Traj_GmxTng.h Traj_GmxTrX.h Traj_GmxXtc.h Traj_Gro.h Traj_Mol2File.h Traj_NcEnsemble.h Traj_PDBfile.h Traj_SDF.h Traj_SQM.h Traj_Tinker.h Traj_XYZ.h TrajectoryFile.h TrajectoryIO.h Vec3.h TrajectoryIO.o : TrajectoryIO.cpp BaseIOtype.h Box.h CoordinateInfo.h FramePtrArray.h Matrix_3x3.h Parallel.h ReplicaDimArray.h TrajectoryIO.h Vec3.h TrajinList.o : TrajinList.cpp ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSet_RemLog.h Dimension.h EnsembleIn.h EnsembleIn_Multi.h EnsembleIn_Single.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h ReplicaInfo.h Residue.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h TrajIOarray.h TrajectoryFile.h TrajectoryIO.h Trajin.h TrajinList.h Trajin_Multi.h Trajin_Single.h TypeNameHolder.h Vec3.h Trajin_Multi.o : Trajin_Multi.cpp ArgList.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h CpptrajStdio.h FileName.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReplicaDimArray.h ReplicaInfo.h Residue.h StringRoutines.h SymbolExporting.h Topology.h TrajFrameCounter.h TrajIOarray.h TrajectoryIO.h Trajin.h Trajin_Multi.h TypeNameHolder.h Vec3.h @@ -387,3 +389,4 @@ Vec3.o : Vec3.cpp Constants.h CpptrajStdio.h Vec3.h ViewRst.o : ViewRst.cpp ActionFrameCounter.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h BondSearch.h Box.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h FileIO.h FileName.h FileTypes.h Frame.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h OutputTrajCommon.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h ParmFile.h ParmIO.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajectoryFile.h Trajout_Single.h TypeNameHolder.h Vec3.h ViewRst.h main.o : main.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Constants.h CoordinateInfo.h Cpptraj.h CpptrajFile.h CpptrajState.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Vec3.h molsurf.o : molsurf.c molsurf.h +vmdplugin/dtrplugin.o : vmdplugin/dtrplugin.cpp vmdplugin/../ByteRoutines.h vmdplugin/dtrplugin.hxx vmdplugin/vmddir.h diff --git a/src/cpptrajfiles b/src/cpptrajfiles index ec7fc7982f..e34ddcd3d8 100644 --- a/src/cpptrajfiles +++ b/src/cpptrajfiles @@ -288,6 +288,7 @@ COMMON_SOURCES= \ FileIO_MpiShared.cpp \ FileIO_Std.cpp \ FileName.cpp \ + FileType.cpp \ FileTypes.cpp \ ForLoop.cpp \ ForLoop_dataSetBlocks.cpp \ @@ -360,6 +361,7 @@ COMMON_SOURCES= \ Traj_CharmmRestart.cpp \ Traj_CIF.cpp \ Traj_Conflib.cpp \ + Traj_DTR.cpp \ Traj_GmxDump.cpp \ Traj_GmxTng.cpp \ Traj_GmxTrX.cpp \ @@ -382,7 +384,8 @@ COMMON_SOURCES= \ Trajout_Single.cpp \ TrajoutList.cpp \ Vec3.cpp \ - ViewRst.cpp + ViewRst.cpp \ + vmdplugin/dtrplugin.cpp CSOURCES= molsurf.c diff --git a/src/vmdplugin/README b/src/vmdplugin/README new file mode 100644 index 0000000000..053c3d5c60 --- /dev/null +++ b/src/vmdplugin/README @@ -0,0 +1 @@ +These are the DTR and DMS plugins (and some of the associated framework) from VMD version 1.9.4a diff --git a/src/vmdplugin/dmsplugin.cpp b/src/vmdplugin/dmsplugin.cpp new file mode 100644 index 0000000000..edc1b83340 --- /dev/null +++ b/src/vmdplugin/dmsplugin.cpp @@ -0,0 +1,919 @@ +/* +Copyright 2009, D. E. Shaw Research, LLC +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +* Redistributions of source code must retain the above copyright + notice, this list of conditions, and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright + notice, this list of conditions, and the following disclaimer in the + documentation and/or other materials provided with the distribution. + +* Neither the name of D. E. Shaw Research, LLC nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ + +#include "molfile_plugin.h" +#include + +#include /* for isspace() */ +#include +#include +#include +#include +#include +#include + + +#ifdef _MSC_VER +#include +static int fs_remove_file( const char * path ) { + DeleteFileA( path ); + return 1; /* FIXME: check for failure */ +} +#else +#include +#include +static int fs_remove_file( const char * path ) { + if (unlink(path)!=0 && errno!=ENOENT) { + fprintf(stderr, "Removing %s failed: %s\n", path, strerror(errno)); + return 0; + } + return 1; +} +#endif + +#ifndef M_PI +#define M_PI (3.1415926535897932385) +#endif + +#ifndef M_PI_2 +#define M_PI_2 (1.5707963267948966192) +#endif + +struct element { + double daltons; + const char* abbreviation; + const char* name; +}; + +/* url = "http://physics.nist.gov/cgi-bin/Elements/elInfo.pl?element=%d&context=noframes"%element */ +static struct element amu[] = { + {1.00794,"H","Hydrogen"}, + {4.002602,"He","Helium"}, + {6.941,"Li","Lithium"}, + {9.012182,"Be","Beryllium"}, + {10.811,"B","Boron"}, + {12.0107,"C","Carbon"}, + {14.0067,"N","Nitrogen"}, + {15.9994,"O","Oxygen"}, + {18.9984032,"F","Fluorine"}, + {20.1797,"Ne","Neon"}, + {22.989770,"Na","Sodium"}, + {24.3050,"Mg","Magnesium"}, + {26.981538,"Al","Aluminum"}, + {28.0855,"Si","Silicon"}, + {30.973761,"P","Phosphorus"}, + {32.065,"S","Sulfur"}, + {35.453,"Cl","Chlorine"}, + {39.0983,"K","Potassium"}, + {39.948,"Ar","Argon"}, + {40.078,"Ca","Calcium"}, + {44.955910,"Sc","Scandium"}, + {47.867,"Ti","Titanium"}, + {50.9415,"V","Vanadium"}, + {51.9961,"Cr","Chromium"}, + {54.938049,"Mn","Manganese"}, + {55.845,"Fe","Iron"}, + {58.6934,"Ni","Nickel"}, + {58.933200,"Co","Cobalt"}, + {63.546,"Cu","Copper"}, + {65.409,"Zn","Zinc"}, + {69.723,"Ga","Gallium"}, + {72.64,"Ge","Germanium"}, + {74.92160,"As","Arsenic"}, + {78.96,"Se","Selenium"}, + {79.904,"Br","Bromine"}, + {83.798,"Kr","Krypton"}, + {85.4678,"Rb","Rubidium"}, + {87.62,"Sr","Strontium"}, + {88.90585,"Y","Yttrium"}, + {91.224,"Zr","Zirconium"}, + {92.90638,"Nb","Niobium"}, + {95.94,"Mo","Molybdenum"}, + {101.07,"Ru","Ruthenium"}, + {102.90550,"Rh","Rhodium"}, + {106.42,"Pd","Palladium"}, + {107.8682,"Ag","Silver"}, + {112.411,"Cd","Cadmium"}, + {114.818,"In","Indium"}, + {118.710,"Sn","Tin"}, + {121.760,"Sb","Antimony"}, + {126.90447,"I","Iodine"}, + {127.60,"Te","Tellurium"}, + {131.293,"Xe","Xenon"}, + {132.90545,"Cs","Cesium"}, + {137.327,"Ba","Barium"}, + {138.9055,"La","Lanthanum"}, + {140.116,"Ce","Cerium"}, + {140.90765,"Pr","Praseodymium"}, + {144.24,"Nd","Neodymium"}, + {150.36,"Sm","Samarium"}, + {151.964,"Eu","Europium"}, + {157.25,"Gd","Gadolinium"}, + {158.92534,"Tb","Terbium"}, + {162.500,"Dy","Dysprosium"}, + {164.93032,"Ho","Holmium"}, + {167.259,"Er","Erbium"}, + {168.93421,"Tm","Thulium"}, + {173.04,"Yb","Ytterbium"}, + {174.967,"Lu","Lutetium"}, + {178.49,"Hf","Hafnium"}, + {180.9479,"Ta","Tantalum"}, + {183.84,"W","Tungsten"}, + {186.207,"Re","Rhenium"}, + {190.23,"Os","Osmium"}, + {192.217,"Ir","Iridium"}, + {195.078,"Pt","Platinum"}, + {196.96655,"Au","Gold"}, + {200.59,"Hg","Mercury"}, + {204.3833,"Tl","Thallium"}, + {207.2,"Pb","Lead"}, + {208.98038,"Bi","Bismuth"}, + {231.03588,"Pa","Protactinium"}, + {232.0381,"Th","Thorium"}, + {238.02891,"U","Uranium"} +}; + +static const int nelements = sizeof(amu)/sizeof(amu[0]); + +static const char * +find_element_by_atomic_number(int target) { + if (target < 1) target=1; + if (target >= nelements) target = nelements; + return amu[target-1].abbreviation; +} + +static int +find_element_by_amu(double target) { + int left = 0; + int right = nelements-1; + int swap; + + /* ----------------------------------------------- + // Knuth's binary search + // ----------------------------------------------- */ + while(left <= right) { + int mid = (left+right)/2; + if (target> amu[mid].daltons) { + left = mid + 1; + } else if (target< amu[mid].daltons) { + right = mid - 1; + } else { + /* Exact match (unlikely) */ + left = right = mid; + return left+1; + } + } + + /* ----------------------------------------------- + // CAUTION: at this point, the meanings of + // left and right are switched (i.e. left >= right, + // see the while() loop above if you don't believe me! + // ----------------------------------------------- */ + swap = left; + left = right; + right = swap; + + if (left < 0) left = right; + if (right > nelements-1) right = left; + + if (target - amu[left].daltons < amu[right].daltons - target) { + return left+1; + } + + return right+1; +} + +static int table_size( sqlite3 * db, const char * tname, + sqlite_int64 * count ) { + sqlite3_stmt *stmt; + char * buf; + + if (!tname) return 0; + buf = (char *)malloc(strlen(tname) + 100); + sprintf(buf, "select count() from '%s'", tname); + if (sqlite3_prepare_v2(db, buf, -1, &stmt, NULL)) { + free(buf); + return 0; + } + if (SQLITE_ROW != sqlite3_step(stmt)) { + sqlite3_finalize(stmt); + free(buf); + return 0; + } + if (count) *count = sqlite3_column_int64(stmt, 0); + sqlite3_finalize(stmt); + return 1; /* success */ +} + +static molfile_plugin_t plugin; +//static molfile_plugin_t append_plugin; + +namespace { + struct Handle { + sqlite3 * db; + int natoms; + int frames_read; + std::vector from, to, gids; + std::vector glue_from, glue_to; + std::vector order; + Handle(sqlite3 * _db, int n) + : db(_db), natoms(n), frames_read(0) {} + + ~Handle() { + if (db && sqlite3_close(db)) + fprintf(stderr, "Error closing db: %s\n", sqlite3_errmsg(db)); + } + }; +} + + + static void *open_file_read( const char *filename, const char *filetype, + int *natoms ) { + + sqlite3 * db; + sqlite_int64 count; + Handle * h; + if (sqlite3_open_v2( filename, &db, SQLITE_OPEN_READONLY, NULL)) { + fprintf(stderr, "Error opening dms at %s: %s", + filename, sqlite3_errmsg(db)); + return NULL; + } + if (!table_size( db, "particle", &count )) return NULL; + h = new Handle(db, count); + *natoms = count; + return h; + } + +static void close_file_read(void *v) { delete (Handle *)v; } + +static int has_nonwhitespace( const char * p ) { + for (; *p; ++p) { + if (!isspace(*p)) return 1; + } + return 0; +} + +#define GET_STRING( col, field ) do { \ + if ( col >= 0 ) { \ + const char * val = (const char *)sqlite3_column_text(stmt, col); \ + if (val) { \ + strncpy( atom->field, val, sizeof(atom->field) ); \ + atom->field[sizeof(atom->field)-1] = '\0'; \ + } \ + } \ +} while(0) + +#define GET_INT( col, field ) do { \ + if ( col >= 0 ) { \ + atom->field = sqlite3_column_int(stmt, col); \ + } \ +} while (0) + +#define GET_DOUBLE( col, field ) do { \ + if ( col >= 0 ) { \ + atom->field = sqlite3_column_double(stmt, col); \ + } \ +} while (0) + + static + int read_structure( void *v, int *optflags, molfile_atom_t *atoms ) { + Handle *h = (Handle *)v; + sqlite3 * db = h->db; + sqlite3_stmt * stmt; + static const char infosql[] = "pragma table_info(particle)"; + static const char loadsql[] = "select * from particle"; + molfile_atom_t * atom = NULL; + + int gid_col = -1; + int name_col = -1; + int anum_col = -1; + int resid_col = -1; + int resname_col = -1; + int chain_col = -1; + int segid_col = -1; + int mass_col = -1; + int charge_col = -1; + int x_col = -1; + int y_col = -1; + int z_col = -1; + int vx_col = -1; + int vy_col = -1; + int vz_col = -1; + int occ_col = -1; + + int ncols=0; + *optflags = 0; + + if (sqlite3_prepare_v2(db, infosql, -1, &stmt, NULL)) { + fprintf(stderr, "Error getting particle info: %s", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + atom = atoms; + while (SQLITE_ROW==sqlite3_step(stmt)) { + const char * name = (const char *)sqlite3_column_text(stmt,1); + + if (!strcmp(name, "anum")) anum_col = ncols; + else if (!strcmp(name, "id")) gid_col = ncols; + else if (!strcmp(name, "name")) name_col = ncols; + else if (!strcmp(name, "resid")) resid_col = ncols; + else if (!strcmp(name, "resname")) resname_col = ncols; + else if (!strcmp(name, "chain")) chain_col = ncols; + else if (!strcmp(name, "segid")) segid_col = ncols; + else if (!strcmp(name, "mass")) mass_col = ncols; + else if (!strcmp(name, "charge")) charge_col = ncols; + else if (!strcmp(name, "x")) x_col = ncols; + else if (!strcmp(name, "y")) y_col = ncols; + else if (!strcmp(name, "z")) z_col = ncols; + else if (!strcmp(name, "vx")) vx_col = ncols; + else if (!strcmp(name, "vy")) vy_col = ncols; + else if (!strcmp(name, "vz")) vz_col = ncols; + else if (!strcmp(name, "occupancy")) occ_col = ncols; + ++ncols; + } + sqlite3_finalize(stmt); + + if (anum_col >= 0) *optflags |= MOLFILE_ATOMICNUMBER; + if (charge_col >= 0) *optflags |= MOLFILE_CHARGE; + if (mass_col >= 0) *optflags |= MOLFILE_MASS; + if (occ_col >= 0) *optflags |= MOLFILE_OCCUPANCY; + + if (gid_col<0) { + fprintf(stderr, "Missing id column in particle table\n"); + return MOLFILE_ERROR; + } + + if (sqlite3_prepare_v2(db, loadsql, -1, &stmt, NULL)) { + fprintf(stderr, "Error loading particle structure: %s", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + while (SQLITE_ROW==sqlite3_step(stmt)) { + memset(atom, 0, sizeof(*atom)); + GET_STRING( name_col, name ); + GET_INT( anum_col, atomicnumber ); + GET_INT( resid_col, resid ); + GET_STRING( resname_col, resname ); + GET_DOUBLE( mass_col, mass ); + GET_DOUBLE( charge_col, charge ); + GET_DOUBLE( occ_col, occupancy ); + GET_STRING( chain_col, chain ); + GET_STRING( segid_col, segid ); + + h->gids.push_back(sqlite3_column_int(stmt, gid_col)); + + /* if the name is all space, and we have an atomic number, use the + element name instead. */ + if (atom->atomicnumber) { + if (!has_nonwhitespace(atom->name)) { + const char *nm = + find_element_by_atomic_number(atom->atomicnumber); + strncpy( atom->name, nm, sizeof(atom->name)); + } + } + ++atom; + } + sqlite3_finalize(stmt); + return MOLFILE_SUCCESS; + } + +static +#if vmdplugin_ABIVERSION > 14 +int read_bonds(void *v, int *nbonds, int **from, int **to, float** order, + int **bondtype, int *nbondtypes, char ***bondtypename) { +#else +int read_bonds(void *v, int *nbonds, int **from, int **to, float** order) { +#endif + Handle *h = (Handle *)v; + sqlite3 * db = h->db; + sqlite_int64 count = 0; + if (table_size(h->db, "bond", &count) && count>0) { + sqlite3_stmt * stmt; + int with_order = 1; + /* Try including order first */ + static const char * sql = "select p0,p1,\"order\" from bond"; + if (sqlite3_prepare_v2(db, sql, -1, &stmt, NULL)) { + /* Fall back to just p0, p1 */ + static const char * sql = "select p0,p1 from bond"; + with_order = 0; + if (sqlite3_prepare_v2(db, sql, -1, &stmt, NULL)) { + fprintf(stderr, "Error reading bonds: %s", sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + } + h->from.resize(count); + h->to.resize(count); + h->order.resize(count); + + *from = &h->from[0]; + *to = &h->to[0]; + *order = &h->order[0]; +#if vmdplugin_ABIVERSION > 14 + *bondtype = NULL; + *nbondtypes = 0; + *bondtypename = NULL; + #endif + + std::map gidmap; + for (unsigned i=0; igids.size(); i++) gidmap[h->gids[i]]=i+1; + + for (int i=0; sqlite3_step(stmt)==SQLITE_ROW; i++) { + std::map::const_iterator iter; + + int from = sqlite3_column_int(stmt,0); + int to = sqlite3_column_int(stmt,1); + h->order[i] = with_order ? sqlite3_column_int(stmt,2) : 1; + if ((iter=gidmap.find(from))==gidmap.end()) { + fprintf(stderr, "Illegal bond in row %d: %d-%d\n", + i+1, from, to); + return MOLFILE_ERROR; + } + h->from[i] = iter->second; + if ((iter=gidmap.find(to))==gidmap.end()) { + fprintf(stderr, "Illegal bond in row %d: %d-%d\n", + i+1, from, to); + return MOLFILE_ERROR; + } + h->to[i] = iter->second; + } + sqlite3_finalize(stmt); + } + *nbonds = count; + return MOLFILE_SUCCESS; +} + +static +int read_fictitious_bonds(void *v, int *nbonds, int **from, int **to) { + Handle *h = (Handle *)v; + sqlite3 * db = h->db; + sqlite_int64 count = 0; + if (table_size(h->db, "glue", &count) && count>0) { + sqlite3_stmt * stmt; + static const char * sql = "select p0,p1 from glue"; + if (sqlite3_prepare_v2(db, sql, -1, &stmt, NULL)) { + fprintf(stderr, "Error reading bonds: %s", sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + h->glue_from.resize(count); + h->glue_to.resize(count); + + *from = &h->glue_from[0]; + *to = &h->glue_to[0]; + + std::map gidmap; + for (unsigned i=0; igids.size(); i++) gidmap[h->gids[i]]=i+1; + + for (int i=0; sqlite3_step(stmt)==SQLITE_ROW; i++) { + std::map::const_iterator iter; + + int from = sqlite3_column_int(stmt,0); + int to = sqlite3_column_int(stmt,1); + if ((iter=gidmap.find(from))==gidmap.end()) { + fprintf(stderr, "Illegal glue in row %d: %d-%d\n", + i+1, from, to); + return MOLFILE_ERROR; + } + h->glue_from[i] = iter->second; + if ((iter=gidmap.find(to))==gidmap.end()) { + fprintf(stderr, "Illegal glue in row %d: %d-%d\n", + i+1, from, to); + return MOLFILE_ERROR; + } + h->glue_to[i] = iter->second; + } + sqlite3_finalize(stmt); + } + *nbonds = count; + return MOLFILE_SUCCESS; +} + + static + double dotprod(const double *x, const double *y) { + return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; + } + + static + int read_global_cell(Handle * h, molfile_timestep_t * ts) { + sqlite3_stmt * stmt; + sqlite_int64 count; + int i,j; + double cell[3][3]; + double cosAB, cosAC, cosBC; + const double * A, * B, * C; + if (!table_size(h->db, "global_cell", &count)) return MOLFILE_SUCCESS; + if (count != 3) { + fprintf(stderr, "Error: expected global_cell size of 3, got %d", + (int)count); + return MOLFILE_ERROR; + } + if (sqlite3_prepare_v2(h->db,"select x,y,z from global_cell", -1, &stmt, NULL)) { + fprintf(stderr, "Error compiling global_cell reader: %s", + sqlite3_errmsg(h->db)); + return MOLFILE_ERROR; + } + for (i=0; i<3; i++) { + sqlite3_step(stmt); + for (j=0; j<3; j++) { + cell[i][j] = sqlite3_column_double(stmt,j); + } + } + sqlite3_finalize(stmt); + + A = cell[0]; + B = cell[1]; + C = cell[2]; + + /* store lengths */ + ts->A = sqrt(dotprod(A,A)); + ts->B = sqrt(dotprod(B,B)); + ts->C = sqrt(dotprod(C,C)); + + if (ts->A==0 || ts->B==0 || ts->C==0) { + ts->alpha = ts->beta = ts->gamma = 90.0; + + } else { + /* compute angles */ + cosAB = dotprod(A,B)/(ts->A * ts->B); + cosAC = dotprod(A,C)/(ts->A * ts->C); + cosBC = dotprod(B,C)/(ts->B * ts->C); + + /* clamp */ + if (cosAB > 1.0) cosAB = 1.0; else if (cosAB < -1.0) cosAB = -1.0; + if (cosAC > 1.0) cosAC = 1.0; else if (cosAC < -1.0) cosAC = -1.0; + if (cosBC > 1.0) cosBC = 1.0; else if (cosBC < -1.0) cosBC = -1.0; + + /* convert to angles using asin to avoid nasty rounding when we are + close to 90 degree angles. */ + ts->alpha = 90.0 - asin(cosBC) * 90.0 / M_PI_2; /* cosBC */ + ts->beta = 90.0 - asin(cosAC) * 90.0 / M_PI_2; /* cosAC */ + ts->gamma = 90.0 - asin(cosAB) * 90.0 / M_PI_2; /* cosAB */ + } + return MOLFILE_SUCCESS; + } + + static + int read_timestep(Handle * h, molfile_timestep_t *ts) { + + sqlite3_stmt * stmt; + float * pos = ts->coords; + float * vel = ts->velocities; + const char * sql = vel ? "select x,y,z,vx,vy,vz from particle" + : "select x,y,z from particle"; + if (sqlite3_prepare_v2( h->db, sql, -1, &stmt, NULL)) { + fprintf(stderr, "Error reading timestep: %s\n", + sqlite3_errmsg(h->db)); + return MOLFILE_ERROR; + } + while (SQLITE_ROW==sqlite3_step(stmt)) { + int i; + for (i=0; i<3; i++) *pos++ = sqlite3_column_double(stmt,i); + if (vel) { + for (i=0; i<3; i++) *vel++ = sqlite3_column_double(stmt,i+3); + } + } + sqlite3_finalize(stmt); + return read_global_cell(h,ts); + } + + + static + int read_next_timestep( void *v, int atoms, molfile_timestep_t *ts) { + Handle *h = (Handle *)v; + if (h->frames_read++) return MOLFILE_EOF; + return read_timestep(h, ts); + } + + static + int read_timestep2(void *v, ssize_t n, molfile_timestep_t *ts) { + Handle *h = (Handle *)v; + if (n!=0) return MOLFILE_EOF; + return read_timestep(h, ts); + } + + static + int read_timestep_metadata( void *v, molfile_timestep_metadata_t *m) { + /* FIXME: assume velocities */ + m->has_velocities = 1; + m->count = 1; /* number of timesteps */ + m->avg_bytes_per_timestep = 10000; /* FIXME */ + return MOLFILE_SUCCESS; + } + + static + void *open_file_write(const char *path, const char *type, int natoms) { + + sqlite3 * db; + Handle *h=NULL; + /* delete existing file so we don't just append */ + if (!fs_remove_file(path)) return NULL; + if (sqlite3_open(path, &db)) { + fprintf(stderr, "Failed opening %s: %s\n", + path, sqlite3_errmsg(db)); + return NULL; + } + h = new Handle(db, natoms); + return h; + } + + static +#if vmdplugin_ABIVERSION > 14 + int write_bonds(void *v, int nbonds, int *from, int *to, float *order, + int *bondtype, int nbondtypes, char **bondtypename) { +#else + int write_bonds(void *v, int nbonds, int *from, int *to, float *order) { +#endif + Handle *h = (Handle *)v; + h->from.resize(nbonds); + h->to.resize(nbonds); + + std::copy(from, from+nbonds, h->from.begin()); + std::copy(to, to+nbonds, h->to.begin()); + if (order) { + h->order.resize(nbonds); + std::copy(order, order+nbonds, h->order.begin()); + } else { + h->order.clear(); + h->order.insert(h->order.begin(), nbonds, 1); + } + return MOLFILE_SUCCESS; + } + + static + int write_structure(void *v, int optflags, const molfile_atom_t *atoms) { + Handle *h = (Handle *)v; + sqlite3 * db = h->db; + sqlite3_stmt * stmt; + int i; + + sqlite3_exec(db, "begin", NULL, NULL, NULL); + if (sqlite3_exec(db, + " create table if not exists particle (\n" + "id integer primary key, \n" + "anum integer,\n" + "x float, y float, z float,\n" + "vx float, vy float, vz float,\n" + "mass float, charge float,\n" + "name text, resname text, resid integer, chain text, segid text)", + NULL, NULL, NULL)) { + fprintf(stderr, "Error creating particle table %s\n", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + if (sqlite3_prepare_v2(db, "insert into particle values (" + "?," /* id */ + "?," /* anum */ + "?,?,?,?,?,?," /* x,y,z, vx,vy,vz */ + "?,?," /* mass, charge */ + "?,?,?,?,?)", /* name, resname, resid, chain, segid */ + -1, &stmt, NULL)) { + fprintf(stderr, "Error compiling insert particle statement %s\n", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + for (i=0; inatoms; i++) { + const molfile_atom_t * a = atoms+i; + int anum=0; + if (optflags & MOLFILE_ATOMICNUMBER) { + anum=a->atomicnumber; + } else if (optflags & MOLFILE_MASS) { + anum=find_element_by_amu(a->mass); + } + + sqlite3_bind_int(stmt, 1, i ); + sqlite3_bind_int(stmt, 2, anum ); + sqlite3_bind_double( stmt, 9, + optflags & MOLFILE_MASS ? a->mass : 0.0 ); + sqlite3_bind_double( stmt, 10, + optflags & MOLFILE_CHARGE ? a->charge : 0.0 ); + sqlite3_bind_text( stmt, 11, a->name, -1, SQLITE_STATIC ); + sqlite3_bind_text( stmt, 12, a->resname, -1, SQLITE_STATIC ); + sqlite3_bind_int( stmt, 13, a->resid ); + sqlite3_bind_text( stmt, 14, a->chain, -1, SQLITE_STATIC ); + sqlite3_bind_text( stmt, 15, a->segid, -1, SQLITE_STATIC ); + + if (sqlite3_step(stmt)!=SQLITE_DONE) { + fprintf(stderr, "Error adding particle: %s\n", + sqlite3_errmsg(db)); + sqlite3_finalize(stmt); + return MOLFILE_ERROR; + } + sqlite3_reset(stmt); + } + sqlite3_finalize(stmt); + + /* bonds */ + if (sqlite3_exec(db, + " create table if not exists bond (\n" + " p0 integer, p1 integer, 'order' integer)", NULL, NULL, NULL )) { + fprintf(stderr, "Error creating bond table: %s\n", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + if (sqlite3_prepare_v2(db, "insert into bond values (?,?,?)", -1, &stmt, NULL)) { + fprintf(stderr, "Error compiling bond insert statement: %s\n", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + for (unsigned i=0; ifrom.size(); i++) { + sqlite3_bind_int( stmt, 1, h->from[i]-1); + sqlite3_bind_int( stmt, 2, h->to[i]-1); + sqlite3_bind_int( stmt, 3, h->order[i]); + if (sqlite3_step(stmt)!=SQLITE_DONE) { + fprintf(stderr, "Error adding bond: %s\n", + sqlite3_errmsg(db)); + sqlite3_finalize(stmt); + return MOLFILE_ERROR; + } + sqlite3_reset(stmt); + } + sqlite3_finalize(stmt); + sqlite3_exec(db, "commit", NULL, NULL, NULL); + return MOLFILE_SUCCESS; + } + + static + void convert_to_homebox( const molfile_timestep_t * ts, double box[9]) { + + double A[3], B[3], C[3]; + + /* Convert VMD's unit cell information */ + double cosBC = sin( ((90 - ts->alpha ) / 180) * M_PI ); + double cosAC = sin( ((90 - ts->beta ) / 180) * M_PI ); + double cosAB = sin( ((90 - ts->gamma ) / 180) * M_PI ); + double sinAB = cos( ((90 - ts->gamma ) / 180) * M_PI ); + + double Ax = ts->A; + double Ay = 0; + double Az = 0; + double Bx = ts->B * cosAB; + double By = ts->B * sinAB; + double Bz = 0; + double Cx,Cy,Cz; + if (sinAB != 0) { + Cx = cosAC; + Cy = (cosBC - cosAC*cosAB) / sinAB; + Cz = sqrt(1-Cx*Cx-Cy*Cy); + Cx *= ts->C; + Cy *= ts->C; + Cz *= ts->C; + } else { + Cx=Cy=Cz=0; + } + A[0] = Ax; A[1] = Ay; A[2] = Az; + B[0] = Bx; B[1] = By; B[2] = Bz; + C[0] = Cx; C[1] = Cy; C[2] = Cz; + + /* put vectors in rows of box */ + box[0] = A[0]; box[1] = A[1]; box[2] = A[2]; + box[3] = B[0]; box[4] = B[1]; box[5] = B[2]; + box[6] = C[0]; box[7] = C[1]; box[8] = C[2]; + } + + static + int write_timestep(void *v, const molfile_timestep_t *ts) { + Handle *h = (Handle *)v; + sqlite3 * db = h->db; + sqlite3_stmt * stmt; + double box[9]; + int i; + if (h->frames_read) return MOLFILE_EOF; + static const char sql[]= + "update particle set x=?,y=?,z=?,vx=?,vy=?,vz=? where id=?"; + if (sqlite3_prepare_v2(db, sql, -1, &stmt, NULL)) { + fprintf(stderr, "Error compiling update positions statement: %s\n", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + sqlite3_exec(db, "begin", NULL, NULL, NULL); + for (i=0; inatoms; i++) { + sqlite3_bind_double( stmt, 1, ts->coords[3*i ] ); + sqlite3_bind_double( stmt, 2, ts->coords[3*i+1] ); + sqlite3_bind_double( stmt, 3, ts->coords[3*i+2] ); + double vx=0, vy=0, vz=0; + if (ts->velocities) { + vx=ts->velocities[3*i]; + vy=ts->velocities[3*i+1]; + vz=ts->velocities[3*i+2]; + } + sqlite3_bind_double( stmt, 4, vx ); + sqlite3_bind_double( stmt, 5, vy ); + sqlite3_bind_double( stmt, 6, vz ); + sqlite3_bind_int( stmt, 7, i ); + if (sqlite3_step(stmt)!=SQLITE_DONE) { + fprintf(stderr, "Error updating position: %s\n", + sqlite3_errmsg(db)); + sqlite3_finalize(stmt); + return MOLFILE_ERROR; + } + sqlite3_reset(stmt); + } + sqlite3_finalize(stmt); + + if (sqlite3_exec(db, + " create table if not exists global_cell (\n" + " id integer primary key, x float, y float, z float)", + NULL, NULL, NULL)) { + fprintf(stderr, "Error creating global cell table %s\n", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + sqlite3_exec(db, "delete from global_cell", NULL, NULL, NULL); + if (sqlite3_prepare_v2(db, "insert into global_cell (x,y,z) values (?,?,?)", + -1, &stmt, NULL)) { + fprintf(stderr, "Error compiling global_cell statement: %s\n", + sqlite3_errmsg(db)); + return MOLFILE_ERROR; + } + convert_to_homebox( ts, box ); + for (i=0; i<3; i++) { + sqlite3_bind_double(stmt, 1, box[3*i ]); + sqlite3_bind_double(stmt, 2, box[3*i+1]); + sqlite3_bind_double(stmt, 3, box[3*i+2]); + if (sqlite3_step(stmt)!=SQLITE_DONE) { + fprintf(stderr, "Error updating global cell: %s\n", + sqlite3_errmsg(db)); + sqlite3_finalize(stmt); + return MOLFILE_ERROR; + } + sqlite3_reset(stmt); + } + sqlite3_finalize(stmt); + sqlite3_exec(db, "commit", NULL, NULL, NULL); + return MOLFILE_SUCCESS; + } + + static void close_file_write(void *v) { delete (Handle *)v; } + + +VMDPLUGIN_API int VMDPLUGIN_init(void) { + memset(&plugin, 0, sizeof(plugin)); + plugin.abiversion = vmdplugin_ABIVERSION; + plugin.type = MOLFILE_PLUGIN_TYPE; + plugin.name = "dms"; + plugin.prettyname = "DESRES Molecular Structure"; + plugin.author = "D.E. Shaw Research, LLC: Justin Gullingsrud"; + plugin.majorv = 1; + plugin.minorv = 0; + plugin.is_reentrant = VMDPLUGIN_THREADSAFE; + + plugin.filename_extension = "dms"; + plugin.open_file_read = open_file_read; + plugin.read_structure = read_structure; + plugin.read_bonds = read_bonds; + plugin.read_next_timestep = read_next_timestep; +#if defined(DESRES_READ_TIMESTEP2) + plugin.read_timestep2 = read_timestep2; + plugin.read_fictitious_bonds = read_fictitious_bonds; +#endif + plugin.close_file_read = close_file_read; + + plugin.read_timestep_metadata = read_timestep_metadata; + + plugin.open_file_write = open_file_write; + plugin.write_structure = write_structure; + plugin.write_bonds = write_bonds; + plugin.write_timestep = write_timestep; + plugin.close_file_write = close_file_write; + + return VMDPLUGIN_SUCCESS; +} + +VMDPLUGIN_API int VMDPLUGIN_register( void *v, vmdplugin_register_cb cb) { + cb( v, (vmdplugin_t *)&plugin); + return VMDPLUGIN_SUCCESS; +} + +VMDPLUGIN_API int VMDPLUGIN_fini(void) { + return VMDPLUGIN_SUCCESS; +} + diff --git a/src/vmdplugin/dtrplugin.cpp b/src/vmdplugin/dtrplugin.cpp new file mode 100644 index 0000000000..9fba377ac2 --- /dev/null +++ b/src/vmdplugin/dtrplugin.cpp @@ -0,0 +1,2381 @@ +// +// Version info for VMD plugin tree: +// $Id: dtrplugin.cxx,v 1.26 2016/11/06 19:25:35 johns Exp $ +// +// Version info for last sync with D. E. Shaw Research: +// //depot/desrad/main/sw/libs/molfile/plugins/dtrplugin.cxx#30 +// + +/* +Copyright 2009, D. E. Shaw Research, LLC +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +* Redistributions of source code must retain the above copyright + notice, this list of conditions, and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright + notice, this list of conditions, and the following disclaimer in the + documentation and/or other materials provided with the distribution. + +* Neither the name of D. E. Shaw Research, LLC nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ +#ifdef ENABLE_DTR +#include "dtrplugin.hxx" + +using namespace desres::molfile; + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "vmddir.h" +//#include "endianswap.h" +#include "../ByteRoutines.h" + +static const char SERIALIZED_VERSION[] = "0006"; +const char * desres::molfile::dtr_serialized_version() { + return SERIALIZED_VERSION; +} + +static bool badversion(const std::string& version) { + return version != SERIALIZED_VERSION; +} + +#ifndef _WIN32 +static const char s_sep = '/'; +#include + +#include /* for htonl */ +#if defined(_AIX) || defined(ANDROID) +#include +#else +#include +#endif + +#if defined(__sun) +#include +#include +#include +#endif + +#else +/// windows version + +#define M_PI (3.1415926535897932385) +#define M_PI_2 (1.5707963267948966192) + +#ifndef S_ISREG +#define S_ISREG(x) (((x) & S_IFMT) == S_IFREG) +#endif + +#ifndef S_ISDIR +#define S_ISDIR(x) (((x) & S_IFMT) == S_IFDIR) +#endif + +#ifndef PRIu64 +#if _MSC_VER < 1400 +#define PRIu64 "ld" +#else +#define PRIu64 "I64u" +#endif +#endif + +static const char s_sep = '\\'; + +#endif + +static const uint32_t magic_timekey = 0x4445534b; +static const uint32_t magic_frame = 0x4445534d; +static const uint32_t s_version = 0x00000100; +static const uint32_t s_irosetta = 0x12345678; +static const float s_frosetta = 1234.5; +static const double s_drosetta = 1234.5e6; +static const uint32_t s_lrosetta_lo = 0x89abcdef; +static const uint32_t s_lrosetta_hi = 0x01234567; +static const uint32_t s_blocksize = 4096; +static const uint32_t s_alignsize = 8; + +namespace { + + const double PEAKmassInAmu = 418.4; + + double sfxp_ulp32flt(int32_t x) { + return ldexp(((double) x),-31); + } + + + uint64_t assemble64( uint32_t lo, uint32_t hi) { + uint64_t hi64 = hi; + return (hi64 << 32) | lo; + } + + double assembleDouble(uint32_t lo, uint32_t hi) { + union { + uint64_t ival; + double dval; + } u; + u.ival = assemble64(lo,hi); + return u.dval; + } + + /// definitions of binary representation of frameset files + + typedef struct { + uint32_t magic; //!< Magic number + uint32_t version; //!< Version of creator + uint32_t framesize_lo; //!< bytes in frame (low) + uint32_t framesize_hi; //!< bytes in frame (high) + } required_header_t; + + //! Header structure within file. + typedef struct { + required_header_t required; //!< 4 word mini-header + + uint32_t size_header_block; //!< Size of this header + uint32_t unused0; //!< not used in current implementation + uint32_t irosetta; //!< 32-bit integer rosetta value. + float frosetta; //!< 32-bit float rosetta + + uint32_t drosetta_lo; //!< 64-bit float rosetta (low) + uint32_t drosetta_hi; //!< 64-bit float rosetta (high) + uint32_t lrosetta_lo; //!< 64-bit integer rosetta (low) + uint32_t lrosetta_hi; //!< 64-bit integer rosetta (high) + + uint32_t endianism; //!< Endianism of writer machine. + uint32_t nlabels; //!< Number of labeled fields. + uint32_t size_meta_block; //!< Number of bytes of meta information (padded) + uint32_t size_typename_block; //!< Number of bytes of typenames (padded) + + uint32_t size_label_block; //!< Number of bytes to store label strings (padded) + uint32_t size_scalar_block; //!< Number of bytes of scalar storage (padded) + uint32_t size_field_block_lo; //!< Number of bytes of field storage (padded) + uint32_t size_field_block_hi; //!< Number of bytes of field storage (padded) + + uint32_t size_crc_block; //!< Size of the trailing CRC field (unused!) + uint32_t size_padding_block; //!< Number of ignored bytes to pad to pagesize boundary. + uint32_t unused1; //!< Not used in current implementation. + uint32_t unused2; //!< Not used in current implementation. + + } header_t; + + typedef struct { + uint32_t type; //!< \brief Typecode for this type. + uint32_t elementsize; //!< \brief Number of bytes in an element + uint32_t count_lo; //!< \brief Number of elements (low) + uint32_t count_hi; //!< \brief Number of elements (high) + } metadisk_t; + + typedef struct key_prologue { + uint32_t magic; /* Magic number for frames */ + uint32_t frames_per_file; /* Number of frames in each file */ + uint32_t key_record_size; /* The size of each key record */ + } key_prologue_t; + + + //// utility routines + + /*! + * Extracts the low 32 bits of a 64 bit integer by masking. + */ + uint32_t lobytes(const uint64_t& x) { + uint32_t mask = 0xffffffff; + return x & mask; + } + + /*! + * Extract the high 32 bits of a 64 bit integer by shifting. + */ + uint32_t hibytes(const uint64_t& x) { + return x >> 32; + } + + /*! + * Extract the low 32 bits of a 64 bit float as an integer. + */ + uint32_t lobytes(const double& x) { + union { + uint64_t ival; + double dval; + } u; + u.dval = x; + return lobytes(u.ival); + } + + /*! + * Extract the high 32 bits of a 64 bit float as an integer. + */ + uint32_t hibytes(const double& x) { + union { + uint64_t ival; + double dval; + } u; + u.dval = x; + return hibytes(u.ival); + } + + /*! + * The byte order associated with this machine. We use + * 1234 for little endian, 4321 for big endian, and + * 3412 for the unlikely PDB endianism. + */ + uint32_t machineEndianism() { +#if __BYTE_ORDER == __LITTLE_ENDIAN + uint32_t byteorder = 1234; +#else +#if __BYTE_ORDER == __BIG_ENDIAN + uint32_t byteorder = 4321; +#else +#ifdef PDB_ENDIAN +#if __BYTE_ORDER == __PDB_ENDIAN + uint32_t byteorder = 3412; +#endif +#endif +#endif +#endif + // If we get a compile error here, then __BYTE_ORDER + // has an unexpected value. + return byteorder; + } + + uint64_t alignInteger( const uint64_t &x, unsigned border) { + return x + (border - x%border)%border; + } + + + /*! + * See RFC 1146 for Fletcher's Checksum (http://tools.ietf.org/html/rfc1146) + */ + uint32_t fletcher( const uint16_t *data, size_t len ) { + uint32_t sum1 = 0xffff, sum2 = 0xffff; + + while (len) { + unsigned tlen = len > 360 ? 360 : len; + len -= tlen; + do { + sum1 += *data++; + sum2 += sum1; + } while (--tlen); + sum1 = (sum1 & 0xffff) + (sum1 >> 16); + sum2 = (sum2 & 0xffff) + (sum2 >> 16); + } + /* Second reduction step to reduce sums to 16 bits */ + sum1 = (sum1 & 0xffff) + (sum1 >> 16); + sum2 = (sum2 & 0xffff) + (sum2 >> 16); + return sum2 << 16 | sum1; + } + + bool isfile(const std::string &name) { + struct stat statbuf; + return (stat(name.c_str(),&statbuf) == 0 && S_ISREG(statbuf.st_mode)); + } + + /*! + * Remove a file or directory. For directories, + * we recurse through subfiles and remove those + * before attempting the ::rmdir(); + */ + void recursivelyRemove(std::string path) { + struct stat statbuf; + + // ----------------------------------------------- + // Only try to unlink if the file exists + // We recurse through directories and unlink + // other files. + // ----------------------------------------------- + +#ifdef _WIN32 + // Use ::stat instead of ::lstat on windows since there are no symlinks + if (stat(path.c_str(),&statbuf) == 0) { +#else + if (::lstat(path.c_str(),&statbuf) == 0) { +#endif + if (!S_ISDIR(statbuf.st_mode)) { + if (::unlink(path.c_str()) != 0) { + throw std::runtime_error(strerror(errno)); + } + } else { + VMDDIR* directory = NULL; + try { + directory = vmd_opendir(path.c_str()); + if (directory) { + // Remove subfiles + char * entry; + while( (entry=vmd_readdir(directory)) != NULL ) { + // Don't unlink . or .. + if (entry[0] == '.') { + if (entry[1] == 0) continue; + if (entry[1] == '.' && entry[2] == 0) continue; + } + recursivelyRemove(path + s_sep + entry); + } + vmd_closedir(directory); + directory = NULL; + + // Remove the actual directory + if (::rmdir(path.c_str()) != 0) { + throw std::runtime_error(strerror(errno)); + } + } + } catch(...) { + if (directory) vmd_closedir(directory); + throw; + } + } + } + } + + + //////// + // CRC + //////// + + typedef uint32_t crc; + + #define POLYNOMIAL 0x04C11DB7 + #define WIDTH (8 * sizeof(crc)) + #define TOPBIT (1 << (WIDTH - 1)) + #define FINAL_XOR_VALUE 0xFFFFFFFF + + crc processByte( crc remainder, char msg ) { + remainder ^= (msg << (WIDTH - 8)); + for (uint8_t bit = 8; bit > 0; --bit) + { + if (remainder & TOPBIT) { + remainder = (remainder << 1) ^ POLYNOMIAL; + } else { + remainder = (remainder << 1); + } + } + return remainder; + } + + crc processBytes(const char *message, int nBytes) { + crc remainder = 0; + for (int byte = 0; byte < nBytes; ++byte) { + remainder = processByte( remainder, message[byte] ); + } + return remainder; + } + + int32_t cksum(const std::string &s) { + ssize_t len = s.size(); + int32_t result = processBytes( s.c_str(), len ); + + for ( ; len; len >>= 8) { + result = processByte( result, len & 0xff ); + } + return result ^ FINAL_XOR_VALUE; + } + +} + +bool Timekeys::init(const std::string& path ) { + std::string timekeys_path = path; + timekeys_path += s_sep; + timekeys_path += "timekeys"; + FILE * fd = fopen( timekeys_path.c_str(), "rb" ); + if (!fd) { + fprintf(stderr, "Could not find timekeys file at %s\n", + timekeys_path.c_str()); + return false; + } + + /* check the magic number */ + key_prologue_t prologue[1]; + if (fread( prologue, sizeof(key_prologue_t), 1, fd )!=1) { + fprintf(stderr, "Failed to read key prologue from %s\n", + timekeys_path.c_str()); + fclose(fd); + return false; + } + prologue->magic = htonl(prologue->magic); + if (prologue->magic != magic_timekey) { + fprintf(stderr, "timekeys magic number %x doesn't match %x\n", + prologue->magic, magic_timekey); + fclose(fd); + return false; + } + + /* get frames per file and key record size */ + prologue->frames_per_file = ntohl( prologue->frames_per_file ); + prologue->key_record_size = ntohl( prologue->key_record_size ); + m_fpf = prologue->frames_per_file; + + /* read all key records */ + fseek(fd, 0, SEEK_END); + off_t keyfile_size = ftello(fd); + size_t nframes = (keyfile_size-sizeof(key_prologue_t))/sizeof(key_record_t); + + keys.resize(nframes); + fseek(fd, sizeof(key_prologue_t), SEEK_SET); + if (fread(&keys[0], sizeof(key_record_t), nframes, fd)!=nframes) { + fprintf(stderr, "Failed to read all timekeys records: %s\n", + strerror(errno)); + fclose(fd); + return false; + } + fclose(fd); + + /* Check that we didn't get zero-length frames; this would be a strong + * indicator of file corruption! */ + int warning_count=0; + size_t i; + for (i=0; i 1e-3) { + if (getenv("DTRPLUGIN_VERBOSE")) fprintf(stderr, + "non-constant time interval at frame %zd\n", i); + return true; + } + /* constant offset */ + if (keys[i].offset() != m_framesize*( i % m_fpf)) { + fprintf(stderr, "unexpected offset for frame %zd\n", i); + return true; + } + } + /* looks good! Don't need the explicit key records anymore */ + keys.clear(); + return true; +} + +key_record_t Timekeys::operator[](uint64_t i) const { + if (i>m_fullsize) throw std::runtime_error("frame index out of range"); + if (keys.size()) return keys.at(i); + + key_record_t timekey; +#if defined(_MSC_VER) + double time = m_first + ((__int64) i)*m_interval; +#else + double time = m_first + i*m_interval; +#endif + uint64_t offset = (i % m_fpf) * m_framesize; + + timekey.time_lo = htonl(lobytes(time)); + timekey.time_hi = htonl(hibytes(time)); + timekey.offset_lo = htonl(lobytes(offset)); + timekey.offset_hi = htonl(hibytes(offset)); + timekey.framesize_lo = htonl(lobytes(m_framesize)); + timekey.framesize_hi = htonl(hibytes(m_framesize)); + return timekey; +} + +namespace { + template + void rawdump(std::ostream& out, const T& v) { + out.write((char *)&v, sizeof(v)); + } + + template + void rawload(std::istream& in, T& v) { + in.read((char *)&v, sizeof(v)); + } +} + +void Timekeys::dump(std::ostream& out) const { + rawdump(out, m_first); + rawdump(out, m_interval); + rawdump(out, m_framesize); + rawdump(out, m_size); + rawdump(out, m_fullsize); + rawdump(out, m_fpf); + rawdump(out, keys.size()); + if (keys.size()) { + out.write((const char *)&keys[0], keys.size()*sizeof(keys[0])); + } +} + +void Timekeys::load(std::istream& in) { + size_t sz; + rawload(in, m_first); + rawload(in, m_interval); + rawload(in, m_framesize); + rawload(in, m_size); + rawload(in, m_fullsize); + rawload(in, m_fpf); + rawload(in, sz); + if (sz) { + keys.resize(sz); + in.read((char *)&keys[0], keys.size()*sizeof(keys[0])); + } +} + +namespace { + struct Blob { + std::string type; + uint64_t count; + const void *data; + bool byteswap; + + Blob() : count(0), data(0) {} + Blob( const std::string &type_, uint64_t count_, const void *data_, + uint32_t frame_endianism ) + : type(type_), count(count_), data(data_), byteswap(false) { + uint32_t my_endianism = machineEndianism(); + if (frame_endianism != my_endianism) { + if ( (frame_endianism==1234 && my_endianism==4321) || + (frame_endianism==4321 && my_endianism==1234) ) { + byteswap=true; + } else { + throw std::runtime_error("Unable to handle frame endianness"); + } + } + } + + std::string str() const { + if (type=="char" && count>0) { + const char *s=(const char *)data; + return std::string(s, s+count); + } + return ""; + } + void get_float(float *buf) const { + if (type=="float") { + memcpy(buf, data, count*sizeof(float)); + } else if (type=="double") { + const double *p = reinterpret_cast(data); + for (uint64_t i=0; i(data); + for (uint64_t i=0; i BlobMap; +} + +static inline std::string addslash(const std::string& s){ + return (s.rbegin()[0] == '/') ? s : s + "/"; +} + +#define DD_RELPATH_MAXLEN (9) +static std::string +DDreldir(const std::string& fname, int ndir1, int ndir2){ + + if( fname.find('/', 0) != std::string::npos ) { + fprintf(stderr, "DDreldir: filename '%s' must not contain '/'\n", + fname.c_str()); + return ""; + } + + uint32_t hash = cksum(fname); + + // uint32_t u1 = ndir1; + // uint32_t u2 = ndir2; + uint32_t d1, d2; + char answer[DD_RELPATH_MAXLEN]; + if(ndir1 > 0){ + d1 = hash%ndir1; + if(ndir2 > 0){ + d2 = (hash/ndir1)%ndir2; + sprintf(answer, "%03x/%03x/", d1, d2); + }else{ + sprintf(answer, "%03x/", d1); + } + }else{ + sprintf(answer, "./"); + } + return std::string(answer); +} + +namespace { + class DDException : public std::runtime_error{ + public: + int eno; + DDException(const std::string &text, int _eno=0) + : std::runtime_error(text + strerror(_eno)), eno(_eno){} + }; +} + +void DDmkdir(const std::string &dirpath, mode_t mode, int ndir1, int ndir2){ + std::string dpslash(addslash(dirpath)); + + mode_t openmode = mode | 0300; // make sure we can write into the directory + if( mkdir(dpslash.c_str(), openmode) < 0 ) + throw DDException("mkdir", errno); + + if( mkdir((dpslash + "not_hashed").c_str(), openmode) < 0 ) + throw DDException("mkdir not_hashed subdirectory", errno); + + FILE *fp = fopen((dpslash + "not_hashed/.ddparams").c_str(), "w"); + if(fp == NULL) + throw DDException("fopen( .ddparams, \"w\" )", errno); + if( fprintf(fp, "%d %d\n", ndir1, ndir2) < 0 ){ + fclose(fp); + throw DDException("fprintf(.ddparams ...)", errno); + } + if( fclose(fp) ) + throw DDException("fclose(.ddparams)", errno); + + for(int i=0; i(mapping); + const header_t *header = reinterpret_cast(base); + if (lenrequired.magic) != magic_frame) { + char buf[256]; + sprintf(buf, "invalid magic number: expected %d, got %d\n", + magic_frame, ntohl(header->required.magic)); + throw std::runtime_error(buf); + } + + uint32_t size_header_block = ntohl(header->size_header_block); + uint32_t frames_endianism = ntohl(header->endianism); + uint32_t frames_nlabels = ntohl(header->nlabels); + uint32_t size_meta_block = ntohl(header->size_meta_block); + uint32_t size_typename_block = ntohl(header->size_typename_block); + uint32_t size_label_block = ntohl(header->size_label_block); + uint32_t size_scalar_block = ntohl(header->size_scalar_block); + uint32_t size_field_block_lo = ntohl(header->size_field_block_lo); + uint32_t size_field_block_hi = ntohl(header->size_field_block_hi); + uint64_t size_field_block = assemble64(size_field_block_lo, + size_field_block_hi); + + uint64_t offset_header_block = 0; + uint64_t offset_meta_block = offset_header_block + size_header_block; + uint64_t offset_typename_block = offset_meta_block + size_meta_block; + uint64_t offset_label_block = offset_typename_block + size_typename_block; + uint64_t offset_scalar_block = offset_label_block + size_label_block; + uint64_t offset_field_block = offset_scalar_block + size_scalar_block; + uint64_t offset_crc_block = offset_field_block + size_field_block; + + const metadisk_t* diskmeta = reinterpret_cast(base+offset_meta_block); + const char* typenames = reinterpret_cast(base+offset_typename_block); + const char* labels = reinterpret_cast(base+offset_label_block); + const char* scalars = reinterpret_cast(base+offset_scalar_block); + const char* fields = reinterpret_cast(base+offset_field_block); + const uint32_t * crc = reinterpret_cast(base+offset_crc_block); + if (*crc != 0) { + uint32_t frame_crc = fletcher(reinterpret_cast(base),offset_crc_block/2); + if (frame_crc != *crc) { + throw std::runtime_error("Checksum did not match"); + } + } + /* More sanity checks */ + if (len types; + while(*typenames) { + if (typenames >= labels) { + fprintf(stderr, "More typenames than labels!\n"); + break; + } + std::string type(typenames); + types.push_back(type); + typenames += type.size()+1; + } + + BlobMap blobs; + + for (size_t ii=0; iiinvmass.resize(natoms); + blob.get_float(&meta->invmass[0]); + } + } + free(meta_mapping); + close(meta_fd); + return meta; +} + +bool StkReader::recognizes(const std::string &path) { + return path.size()>4 && + path.substr(path.size()-4)==".stk" && + isfile(path); +} + +StkReader::StkReader(DtrReader *reader) { + dtr=reader->path(); + framesets.push_back(reader); + curframeset=0; +} + +bool StkReader::init(const std::string &path, int * changed) { + curframeset=0; + dtr=path; + + if (changed) *changed = 0; + /* process all the lines in the stk file */ + std::vector fnames; + std::ifstream input(path.c_str()); + if (!input) { + fprintf(stderr, "Cannot open '%s' for reading\n", path.c_str()); + return false; + } + std::string fname; + /* instantiate all the dtr readers */ + while (std::getline(input, fname)) { + fnames.push_back(fname); + } + if (!fnames.size()) { + fprintf(stderr, "Empty stk file\n"); + return false; + } + if (framesets.size()) { + /* reloading an stk. Find the dtrs that match the ones we already have, + * and discard the rest */ + unsigned i=0; /* i will become the index of the last dtr that we've + already loaded */ + for (; ipath()) break; + if (getenv("DTRPLUGIN_VERBOSE")) + fprintf(stderr, "StkReader: Reusing dtr at %s\n", + fnames[i].c_str()); + } + /* delete any remaining framesets */ + for (unsigned j=i; jkeys.restore_full_size(); + } + } + + /* instantiate dtr readers */ + for (unsigned i=0; i0) { + const DtrReader * first = framesets[0]; + /* reuse information from earlier readers */ + reader->natoms = first->natoms; + reader->with_velocity = first->with_velocity; + reader->set_meta(first->get_meta()); + } + if (!reader->init(fnames[i], NULL)) { + delete reader; + fprintf(stderr, "Failed opening frameset at %s\n", fnames[i].c_str()); + return false; + } + if (changed) *changed += 1; + framesets.push_back(reader); + if (i==0) this->with_velocity = reader->with_velocity; + } + + natoms=framesets[0]->natoms; + + // now remove overlaps + while (!framesets.back()->size()) { + delete framesets.back(); + framesets.pop_back(); + } + if (framesets.size()) { + double first=framesets.back()->keys[0].time(); + size_t i=framesets.size()-1; + while (i--) { + /* find out how many frames to keep in frameset[i] */ + Timekeys& cur = framesets[i]->keys; + size_t n = cur.size(); + while (n && cur[n-1].time() >= first) --n; + cur.truncate( n ); + if (cur.size()) { + double c0t = cur[0].time(); + first = (first < c0t) ? first : c0t; + } + } + } + return true; +} + +ssize_t StkReader::size() const { + ssize_t result=0; + for (size_t i=0; ikeys.size(); + return result; +} + +int StkReader::next(molfile_timestep_t *ts) { + int rc=MOLFILE_EOF; + while (curframeset < framesets.size() && + (rc=framesets[curframeset]->next(ts))==MOLFILE_EOF) { + ++curframeset; + } + return rc; +} + +const DtrReader * StkReader::component(ssize_t &n) const { + for (size_t i=0; isize(); + if (n < size) return framesets[i]; + n -= size; + } + return NULL; +} + +int StkReader::frame(ssize_t n, molfile_timestep_t *ts) const { + const DtrReader *comp = component(n); + if (!comp) return MOLFILE_EOF; + return comp->frame(n, ts); +} + +StkReader::~StkReader() { + for (size_t i=0; i0 && natoms==0) { + if (getenv("DTRPLUGIN_VERBOSE")) { + fprintf(stderr, "reading first frame to get atom count\n"); + } + std::string fname=::framefile(dtr, 0, keys.framesperfile(), + ndir1(), ndir2()); + int fd = open(fname.c_str(), O_RDONLY|O_BINARY); + ssize_t framesize=0; + unsigned i; + void *mapping = read_file( fd, 0, framesize ); + if (mapping==NULL) { + fprintf(stderr, "Failed to find frame at %s\n", fname.c_str()); + close(fd); + return false; + } + BlobMap blobs; + try { + blobs = read_frame(mapping, framesize); + } + catch (std::exception &e) { + fprintf(stderr, "Warning: reading first frame failed, %s", e.what()); + } + with_momentum = blobs.find("MOMENTUM")!=blobs.end(); + + // I'm aware of three possible sources of positions: + // "POSN" (the original frameset format) + // "POSITION" (the wrapped frameset formats) + // "POS" (anton trajectories) + const char *posnames[] = { "POSN", "POSITION", "POS" }; + for (i=0; i<3; i++) { + if (blobs.find(posnames[i])!=blobs.end()) { + natoms = blobs[posnames[i]].count / 3; + break; + } + } + // similar for velocities + const char *velnames[] = { "MOMENTUM", "VELOCITY" }; + for (i=0; i<2; i++) { + if (blobs.find(velnames[i])!=blobs.end()) { + with_velocity=true; + break; + } + } + free(mapping); + close(fd); + } + + if (natoms>0 && meta==NULL && owns_meta==false) { + meta=read_meta( dtr + s_sep + "metadata",natoms, with_momentum ); + owns_meta = true; + } + + /* we always reread the timekeys */ + if (changed) *changed = 1; + return true; + +} + +ssize_t DtrReader::times(ssize_t start, ssize_t count, double *t) const { + ssize_t remaining = keys.size()-start; + count = (count < remaining) ? count : remaining; + for (ssize_t j=0; jsize(); + if (starttimes(start, count, t+nread); + nread += sz; + count -= sz; + start=0; + if (!count) break; + } + return nread; +} + +static double dotprod(const double *x, const double *y) { + return x[0]*y[0] + x[1]*y[1] + x[2]*y[2]; +} + +static void read_homebox( const double *box, + molfile_timestep_t *ts ) { + + ts->A = ts->B = ts->C = 0; + + double A[3] = { box[0], box[3], box[6] }; + double B[3] = { box[1], box[4], box[7] }; + double C[3] = { box[2], box[5], box[8] }; + + // store lengths + ts->A = sqrt(dotprod(A,A)); + ts->B = sqrt(dotprod(B,B)); + ts->C = sqrt(dotprod(C,C)); + + if (ts->A == 0 || ts->B == 0 || ts->C == 0) { + + ts->alpha = ts->beta = ts->gamma = 90; + + } else { + + // compute angles + double cosAB = dotprod(A,B)/(ts->A * ts->B); + double cosAC = dotprod(A,C)/(ts->A * ts->C); + double cosBC = dotprod(B,C)/(ts->B * ts->C); + + // clamp + if (cosAB > 1.0) cosAB = 1.0; else if (cosAB < -1.0) cosAB = -1.0; + if (cosAC > 1.0) cosAC = 1.0; else if (cosAC < -1.0) cosAC = -1.0; + if (cosBC > 1.0) cosBC = 1.0; else if (cosBC < -1.0) cosBC = -1.0; + + // convert to angles using asin to avoid nasty rounding when we are + // close to 90 degree angles. + ts->alpha = 90.0 - asin(cosBC) * 90.0 / M_PI_2; /* cosBC */ + ts->beta = 90.0 - asin(cosAC) * 90.0 / M_PI_2; /* cosAC */ + ts->gamma = 90.0 - asin(cosAB) * 90.0 / M_PI_2; /* cosAB */ + } +} + +void write_homebox( const molfile_timestep_t * ts, + float * box ) { + + double A[3], B[3], C[3]; + + // Convert VMD's unit cell information + double cosBC = sin( ((90 - ts->alpha ) / 180) * M_PI ); + double cosAC = sin( ((90 - ts->beta ) / 180) * M_PI ); + double cosAB = sin( ((90 - ts->gamma ) / 180) * M_PI ); + double sinAB = cos( ((90 - ts->gamma ) / 180) * M_PI ); + + double Ax = ts->A; + double Ay = 0; + double Az = 0; + double Bx = ts->B * cosAB; + double By = ts->B * sinAB; + double Bz = 0; + double Cx,Cy,Cz; + if (sinAB != 0) { + Cx = cosAC; + Cy = (cosBC - cosAC*cosAB) / sinAB; + Cz = sqrt(1-Cx*Cx-Cy*Cy); + Cx *= ts->C; + Cy *= ts->C; + Cz *= ts->C; + } else { + Cx=Cy=Cz=0; + } + A[0] = Ax; A[1] = Ay; A[2] = Az; + B[0] = Bx; B[1] = By; B[2] = Bz; + C[0] = Cx; C[1] = Cy; C[2] = Cz; + + // put vectors in column of homebox + box[0] = A[0]; box[3] = A[1]; box[6] = A[2]; + box[1] = B[0]; box[4] = B[1]; box[7] = B[2]; + box[2] = C[0]; box[5] = C[1]; box[8] = C[2]; +} + +static int handle_wrapped_v2( + BlobMap &blobs, + uint32_t natoms, + bool with_velocity, + molfile_timestep_t *ts ) { + + // just read POSITION in either single or double precision + if (blobs.find("POSITION")==blobs.end()) { + fprintf(stderr, "ERROR, Missing POSITION field in frame\n"); + return MOLFILE_ERROR; + } + Blob pos=blobs["POSITION"]; + if (pos.count != 3*natoms) { + fprintf(stderr, "ERROR, Expected %d elements in POSITION; got %" PRIu64 "\n", + 3*natoms, pos.count); + return MOLFILE_ERROR; + } + pos.get_float(ts->coords); + + if (with_velocity && ts->velocities && blobs.find("VELOCITY")!=blobs.end()) { + Blob vel=blobs["VELOCITY"]; + if (vel.count != 3*natoms) { + fprintf(stderr, "ERROR, Expected %d elements in VELOCITY; got %" PRIu64 "\n", + 3*natoms, vel.count); + return MOLFILE_ERROR; + } + vel.get_float(ts->velocities); + } + + if (blobs.find("UNITCELL")!=blobs.end()) { + double box[9]; + blobs["UNITCELL"].get_double(box); + read_homebox( box, ts ); + } + +#if defined(DESRES_READ_TIMESTEP2) + if (blobs.find("ENERGY")!=blobs.end()) { + blobs["ENERGY"].get_double(&ts->total_energy); + } + + if (blobs.find("POT_ENERGY")!=blobs.end()) { + blobs["POT_ENERGY"].get_double(&ts->potential_energy); + } + + if (blobs.find("KIN_ENERGY")!=blobs.end()) { + blobs["KIN_ENERGY"].get_double(&ts->kinetic_energy); + } + + if (blobs.find("EX_ENERGY")!=blobs.end()) { + blobs["EX_ENERGY"].get_double(&ts->extended_energy); + } + + if (blobs.find("PRESSURE")!=blobs.end()) { + blobs["PRESSURE"].get_double(&ts->pressure); + } + + if (blobs.find("TEMPERATURE")!=blobs.end()) { + blobs["TEMPERATURE"].get_double(&ts->temperature); + } +#endif + + return MOLFILE_SUCCESS; +} + +namespace { + + inline void + compute_center(int partition, + int nx, int ny, int nz, + float b0, float b1, float b2, + float b3, float b4, float b5, + float b6, float b7, float b8, + float* cx, float* cy, float* cz) { + double nu, nv, nw, mu, mv, mw; + double xc, yc, zc; + + // ----------------------------------------------- + // Map the partition number to its "mesh" position + // (see define_mesh_collective in topology.c) + // ----------------------------------------------- + int hmx = partition; + int hmy = hmx / nx; /* y = y + ny*( z + nz*r ) */ + int hmz = hmy / ny; /* z = z + nz*r */ + hmx -= hmy * nx; /* x = x */ + hmy -= hmz * ny; /* y = y */ + + nu = (double)nx; + nv = (double)ny; + nw = (double)nz; + + // ----------------------------------------------- + // Code adapted from configure_global_cell in + // topology.c + // ----------------------------------------------- + mu = -0.5*(nu-1) + (double)hmx; + mv = -0.5*(nv-1) + (double)hmy; + mw = -0.5*(nw-1) + (double)hmz; + + // We used to do FORCE_PRECISION(xc,float) here, but that + // seems unnecessary in the context of trajectory writing. + xc = b0*mu + b1*mv + b2*mw; + yc = b3*mu + b4*mv + b5*mw; + zc = b6*mu + b7*mv + b8*mw; + + *cx = xc; + *cy = yc; + *cz = zc; + } + + inline int + posn_momentum_v_1(int32_t nx, int32_t ny, int32_t nz, + uint64_t nparticles, + const double * home_box, + const uint32_t* gid, + const uint32_t* npp, + const float * rm, // reciprocal mass + const float* posn, const float* momentum, + /* returns */ + float *position, float *velocity, double *box) { + + // bounding box is a straight multiple of the home box + if (box) { + box[0] = home_box[0]*nx; + box[1] = home_box[1]*ny; + box[2] = home_box[2]*nz; + + box[3] = home_box[3]*nx; + box[4] = home_box[4]*ny; + box[5] = home_box[5]*nz; + + box[6] = home_box[6]*nx; + box[7] = home_box[7]*ny; + box[8] = home_box[8]*nz; + } + + + int partition = 0; + int remaining = 0; + float cx = 0; + float cy = 0; + float cz = 0; + float ux = home_box[0]; + float vx = home_box[1]; + float wx = home_box[2]; + float uy = home_box[3]; + float vy = home_box[4]; + float wy = home_box[5]; + float uz = home_box[6]; + float vz = home_box[7]; + float wz = home_box[8]; + + for(uint64_t i=0; i= nparticles) { + fprintf(stderr, "non-contiguous particles\n"); + return MOLFILE_ERROR; + } + + if (posn) { + float x = posn[3*i+0]; + float y = posn[3*i+1]; + float z = posn[3*i+2]; + + position[3*id+0] = ux*x + vx*y + wx*z + cx; + position[3*id+1] = uy*x + vy*y + wy*z + cy; + position[3*id+2] = uz*x + vz*y + wz*z + cz; + } + + if (velocity && momentum && rm) { + velocity[3*id+0] = momentum[3*i+0]*rm[id]; + velocity[3*id+1] = momentum[3*i+1]*rm[id]; + velocity[3*id+2] = momentum[3*i+2]*rm[id]; + } else if (velocity) { + velocity[3*id+0] = 0.0; + velocity[3*id+1] = 0.0; + velocity[3*id+2] = 0.0; + } + --remaining; + } + return MOLFILE_SUCCESS; + } +} + +static int handle_posn_momentum_v1( + BlobMap &blobs, + uint32_t natoms, + bool with_velocity, + const float * rmass, + molfile_timestep_t *ts ) { + + int32_t nx, ny, nz; + double home_box[9], box[9]; + blobs["HOME_BOX"].get_double(home_box); + blobs["NX"].get_int32(&nx); + blobs["NY"].get_int32(&ny); + blobs["NZ"].get_int32(&nz); + + std::vector gid, npp; + std::vector pos, mtm; + Blob gidblob=blobs["GID"]; + Blob nppblob=blobs["NPP"]; + Blob posblob=blobs["POSN"]; + Blob mtmblob=blobs["MOMENTUM"]; + + if (gidblob.count != natoms) { + fprintf(stderr, "Missing GID field\n"); + return MOLFILE_ERROR; + } + if (posblob.count != 3*natoms) { + fprintf(stderr, "Missing POSN field\n"); + return MOLFILE_ERROR; + } + gid.resize(gidblob.count); + npp.resize(nppblob.count); + pos.resize(posblob.count); + mtm.resize(mtmblob.count); + + gidblob.get_uint32(&gid[0]); + nppblob.get_uint32(&npp[0]); + posblob.get_float(&pos[0]); + + if (rmass && with_velocity) mtmblob.get_float(&mtm[0]); + + posn_momentum_v_1( nx, ny, nz, natoms, home_box, + &gid[0], &npp[0], rmass, + &pos[0], + &mtm[0], + ts->coords, + ts->velocities, + box ); + + read_homebox( box, ts ); + return MOLFILE_SUCCESS; +} + +static int handle_wrapped_v1( + BlobMap &blobs, + uint32_t natoms, + bool with_velocity, + molfile_timestep_t *ts ) { + + { + // homebox + double home_box[9], box[9]; + int32_t nx, ny, nz; + blobs["HOME_BOX"].get_double(home_box); + blobs["NX"].get_int32(&nx); + blobs["NY"].get_int32(&ny); + blobs["NZ"].get_int32(&nz); + box[0] = home_box[0]*nx; + box[1] = home_box[1]*ny; + box[2] = home_box[2]*nz; + + box[3] = home_box[3]*nx; + box[4] = home_box[4]*ny; + box[5] = home_box[5]*nz; + + box[6] = home_box[6]*nx; + box[7] = home_box[7]*ny; + box[8] = home_box[8]*nz; + read_homebox( box, ts ); + } + + Blob posblob=blobs["POSN"]; + Blob velblob=blobs["VELOCITY"]; + + // get positions + if (posblob.count != 3*natoms) { + fprintf(stderr, "Missing POSN field\n"); + return MOLFILE_ERROR; + } + posblob.get_float(ts->coords); + + // if required, get velocities + if (ts->velocities && velblob.count > 0) { + if (velblob.count != 3*natoms) { + fprintf(stderr, "VELOCITY field has %" PRIu64 " values; expected %d\n", + velblob.count, 3*natoms); + return MOLFILE_ERROR; + } + velblob.get_float(ts->velocities); + } + return MOLFILE_SUCCESS; +} + +static int handle_anton_sfxp_v3( + BlobMap &blobs, + uint32_t natoms, + bool with_velocity, + const float * rmass, + molfile_timestep_t *ts ) { + + if (!rmass) { + fprintf(stderr, "Cannot read anton_sfxp_v3 frame without rmass\n"); + return MOLFILE_ERROR; + } + + double positionScale=0, momentumScale=0; + // position scale... + { + Blob blob = blobs["POSITIONSCALE"]; + if (blob.count != 1) { + fprintf(stderr, "Missing POSITIONSCALE field\n"); + return MOLFILE_ERROR; + } + blob.get_double(&positionScale); + } + // momentum scale + if (ts->velocities) { + Blob blob = blobs["MOMENTUMSCALE"]; + if (blob.count != 1) { + fprintf(stderr, "Missing MOMENTUMSCALE field\n"); + return MOLFILE_ERROR; + } + blob.get_double(&momentumScale); + momentumScale *= PEAKmassInAmu; + } + + // box + { + double box[9] = { 0,0,0, 0,0,0, 0,0,0 }; + uint32_t anton_box[3]; + Blob boxblob = blobs["BOX"]; + if (boxblob.count != 3) { + fprintf(stderr, "Missing BOX field\n"); + return MOLFILE_ERROR; + } + boxblob.get_uint32(anton_box); + box[0] = sfxp_ulp32flt(anton_box[0])*positionScale; + box[4] = sfxp_ulp32flt(anton_box[1])*positionScale; + box[8] = sfxp_ulp32flt(anton_box[2])*positionScale; + read_homebox( box, ts ); + } + + // velocities + std::vector vel; + if (ts->velocities) { + Blob velblob = blobs["MOMENTUM"]; + if (velblob.count != 3*natoms) { + fprintf(stderr, "Missing MOMENTUM field\n"); + return MOLFILE_ERROR; + } + vel.resize(3*natoms); + velblob.get_int32(&vel[0]); + } + + // positions + std::vector pos(3*natoms); + { + Blob posblob = blobs["POS"]; + if (posblob.count != 3*natoms) { + fprintf(stderr, "Missing POS field\n"); + return MOLFILE_ERROR; + } + posblob.get_int32(&pos[0]); + } + // convert and read into supplied storage + for (unsigned i=0; icoords[3*i ] = sfxp_ulp32flt(pos[3*i+0])*positionScale; + ts->coords[3*i+1] = sfxp_ulp32flt(pos[3*i+1])*positionScale; + ts->coords[3*i+2] = sfxp_ulp32flt(pos[3*i+2])*positionScale; + if (ts->velocities) { + const double rm = rmass[i] * momentumScale; // includes PEAKmassInAmu + ts->velocities[3*i ] = (float)(rm * sfxp_ulp32flt(vel[3*i ])); + ts->velocities[3*i+1] = (float)(rm * sfxp_ulp32flt(vel[3*i+1])); + ts->velocities[3*i+2] = (float)(rm * sfxp_ulp32flt(vel[3*i+2])); + } + } + return MOLFILE_SUCCESS; +} + +int DtrReader::next(molfile_timestep_t *ts) { + + if (eof()) return MOLFILE_EOF; + if (!ts) { + ++m_curframe; + return MOLFILE_SUCCESS; + } + ssize_t iframe = m_curframe; + ++m_curframe; + return frame(iframe, ts); +} + +int DtrReader::ndir1() const { + if (m_ndir1<0) DDgetparams(dtr, &m_ndir1, &m_ndir2); + return m_ndir1; +} + +int DtrReader::ndir2() const { + if (m_ndir2<0) DDgetparams(dtr, &m_ndir1, &m_ndir2); + return m_ndir2; +} + +int DtrReader::frame(ssize_t iframe, molfile_timestep_t *ts) const { + int rc = MOLFILE_SUCCESS; + { + off_t offset=0; + ssize_t framesize=0; + if (framesperfile() != 1) { + offset = assemble64( ntohl(keys[iframe].offset_lo), + ntohl(keys[iframe].offset_hi) ); + framesize = assemble64( ntohl(keys[iframe].framesize_lo), + ntohl(keys[iframe].framesize_hi) ); + + } + ts->physical_time = keys[iframe].time(); + std::string fname=::framefile(dtr, iframe, framesperfile(), ndir1(), ndir2()); + int fd = open(fname.c_str(), O_RDONLY|O_BINARY); + if (fd<0) return MOLFILE_EOF; + void *mapping = read_file( fd, offset, framesize ); + if (mapping==NULL) { + close(fd); + return MOLFILE_EOF; + } + + rc = frame_from_bytes( mapping, framesize, ts ); + + free(mapping); + close(fd); + } + return rc; +} + +int DtrReader::frame_from_bytes(const void *buf, uint64_t len, + molfile_timestep_t *ts) const { + + BlobMap blobs; + try { + blobs = read_frame(buf, len); + } + catch (std::exception &e) { + fprintf(stderr, "Reading frame failed: %s\n", e.what()); + return MOLFILE_ERROR; + } + + const float * rmass = meta && meta->invmass.size() ? + &meta->invmass[0] : NULL; + + // Now, dispatch to routines based on format + std::string format = blobs["FORMAT"].str(); + if (format=="WRAPPED_V_2" || format == "DBL_WRAPPED_V_2") { + return handle_wrapped_v2(blobs, natoms, with_velocity, ts); + + } else if (format=="POSN_MOMENTUM_V_1" || format=="DBL_POSN_MOMENTUM_V_1") { + return handle_posn_momentum_v1(blobs, natoms, with_velocity, rmass, ts); + + } else if (format=="WRAPPED_V_1" || format == "DBL_WRAPPED_V_1") { + return handle_wrapped_v1(blobs, natoms, with_velocity, ts); + + } else if (format=="ANTON_SFXP_V3") { + return handle_anton_sfxp_v3(blobs, natoms, with_velocity, rmass, ts); + } + fprintf(stderr, "ERROR, can't handle format %s\n", format.c_str()); + return MOLFILE_ERROR; +} + + +namespace { + struct meta_t { + std::string label; + std::string typecode; + uint32_t elementsize; + uint64_t count; + const char *bytes; + meta_t() {} + meta_t(const std::string &l, const std::string &t, uint32_t e, uint32_t c, + const void *b) + : label(l), typecode(t), elementsize(e), count(c), + bytes(reinterpret_cast(b)) {} + }; + typedef std::vector MetaList; + + uint64_t typename_size(const MetaList &meta) { + // just the set of distinct types + uint64_t sz=0; + typedef std::set Typemap; + Typemap types; + for (MetaList::const_iterator m=meta.begin(); m!=meta.end(); ++m) + types.insert(m->typecode); + for (Typemap::const_iterator s=types.begin(); s!=types.end();++s) + sz += s->size() + 1; + sz += 1; + return alignInteger(sz, s_alignsize); + } + + uint64_t label_size(const MetaList &meta) { + uint64_t sz=0; + for (MetaList::const_iterator m=meta.begin(); m!=meta.end(); ++m) + sz += m->label.size() + 1; + sz += 1; + return alignInteger(sz, s_alignsize); + } + + uint64_t scalar_size(const MetaList &meta) { + uint64_t sz=0; + for (MetaList::const_iterator m=meta.begin(); m!=meta.end(); ++m) + if (m->count <= 1) + sz += alignInteger( m->elementsize * m->count, s_alignsize ); + return sz; + } + uint64_t field_size(const MetaList &meta) { + uint64_t sz=0; + for (MetaList::const_iterator m=meta.begin(); m!=meta.end(); ++m) + if (m->count > 1) + sz += alignInteger( m->elementsize * m->count, s_alignsize ); + return sz; + } + + void construct_frame( const std::vector& meta, + std::vector& bytes ) { + uint64_t offset_header_block = 0; + uint64_t size_header_block = + alignInteger( sizeof(header_t), s_alignsize ); + + uint64_t offset_meta_block = offset_header_block + size_header_block; + uint64_t size_meta_block = + alignInteger( meta.size()*sizeof(metadisk_t), s_alignsize ); + + uint64_t offset_typename_block = offset_meta_block + size_meta_block; + uint64_t size_typename_block = typename_size(meta); + + uint64_t offset_label_block = offset_typename_block + size_typename_block; + uint64_t size_label_block = label_size(meta); + + uint64_t offset_scalar_block = offset_label_block + size_label_block; + uint64_t size_scalar_block = scalar_size(meta); + + uint64_t offset_field_block = offset_scalar_block + size_scalar_block; + uint64_t size_field_block = field_size(meta); + + uint64_t offset_crc_block = offset_field_block + size_field_block; + uint64_t size_crc_block = sizeof(uint32_t); + + uint64_t offset_padding_block = offset_crc_block + size_crc_block; + uint64_t size_padding_block = + alignInteger(offset_padding_block,s_blocksize) - offset_padding_block; + + uint64_t framesize = offset_padding_block + size_padding_block; + + // construct the frame + bytes.resize(framesize); + char * base = &bytes[0]; + memset( base, 0, framesize ); + + header_t *header = reinterpret_cast(base+offset_header_block); + metadisk_t* diskmeta = reinterpret_cast(base+offset_meta_block); + char* typenames = reinterpret_cast(base+offset_typename_block); + char* labels = reinterpret_cast(base+offset_label_block); + char* scalars = reinterpret_cast(base+offset_scalar_block); + char* fields = reinterpret_cast(base+offset_field_block); + uint32_t* crc = reinterpret_cast(base+offset_crc_block); + //char* padding = reinterpret_cast(base+offset_padding_block); + + /*** header ***/ + memset(header,0,sizeof(header_t)); + header->required.magic = htonl(magic_frame); + header->required.version = htonl(s_version); + + header->required.framesize_lo = htonl(lobytes(framesize)); + header->required.framesize_hi = htonl(hibytes(framesize)); + + header->size_header_block = htonl(size_header_block); + header->unused0 = 0; + uint64_t lrosetta = assemble64(s_lrosetta_lo,s_lrosetta_hi); + header->irosetta = s_irosetta; + header->frosetta = s_frosetta; + + header->drosetta_lo = lobytes(s_drosetta); + header->drosetta_hi = hibytes(s_drosetta); + header->lrosetta_lo = lobytes(lrosetta); + header->lrosetta_hi = hibytes(lrosetta); + + header->endianism = htonl(machineEndianism()); + header->nlabels = htonl(meta.size()); + header->size_meta_block = htonl(size_meta_block); + header->size_typename_block = htonl(size_typename_block); + + header->size_label_block = htonl(size_label_block); + header->size_scalar_block = htonl(size_scalar_block); + header->size_field_block_lo = htonl(lobytes(size_field_block)); + header->size_field_block_hi = htonl(hibytes(size_field_block)); + + header->size_crc_block = htonl(size_crc_block); + header->size_padding_block = htonl(size_padding_block); + header->unused1 = 0; + header->unused2 = 0; + + std::map typemap; + + for (MetaList::const_iterator m=meta.begin(); m!=meta.end(); ++m) { + + if (typemap.find(m->typecode)==typemap.end()) { + unsigned code=typemap.size(); + typemap[m->typecode]=code; + typenames=std::copy(m->typecode.begin(), m->typecode.end(), typenames); + *typenames++ = 0; + } + + diskmeta->type = htonl( typemap[m->typecode] ); + diskmeta->elementsize = htonl( m->elementsize ); + diskmeta->count_lo = htonl( lobytes( m->count )); + diskmeta->count_hi = htonl( hibytes( m->count )); + diskmeta++; + + labels=std::copy(m->label.begin(), m->label.end(), labels); + *labels++ = 0; + + uint64_t nbytes = m->count*m->elementsize; + if (m->count <= 1) { + memcpy( scalars, m->bytes, nbytes ); + scalars += alignInteger( nbytes, s_alignsize ); + } else { + memcpy( fields, m->bytes, nbytes ); + fields += alignInteger( nbytes, s_alignsize ); + } + } + *crc = fletcher(reinterpret_cast(base),offset_crc_block/2); + } +} + +void write_all( int fd, const char * buf, ssize_t count ) { + while (count) { + ssize_t n = ::write(fd, buf, count); + if (n<0) { + if (errno==EINTR) continue; + throw std::runtime_error(strerror(errno)); + } + buf += n; + count -= n; + } +} + +bool DtrWriter::init(const std::string &path) { + + dtr=path; + try { + m_directory=path; + char cwd[4096]; + + while(m_directory.size() > 0 && m_directory[m_directory.size()-1] == s_sep) { + m_directory.erase(m_directory.size()-1); + } + + if ( m_directory[0] != s_sep) { + if (! ::getcwd(cwd,sizeof(cwd))) { + throw std::runtime_error(strerror(errno)); + } + m_directory = std::string(cwd) + s_sep + m_directory; + } + + recursivelyRemove(m_directory); + ::DDmkdir(m_directory,0777,0, 0); + + // craft an empty metadata frame + std::vector meta; + std::vector bytes; + construct_frame( meta, bytes ); + + { + std::string metadata_file = m_directory + s_sep + "metadata"; + FILE *fd = fopen(metadata_file.c_str(), "wb"); + fwrite( &bytes[0], bytes.size(), 1, fd ); + fclose(fd); + } + + // start writing timekeys file */ + std::string timekeys_path = dtr + s_sep + "timekeys"; + timekeys_file = fopen( timekeys_path.c_str(), "wb" ); + if (!timekeys_file) { + fprintf(stderr, "Opening timekeys failed: %s\n", strerror(errno)); + return false; + } else { + key_prologue_t prologue[1]; + prologue->magic = htonl(magic_timekey); + prologue->frames_per_file = htonl(frames_per_file); + prologue->key_record_size = htonl(sizeof(key_record_t)); + fwrite( prologue, sizeof(key_prologue_t), 1, timekeys_file ); + } + } + catch (std::exception &e) { + fprintf(stderr, "%s\n", e.what()); + return false; + } + return true; +} + +int DtrWriter::next(const molfile_timestep_t *ts) { + + try { + static const char *format = "WRAPPED_V_2"; + static const char *title = "written by VMD"; + float box[9]; + write_homebox( ts, box ); + + double time = ts->physical_time; + + /* require increasing times (DESRESCode#1053) */ + if (last_time != HUGE_VAL && time <= last_time) { + fprintf(stderr, + "dtrplugin: framesets require increasing times. previous %e, current %e\n", + last_time, time); + return MOLFILE_ERROR; + } + + std::vector meta; + meta.push_back( + meta_t( "FORMAT", "char", sizeof(char), strlen(format), format )); + meta.push_back( + meta_t( "TITLE", "char", sizeof(char), strlen(title), title )); + meta.push_back( + meta_t( "CHEMICAL_TIME", "double", sizeof(double), 1, &time)); + meta.push_back( + meta_t( "UNITCELL", "float", sizeof(float), 9, box )); + meta.push_back( + meta_t( "POSITION", "float", sizeof(float), 3*natoms, ts->coords )); + if (ts->velocities) meta.push_back( + meta_t( "VELOCITY", "float", sizeof(float), 3*natoms, ts->velocities )); +#if defined(DESRES_READ_TIMESTEP2) + meta.push_back( + meta_t( "ENERGY", "double", sizeof(double), 1, &ts->total_energy )); + meta.push_back( + meta_t( "POT_ENERGY", "double", sizeof(double), 1, &ts->potential_energy )); + meta.push_back( + meta_t( "KIN_ENERGY", "double", sizeof(double), 1, &ts->kinetic_energy )); + meta.push_back( + meta_t( "EX_ENERGY", "double", sizeof(double), 1, &ts->extended_energy )); + meta.push_back( + meta_t( "TEMPERATURE", "double", sizeof(double), 1, &ts->temperature)); + meta.push_back( + meta_t( "PRESSURE", "double", sizeof(double), 1, &ts->pressure)); +#endif + + std::vector base; + construct_frame(meta, base); + uint64_t framesize = base.size(); + uint64_t keys_in_file = nwritten % frames_per_file; + + if (!keys_in_file) { + if (frame_fd>0) ::close(frame_fd); + + framefile_offset = 0; + std::string filepath=framefile(dtr, nwritten, frames_per_file, 0, 0); + frame_fd = open(filepath.c_str(),O_WRONLY|O_CREAT|O_BINARY,0666); + if (frame_fd<0) throw std::runtime_error(strerror(errno)); + } + + // write the data to disk + write_all( frame_fd, &base[0], framesize ); + + // add an entry to the keyfile list + key_record_t timekey; + timekey.time_lo = htonl(lobytes(time)); + timekey.time_hi = htonl(hibytes(time)); + timekey.offset_lo = htonl(lobytes(framefile_offset)); + timekey.offset_hi = htonl(hibytes(framefile_offset)); + timekey.framesize_lo = htonl(lobytes(framesize)); + timekey.framesize_hi = htonl(hibytes(framesize)); + + if (fwrite(&timekey, sizeof(timekey), 1, timekeys_file)!=1) { + fprintf(stderr, "Writing timekey failed\n"); + return MOLFILE_ERROR; + } + +#if defined(_MSC_VER) || defined(__MINGW32__) + _commit(frame_fd); +#else + fsync(frame_fd); +#endif + fflush(timekeys_file); +#if defined(_MSC_VER) || defined(__MINGW32__) + _commit(fileno(timekeys_file)); +#else + fsync(fileno(timekeys_file)); +#endif + + ++nwritten; + framefile_offset += framesize; + } + catch (std::exception &e) { + fprintf(stderr, "Write failed: %s\n",e.what()); + return MOLFILE_ERROR; + } + return MOLFILE_SUCCESS; +} + +DtrWriter::~DtrWriter() { + if (frame_fd>0) ::close(frame_fd); + if (timekeys_file) fclose(timekeys_file); +} + +/* compressed form: first write a -1 to indicate the new format. Then + * write number of ranges n, followed by (start, count) pairs. */ +std::ostream& operator<<(std::ostream& out, const metadata_t& meta) { + out << meta.invmass.size() << ' '; + if (meta.invmass.size()) { + out.write( (const char *)&meta.invmass[0], meta.invmass.size()*sizeof(meta.invmass[0])); + } + return out; +} + +std::istream& operator>>(std::istream& in, metadata_t& meta) { + uint32_t sz; + char c; + in >> sz; + in.get(c); + meta.invmass.resize(sz); + if (sz) { + in.read((char *)&meta.invmass[0], sz*sizeof(meta.invmass[0])); + } + return in; +} + +std::ostream& DtrReader::dump(std::ostream &out) const { + bool has_meta = meta!=NULL; + out << SERIALIZED_VERSION << ' ' + << dtr << ' ' + << natoms << ' ' + << with_velocity << ' ' + << owns_meta << ' ' + << has_meta << ' '; + if (owns_meta && has_meta) { + out << *meta; + } + /* write raw m_ndir values so that we don't read them from .ddparams + * if they haven't been read yet */ + out << m_ndir1 << ' ' + << m_ndir2 << ' '; + keys.dump(out); + return out; +} +std::istream& DtrReader::load(std::istream &in) { + char c; + bool has_meta; + std::string version; + in >> version; + if (badversion(version)) { + fprintf(stderr, "Bad version string\n"); + in.setstate( std::ios::failbit ); + return in; + } + in >> dtr + >> natoms + >> with_velocity + >> owns_meta + >> has_meta; + if (owns_meta && has_meta) { + delete meta; + meta = new metadata_t; + in.get(c); + in >> *meta; + } + in >> m_ndir1 + >> m_ndir2; + in.get(c); + keys.load(in); + return in; +} + +std::ostream& StkReader::dump(std::ostream &out) const { + out << dtr << ' ' + << framesets.size() << ' '; + for (size_t i=0; idump(out); + return out; +} + +std::istream& StkReader::load(std::istream &in) { + in >> dtr; + size_t size; in >> size; framesets.resize(size); + char c; in.get(c); + with_velocity=false; + for (size_t i=0; iload(in); + if (i>0) framesets[i]->set_meta(framesets[0]->get_meta()); + if (i==0) with_velocity=framesets[i]->with_velocity; + } + if (framesets.size()) natoms=framesets[0]->natoms; + return in; +} +#ifdef ENABLE_VMD_PLUGIN +/////////////////////////////////////////////////////////////////// +// +// Plugin Interface +// +// //////////////////////////////////////////////////////////////// + +static void * +open_file_read( const char *filename, const char *filetype, int *natoms ) { + + FrameSetReader *h = NULL; + std::string fname; + + // check for .stk file + if (StkReader::recognizes(filename)) { + h = new StkReader; + + } else { + h = new DtrReader; + // check for "clickme.dtr" + fname=filename; + std::string::size_type pos = fname.rfind( "clickme.dtr" ); + if (pos != std::string::npos) { + fname.resize(pos); + filename = fname.c_str(); + } + } + + if (!h->init(filename)) { + delete h; + return NULL; + } + *natoms = h->natoms; + return h; +} + +static int +read_timestep_metadata(void *v, molfile_timestep_metadata *m) { + FrameSetReader* handle = reinterpret_cast(v); + m->has_velocities = handle->has_velocities(); + m->count = handle->size(); + return MOLFILE_SUCCESS; +} + +static int read_next_timestep(void *v, int natoms, molfile_timestep_t *ts) { + FrameSetReader *h = reinterpret_cast(v); + return h->next(ts); +} + +#if defined(DESRES_READ_TIMESTEP2) +static int read_timestep2(void *v, ssize_t n, molfile_timestep_t *ts) { + FrameSetReader *h = reinterpret_cast(v); + return h->frame(n, ts); +} + +static molfile_ssize_t read_times(void *v, + molfile_ssize_t start, + molfile_ssize_t count, + double * times) { + FrameSetReader *h = reinterpret_cast(v); + return h->times(start, count, times); +} +#endif + +static void close_file_read( void *v ) { + FrameSetReader *h = reinterpret_cast(v); + delete h; +} + +static void *open_file_write(const char *path, const char *type, int natoms) { + DtrWriter *h = new DtrWriter(natoms); + if (!h->init(path)) { + delete h; + h=NULL; + } + return h; +} + +static int write_timestep(void *v, const molfile_timestep_t *ts) { + DtrWriter *h = reinterpret_cast(v); + return h->next(ts); +} + +static void close_file_write( void * v ) { + DtrWriter *h = reinterpret_cast(v); + delete h; +} + + +static molfile_plugin_t desmond; + +VMDPLUGIN_API int VMDPLUGIN_init (void) { + /* Plugin for desmond trajectory files */ + ::memset(&desmond,0,sizeof(desmond)); + desmond.abiversion = vmdplugin_ABIVERSION; + desmond.type = MOLFILE_PLUGIN_TYPE; + desmond.name = "dtr"; + desmond.prettyname = "DESRES Trajectory"; + desmond.author = "D.E. Shaw Research"; + desmond.majorv = 4; + desmond.minorv = 1; + desmond.is_reentrant = VMDPLUGIN_THREADUNSAFE; + + desmond.filename_extension = "dtr,dtr/,stk,atr,atr/"; + desmond.open_file_read = open_file_read; + desmond.read_timestep_metadata = read_timestep_metadata; + desmond.read_next_timestep = read_next_timestep; +#if defined(DESRES_READ_TIMESTEP2) + desmond.read_timestep2 = read_timestep2; + desmond.read_times = read_times; +#endif + desmond.close_file_read = close_file_read; + + desmond.open_file_write = open_file_write; + desmond.write_timestep = write_timestep; + desmond.close_file_write = close_file_write; + + return VMDPLUGIN_SUCCESS; +} + +VMDPLUGIN_API int VMDPLUGIN_register(void *v, vmdplugin_register_cb cb) { + cb(v,reinterpret_cast(&desmond)); + return VMDPLUGIN_SUCCESS; +} + +VMDPLUGIN_API int VMDPLUGIN_fini(void) { + return VMDPLUGIN_SUCCESS; +} +#endif /* ENABLE_VMD_PLUGIN */ + +#if defined(TEST_DTRPLUGIN) + +int main(int argc, char *argv[]) { + + /* check input arguments */ + if (argc != 3) { + fprintf(stderr, "Usage: %s input.dtr output.dtr\n", argv[0]); + return 1; + } + int natoms; + + /* read in all frames */ + void *handle = open_file_read( argv[1], "dtr", &natoms); + printf("got %d atoms\n", natoms); + molfile_timestep_t ts[1]; + std::vector timesteps; + do { + timesteps.resize( timesteps.size()+3*natoms ); + ts->coords = ×teps[ timesteps.size() - 3*natoms ]; + ts->velocities = NULL; + } while (read_next_timestep(handle, natoms, ts)==MOLFILE_SUCCESS); + int nframes = timesteps.size()/(3*natoms) - 1; + + printf("read %d steps\n", nframes ); + close_file_read(handle); + + /* write out all frames */ + handle = open_file_write( argv[2], "dtr", natoms); + if (!handle) return 1; + ts->coords = ×teps[0]; + for (int i=0; icoords += 3*natoms; + } + printf("wrote %d steps\n", nframes); + close_file_write(handle); + + /* now try to read it back in */ + int new_natoms; + handle = open_file_read( argv[2], "dtr", &new_natoms); + if (handle) return 1; + if (new_natoms != natoms) { + fprintf(stderr, "number of atoms changed: %d -> %d\n", natoms, new_natoms); + return 1; + } + close_file_read(handle); + return 0; +} + +#endif + +#if defined(TEST_DTR_DUMP) +int main(int argc, char *argv[]) { + + StkReader src, dst; + src.init(argv[1]); + + std::ostringstream out; + src.dump(out); + assert(out); + + std::istringstream in(out.str()); + dst.load(in); + assert(in); + + return 0; +} +#endif +#endif /* ENABLE_DTR */ diff --git a/src/vmdplugin/dtrplugin.hxx b/src/vmdplugin/dtrplugin.hxx new file mode 100644 index 0000000000..edfce34177 --- /dev/null +++ b/src/vmdplugin/dtrplugin.hxx @@ -0,0 +1,343 @@ +// +// Version info for VMD plugin tree: +// $Id: dtrplugin.hxx,v 1.5 2016/11/06 19:16:52 johns Exp $ +// +// Version info for last sync with D. E. Shaw Research: +// //depot/desrad/main/sw/libs/molfile/plugins/dtrplugin.hxx#13 +// + +/* +Copyright 2009, D. E. Shaw Research, LLC +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + +* Redistributions of source code must retain the above copyright + notice, this list of conditions, and the following disclaimer. + +* Redistributions in binary form must reproduce the above copyright + notice, this list of conditions, and the following disclaimer in the + documentation and/or other materials provided with the distribution. + +* Neither the name of D. E. Shaw Research, LLC nor the names of its + contributors may be used to endorse or promote products derived from + this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS +"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR +A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, +SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, +DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE +OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. + +*/ +#ifdef ENABLE_DTR +#ifndef MOLFILE_DTRPLUGIN_HXX +#define MOLFILE_DTRPLUGIN_HXX + +#if defined(_MSC_VER) +#ifndef DESRES_WIN32 +#define DESRES_WIN32 +#endif +#endif + +#ifndef __STDC_FORMAT_MACROS +#define __STDC_FORMAT_MACROS +#endif + +#include +#include +#ifdef DESRES_WIN32 +#include +#include +#include +#include + +typedef int int32_t; +typedef unsigned char uint8_t; +typedef unsigned int uint32_t; +#if 1 +typedef unsigned __int64 uint64_t; // This also works with MVSC6 +#else +typedef unsigned long long uint64_t; +#endif +typedef unsigned short uint16_t; +#if !defined(__MINGW32__) +typedef unsigned int ssize_t; +typedef int mode_t; +#endif +#define mkdir(a,b) _mkdir(a) +#define rmdir(a) _rmdir(a) +#define ftello(a) ftell(a) +#else +#define O_BINARY 0 +#include +#include +#include +#include +#endif + +#include "molfile_plugin.h" + +#include +#include +#include + + +namespace desres { namespace molfile { + + const char * dtr_serialized_version(); + + struct key_record_t { + uint32_t time_lo; /* Time associated with this file (low bytes). */ + uint32_t time_hi; /* Time associated with this file (high bytes). */ + uint32_t offset_lo; /* Zero in the 1 frame/file case (low bytes) */ + uint32_t offset_hi; /* Zero in the 1 frame/file case (high bytes) */ + uint32_t framesize_lo; /* Number of bytes in frame (low bytes) */ + uint32_t framesize_hi; /* Number of bytes in frame (high bytes) */ + + double time() const; + uint64_t offset() const; + uint64_t size() const; + + }; + + class Timekeys { + + double m_first; /* last time in the key list */ + double m_interval; /* representative time interval between keys */ + uint64_t m_framesize; /* size of all frames */ + size_t m_size; /* number of non-overlapping frames */ + size_t m_fullsize; /* total number of frames */ + uint32_t m_fpf; /* frames per file */ + + /* storage for keys until compressed, or if not compressible */ + std::vector keys; + + public: + Timekeys() + : m_first(0), m_interval(0), m_framesize(0), + m_size(0), m_fullsize(0), m_fpf(0) {} + + bool init( const std::string& path ); + + uint32_t framesperfile() const { return m_fpf; } + + bool is_compact() const { return keys.size()==0; } + + size_t size() const { return m_size; } + size_t full_size() const { return m_fullsize; } + + void truncate( size_t nframes ) { m_size = nframes; } + + void restore_full_size() { m_size = m_fullsize; } + + key_record_t operator[](uint64_t i) const; + + void dump( std::ostream& out ) const; + void load( std::istream& in ); + }; + + class DtrReader; + + class FrameSetReader { + protected: + std::string dtr; + + public: + uint32_t natoms; + bool with_velocity; + + FrameSetReader() + : natoms(0), with_velocity(false) + {} + + virtual ~FrameSetReader() {} + + bool has_velocities() const { return with_velocity; } + + const std::string &path() const { return dtr; } + + // initialize all members from frameset at given path. Return success. + // If changed is provided, set to true/false if timekeys were/were not + // reloaded. + virtual bool init(const std::string &path, int * changed = NULL) = 0; + + // number of frames + virtual ssize_t size() const = 0; + + // read the next frame. If ts is NULL, skip it. If there are no more + // frames to read, return MOLFILE_EOF. Otherwise, fill in coordinates + // in ts->coords, velocities in ts->velocities if ts->velocities is + // non-NULL, and return MOLFILE_SUCCESS if all goes well. + virtual int next(molfile_timestep_t *ts) = 0; + + // Get the DtrReader component corresponding to frame n. Change + // n to the local index within the returned dtr. + virtual const DtrReader * component(ssize_t &n) const = 0; + + // number of framesets + virtual ssize_t nframesets() const = 0; + + // nth frameset + virtual const DtrReader * frameset(ssize_t n) const = 0; + + // read a specified frame. + virtual int frame(ssize_t n, molfile_timestep_t *ts) const = 0; + + // read up to count times beginning at index start into the provided space; + // return the number of times actually read. + virtual ssize_t times(ssize_t start, ssize_t count, double * times) const = 0; + + virtual std::ostream& dump(std::ostream &out) const = 0; + virtual std::istream& load(std::istream &in) = 0; + }; + + struct metadata_t { + std::vector invmass; + }; + + class DtrReader : public FrameSetReader { + mutable int m_ndir1; + mutable int m_ndir2; + ssize_t m_curframe; + + metadata_t * meta; + bool owns_meta; + + bool eof() const { return m_curframe >= (ssize_t)keys.size(); } + + public: + DtrReader() + : m_ndir1(-1), m_ndir2(-1), m_curframe(0), + meta(NULL), owns_meta(false) {} + + virtual ~DtrReader() { + set_meta(NULL); + } + + Timekeys keys; + + // lazy-loaded DDparams + int ndir1() const; + int ndir2() const; + + /* meta and owns_meta are initially set to NULL and false. In init(), + * if meta is NULL and owns_meta is false, we try to read the meta + * from the metadata frame (an expensive operation), in which case + * owns_meta becomes true, whether or not meta was actually read. + * Otherwise, we leave those values alone. In the destructor, we delete + * meta if we own it. The StkReader class can share the meta between + * DtrReader instances in the following way: if the meta it has to share + * is NULL, it should set owns_meta to true in the DtrReader instances, + * so that the DtrReaders don't keep searching for their own meta. If + * the meta it has to share is non-NULL, it should set owns_meta to + * false so that the meta pointer doesn't get double-freed. + */ + metadata_t * get_meta() const { return meta; } + void set_meta(metadata_t * ptr) { + if (meta && owns_meta) delete meta; + if (ptr) { + meta = ptr; + owns_meta = false; + } else { + meta = NULL; + owns_meta = true; + } + } + + ssize_t curframe() const { return m_curframe; } + + uint32_t framesperfile() const { return keys.framesperfile(); } + + virtual bool init(const std::string &path, int * changed=NULL); + virtual ssize_t size() const { return keys.size(); } + virtual ssize_t times(ssize_t start, ssize_t count, double * times) const; + + virtual int next(molfile_timestep_t *ts); + + virtual const DtrReader * component(ssize_t &n) const { + return this; + } + + virtual ssize_t nframesets() const { return 1; } + virtual const DtrReader * frameset(ssize_t n) const { + if (n!=0) throw std::runtime_error("bad index"); + return this; + } + + virtual int frame(ssize_t n, molfile_timestep_t *ts) const; + + + // path for frame at index. Empty string on not found. + std::string framefile(ssize_t n) const; + + // read a frame from supplied bytes + int frame_from_bytes( const void *buf, uint64_t len, + molfile_timestep_t *ts ) const; + + std::ostream& dump(std::ostream &out) const; + std::istream& load(std::istream &in); + }; + + struct DtrWriter { + std::string dtr; + std::string m_directory; + const uint32_t natoms; + int frame_fd; // for writing + uint32_t frames_per_file; + uint64_t framefile_offset; + uint64_t nwritten; + double last_time; + FILE * timekeys_file; + + explicit DtrWriter(uint32_t natoms_) + : natoms(natoms_), frame_fd(0), frames_per_file(256), framefile_offset(0), + nwritten(0), last_time(HUGE_VAL), timekeys_file(NULL) + {} + + ~DtrWriter(); + + // initialize for writing at path + bool init(const std::string &path); + + // write another frame. Return MOLFILE_SUCCESS or MOLFILE_ERROR. + int next(const molfile_timestep_t *ts); + }; + + class StkReader : public FrameSetReader { + std::vector framesets; + size_t curframeset; + + public: + StkReader() : curframeset(0) {} + StkReader(DtrReader *reader); + ~StkReader(); + virtual bool init(const std::string &path, int * changed=NULL); + virtual ssize_t size() const; + virtual ssize_t times(ssize_t start, ssize_t count, double * times) const; + virtual int next(molfile_timestep_t *ts); + virtual int frame(ssize_t n, molfile_timestep_t *ts) const; + virtual const DtrReader * component(ssize_t &n) const; + + virtual ssize_t nframesets() const { return framesets.size(); } + virtual const DtrReader * frameset(ssize_t n) const { + return framesets.at(n); + } + + static bool recognizes(const std::string &path); + + std::ostream& dump(std::ostream &out) const; + std::istream& load(std::istream &in); + }; +} } + +#endif +#endif /* ENABLE_DTR */ diff --git a/src/vmdplugin/molfile_plugin.h b/src/vmdplugin/molfile_plugin.h new file mode 100644 index 0000000000..b1cb0a5593 --- /dev/null +++ b/src/vmdplugin/molfile_plugin.h @@ -0,0 +1,990 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2006 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ + +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: molfile_plugin.h,v $ + * $Author: johns $ $Locker: $ $State: Exp $ + * $Revision: 1.112 $ $Date: 2019/10/17 06:12:24 $ + * + ***************************************************************************/ + +/** @file + * API for C extensions to define a way to load structure, coordinate, + * trajectory, and volumetric data files + */ + +#ifndef MOL_FILE_PLUGIN_H +#define MOL_FILE_PLUGIN_H + +#include "vmdplugin.h" + +#if defined(DESRES_READ_TIMESTEP2) +/* includes needed for large integer types used for frame counts */ +#include +typedef ssize_t molfile_ssize_t; /**< for frame counts */ +#endif + +/** + * Define a common plugin type to be used when registering the plugin. + */ +#define MOLFILE_PLUGIN_TYPE "mol file reader" + +/** + * File converter plugins use the same API but register under a different + * type so that regular file readers can have priority. + */ +#define MOLFILE_CONVERTER_PLUGIN_TYPE "mol file converter" + +/* File plugin symbolic constants for better code readability */ +#define MOLFILE_SUCCESS 0 /**< succeeded in reading file */ +#define MOLFILE_EOF -1 /**< end of file */ +#define MOLFILE_ERROR -1 /**< error reading/opening a file */ +#define MOLFILE_NOSTRUCTUREDATA -2 /**< no structure data in this file */ + +#define MOLFILE_NUMATOMS_UNKNOWN -1 /**< unknown number of atoms */ +#define MOLFILE_NUMATOMS_NONE 0 /**< no atoms in this file type */ + +/** + * Maximum string size macro + */ +#define MOLFILE_BUFSIZ 81 /**< maximum chars in string data */ +#define MOLFILE_BIGBUFSIZ 4096 /**< maximum chars in long strings */ + +#define MOLFILE_MAXWAVEPERTS 25 /**< maximum number of wavefunctions + * per timestep */ + +/** + * Hard-coded direct-I/O page size constants for use by both VMD + * and the plugins that want to use direct, unbuffered I/O for high + * performance with SSDs etc. We use two constants to define the + * range of hardware page sizes that we can support, so that we can + * add support for larger 8KB or 16KB page sizes in the future + * as they become more prevalent in high-end storage systems. + * + * At present, VMD uses a hard-coded 4KB page size to reduce memory + * fragmentation, but these constants will make it easier to enable the + * use of larger page sizes in the future if it becomes necessary. + */ +#define MOLFILE_DIRECTIO_MIN_BLOCK_SIZE 4096 +#define MOLFILE_DIRECTIO_MAX_BLOCK_SIZE 4096 + + +/** + * File level comments, origin information, and annotations. + */ +typedef struct { + char database[81]; /**< database of origin, if any */ + char accession[81]; /**< database accession code, if any */ + char date[81]; /**< date/time stamp for this data */ + char title[81]; /**< brief title for this data */ + int remarklen; /**< length of remarks string */ + char *remarks; /**< free-form remarks about data */ +} molfile_metadata_t; + + +/* + * Struct for specifying atoms in a molecular structure. The first + * six components are required, the rest are optional and their presence is + * indicating by setting the corresponding bit in optsflag. When omitted, + * the application (for read_structure) or plugin (for write_structure) + * must be able to supply default values if the missing parameters are + * part of its internal data structure. + * Note that it is not possible to specify coordinates with this structure. + * This is intentional; all coordinate I/O is done with the read_timestep and + * write_timestep functions. + */ + +/** + * Per-atom attributes and information. + */ +typedef struct { + /* these fields absolutely must be set or initialized to empty */ + char name[16]; /**< required atom name string */ + char type[16]; /**< required atom type string */ + char resname[8]; /**< required residue name string */ + int resid; /**< required integer residue ID */ + char segid[8]; /**< required segment name string, or "" */ +#if 0 && vmdplugin_ABIVERSION > 10000 + /* The new PDB file formats allows for much larger structures, */ + /* which can therefore require longer chain ID strings. The */ + /* new PDBx/mmCIF file formats do not have length limits on */ + /* fields, so PDB chains could be arbitrarily long strings */ + /* in such files. At present, we know we need at least 3-char */ + /* chains for existing PDBx/mmCIF files. */ + char chain[4]; /**< required chain name, or "" */ +#else + char chain[2]; /**< required chain name, or "" */ +#endif + /* rest are optional; use optflags to specify what's present */ + char altloc[2]; /**< optional PDB alternate location code */ + char insertion[2]; /**< optional PDB insertion code */ + float occupancy; /**< optional occupancy value */ + float bfactor; /**< optional B-factor value */ + float mass; /**< optional mass value */ + float charge; /**< optional charge value */ + float radius; /**< optional radius value */ + int atomicnumber; /**< optional element atomic number */ + +#if 0 + char complex[16]; + char assembly[16]; + int qmregion; + int qmregionlink; + int qmlayer; + int qmlayerlink; + int qmfrag; + int qmfraglink; + string qmecp; + int qmadapt; + int qmect; /**< boolean */ + int qmparam; + int autoparam; +#endif + +#if defined(DESRES_CTNUMBER) + int ctnumber; /**< mae ct block, 0-based, including meta */ +#endif +} molfile_atom_t; + +/*@{*/ +/** Plugin optional data field availability flag */ +#define MOLFILE_NOOPTIONS 0x0000 /**< no optional data */ +#define MOLFILE_INSERTION 0x0001 /**< insertion codes provided */ +#define MOLFILE_OCCUPANCY 0x0002 /**< occupancy data provided */ +#define MOLFILE_BFACTOR 0x0004 /**< B-factor data provided */ +#define MOLFILE_MASS 0x0008 /**< Atomic mass provided */ +#define MOLFILE_CHARGE 0x0010 /**< Atomic charge provided */ +#define MOLFILE_RADIUS 0x0020 /**< Atomic VDW radius provided */ +#define MOLFILE_ALTLOC 0x0040 /**< Multiple conformations present */ +#define MOLFILE_ATOMICNUMBER 0x0080 /**< Atomic element number provided */ +#define MOLFILE_BONDSSPECIAL 0x0100 /**< Only non-standard bonds provided */ +#if defined(DESRES_CTNUMBER) +#define MOLFILE_CTNUMBER 0x0200 /**< ctnumber provided */ +#endif +#define MOLFILE_BADOPTIONS 0xFFFFFFFF /**< Detect badly behaved plugins */ + +/*@}*/ + +/*@{*/ +/** Flags indicating availability of optional data fields + * for QM timesteps + */ +#define MOLFILE_QMTS_NOOPTIONS 0x0000 /**< no optional data */ +#define MOLFILE_QMTS_GRADIENT 0x0001 /**< energy gradients provided */ +#define MOLFILE_QMTS_SCFITER 0x0002 +/*@}*/ + +typedef struct molfile_timestep_metadata { + unsigned int count; /**< total # timesteps; -1 if unknown */ + unsigned int avg_bytes_per_timestep; /** bytes per timestep */ + int has_velocities; /**< if timesteps have velocities */ +} molfile_timestep_metadata_t; + + +/* + * + * Per-timestep atom coordinates, velocities, forces, energies, + * and periodic cell information + * + */ + +#if 0 +/** + * Periodically stored energies of various kinds + */ +typedef struct { + int energyflags; // XXX indicate use and semantics of data fields + + double total_energy; // XXX these copied from DESRES_READ_TIMESTEP2 case + double potential_energy; + double kinetic_energy; + double extended_energy; + double force_energy; + double total_pressure; + + // Alchemical free energy methods need to store individual energy samples. + // We don't really want pre-averaged quantities, they lead to problems later. + double lambda; // data gen sim parm: + aux scheduling + soft core parm + double temperature; // temp set by thermostat + // either we use deltaU, or we store U and Uprime... + double deltaU; // Ulambda - Ulambdaprime + double Ulambda; // U for lambda + double Ulambdaprime; // U for lambdaprime + + // IDWS methods + // XXX both values of lambdaprime and required de-interleaving information + + // Replica exchange methods + // XXX rectangle sample params? + + // REST2 method + // XXX + +} molfile_energies_t; +#endif + + +/** + * Per-timestep atom coordinates, velocities, time, energies + * and periodic cell info + */ +typedef struct { + float *coords; /**< coordinates of all atoms, arranged xyzxyzxyz */ + float *velocities; /**< space for velocities of all atoms; same layout */ + /**< NULL unless has_velocities is set */ + + /*@{*/ + /** + * Unit cell specification of the form A, B, C, alpha, beta, gamma. + * notes: A, B, C are side lengths of the unit cell + * alpha = angle between b and c + * beta = angle between a and c + * gamma = angle between a and b + */ + float A, B, C, alpha, beta, gamma; + /*@}*/ + + double physical_time; /**< physical time point associated with this frame */ + +#if defined(DESRES_READ_TIMESTEP2) + /* HACK to support generic trajectory information */ + double total_energy; + double potential_energy; + double kinetic_energy; + double extended_energy; + double force_energy; + double total_pressure; +#endif + +} molfile_timestep_t; + + +/** + * Metadata for volumetric datasets, read initially and used for subsequent + * memory allocations and file loading. + */ +typedef struct { + char dataname[256]; /**< name of volumetric data set */ + float origin[3]; /**< origin: origin of volume (x=0, y=0, z=0 corner */ + + /* + * x/y/z axis: + * These the three cell sides, providing both direction and length + * (not unit vectors) for the x, y, and z axes. In the simplest + * case, these would be <0,size,0> and <0,0,size) for + * an orthogonal cubic volume set. For other cell shapes these + * axes can be oriented non-orthogonally, and the parallelpiped + * may have different side lengths, not just a cube/rhombus. + */ + float xaxis[3]; /**< direction (and length) for X axis */ + float yaxis[3]; /**< direction (and length) for Y axis */ + float zaxis[3]; /**< direction (and length) for Z axis */ + + /* + * x/y/z size: + * Number of grid cells along each axis. This is _not_ the + * physical size of the box, this is the number of voxels in each + * direction, independent of the shape of the volume set. + */ + int xsize; /**< number of grid cells along the X axis */ + int ysize; /**< number of grid cells along the Y axis */ + int zsize; /**< number of grid cells along the Z axis */ + +#if vmdplugin_ABIVERSION > 16 + int has_scalar; /**< flag indicating presence of scalar volume */ + int has_gradient; /**< flag indicating presence of vector volume */ + int has_variance; /**< flag indicating presence of variance map */ +#endif + int has_color; /**< flag indicating presence of voxel color data */ +} molfile_volumetric_t; + + +#if vmdplugin_ABIVERSION > 16 +/** + * Volumetric dataset read/write structure with both flag/parameter sets + * and VMD-allocated pointers for fields to be used by the plugin. + */ +typedef struct { + int setidx; /**< volumetric dataset index to load/save */ + float *scalar; /**< scalar density/potential field data */ + float *gradient; /**< gradient vector field */ + float *variance; /**< variance map indicating signal/noise */ + float *rgb3f; /**< RGB floating point color texture map */ + unsigned char *rgb3u; /**< RGB unsigned byte color texture map */ +} molfile_volumetric_readwrite_t; +#endif + + +/************************************************************** + ************************************************************** + **** **** + **** Data structures for QM files **** + **** **** + ************************************************************** + **************************************************************/ + +/* macros for the convergence status of a QM calculation. */ +#define MOLFILE_QMSTATUS_UNKNOWN -1 /* don't know yet */ +#define MOLFILE_QMSTATUS_OPT_CONV 0 /* optimization converged */ +#define MOLFILE_QMSTATUS_SCF_NOT_CONV 1 /* SCF convergence failed */ +#define MOLFILE_QMSTATUS_OPT_NOT_CONV 2 /* optimization not converged */ +#define MOLFILE_QMSTATUS_FILE_TRUNCATED 3 /* file was truncated */ + +/* macros describing the SCF method (SCFTYP in GAMESS) */ +#define MOLFILE_SCFTYPE_UNKNOWN -1 /* no info about the method */ +#define MOLFILE_SCFTYPE_NONE 0 /* calculation didn't make use of SCF */ +#define MOLFILE_SCFTYPE_RHF 1 /* restricted Hartree-Fock */ +#define MOLFILE_SCFTYPE_UHF 2 /* unrestricted Hartree-Fock */ +#define MOLFILE_SCFTYPE_ROHF 3 /* restricted open-shell Hartree-Fock */ +#define MOLFILE_SCFTYPE_GVB 4 /* generalized valence bond orbitals */ +#define MOLFILE_SCFTYPE_MCSCF 5 /* multi-configuration SCF */ +#define MOLFILE_SCFTYPE_FF 6 /* classical force-field based sim. */ + +/* macros describing the type of calculation (RUNTYP in GAMESS) */ +#define MOLFILE_RUNTYPE_UNKNOWN 0 /* single point run */ +#define MOLFILE_RUNTYPE_ENERGY 1 /* single point run */ +#define MOLFILE_RUNTYPE_OPTIMIZE 2 /* geometry optimization */ +#define MOLFILE_RUNTYPE_SADPOINT 3 /* saddle point search */ +#define MOLFILE_RUNTYPE_HESSIAN 4 /* Hessian/frequency calculation */ +#define MOLFILE_RUNTYPE_SURFACE 5 /* potential surface scan */ +#define MOLFILE_RUNTYPE_GRADIENT 6 /* energy gradient calculation */ +#define MOLFILE_RUNTYPE_MEX 7 /* minimum energy crossing */ +#define MOLFILE_RUNTYPE_DYNAMICS 8 /* Any type of molecular dynamics + * e.g. Born-Oppenheimer, Car-Parinello, + * or classical MD */ +#define MOLFILE_RUNTYPE_PROPERTIES 9 /* Properties were calculated from a + * wavefunction that was read from file */ + + +/** + * Sizes of various QM-related, timestep independent data arrays + * which must be allocated by the caller (VMD) so that the plugin + * can fill in the arrays with data. + */ +typedef struct { + /* hessian data */ + int nimag; /**< number of imaginary modes */ + int nintcoords; /**< number internal coordinates */ + int ncart; /**< number cartesian coordinates */ + + /* orbital/basisset data */ + int num_basis_funcs; /**< number of uncontracted basis functions in basis array */ + int num_basis_atoms; /**< number of atoms in basis set */ + int num_shells; /**< total number of atomic shells */ + int wavef_size; /**< size of the wavefunction + * i.e. size of secular eq. or + * # of cartesian contracted + * gaussian basis functions */ + + /* everything else */ + int have_sysinfo; + int have_carthessian; /**< hessian in cartesian coords available */ + int have_inthessian; /**< hessian in internal coords available */ + int have_normalmodes; /**< normal modes available */ +} molfile_qm_metadata_t; + + +/** + * QM run info. Parameters that stay unchanged during a single file. + */ +typedef struct { + int nproc; /**< number of processors used. */ + int memory; /**< amount of memory used in Mbyte. */ + int runtype; /**< flag indicating the calculation method. */ + int scftype; /**< SCF type: RHF, UHF, ROHF, GVB or MCSCF wfn. */ + int status; /**< indicates wether SCF and geometry optimization + * have converged properly. */ + int num_electrons; /**< number of electrons. XXX: can be fractional in some DFT codes */ + int totalcharge; /**< total charge of system. XXX: can be fractional in some DFT codes */ + int num_occupied_A; /**< number of occupied alpha orbitals */ + int num_occupied_B; /**< number of occupied beta orbitals */ + + double *nuc_charge; /**< array(natom) containing the nuclear charge of atom i */ + + char basis_string[MOLFILE_BUFSIZ]; /**< basis name as "nice" string. */ + char runtitle[MOLFILE_BIGBUFSIZ]; /**< title of run. */ + char geometry[MOLFILE_BUFSIZ]; /**< type of provided geometry, XXX: remove? + * e.g. UNIQUE, ZMT, CART, ... */ + char version_string[MOLFILE_BUFSIZ]; /**< QM code version information. */ +} molfile_qm_sysinfo_t; + + +/** + * Data for QM basis set + */ +typedef struct { + int *num_shells_per_atom; /**< number of shells per atom */ + int *num_prim_per_shell; /**< number of shell primitives shell */ + + float *basis; /**< contraction coeffients and exponents for + * the basis functions in the form + * {exp(1), c-coeff(1), exp(2), c-coeff(2), ...}; + * array size = 2*num_basis_funcs + * The basis must NOT be normalized. */ + int *atomic_number; /**< atomic numbers (chem. element) of atoms in basis set */ + int *angular_momentum; /**< 3 ints per wave function coefficient do describe the + * cartesian components of the angular momentum. + * E.g. S={0 0 0}, Px={1 0 0}, Dxy={1 1 0}, or Fyyz={0 2 1}. + */ + int *shell_types; /**< type for each shell in basis */ +} molfile_qm_basis_t; + + +/** + * Data from QM Hessian/normal mode runs + * + * A noteworthy comment from one of Axel's emails: + * The molfile_qm_hessian_t, I'd rename to molfile_hessian_t (one + * can do vibrational analysis without QM) and would make this a + * completely separate entity. This could then be also used to + * read in data from, say, principal component analysis or normal + * mode analysis and VMD could contain code to either project a + * trajectory on the contained eigenvectors or animate them and + * so on. There is a bunch of possible applications... + */ +typedef struct { + double *carthessian; /**< hessian matrix in cartesian coordinates (ncart)*(ncart) + * as a single array of doubles (row(1), ...,row(natoms)) */ + int *imag_modes; /**< list(nimag) of imaginary modes */ + double *inthessian; /**< hessian matrix in internal coordinates + * (nintcoords*nintcoords) as a single array of + * doubles (row(1), ...,row(nintcoords)) */ + float *wavenumbers; /**< array(ncart) of wavenumbers of normal modes */ + float *intensities; /**< array(ncart) of intensities of normal modes */ + float *normalmodes; /**< matrix(ncart*ncart) of normal modes */ +} molfile_qm_hessian_t; + + +/** + * QM related information that is timestep independent + */ +typedef struct { + molfile_qm_sysinfo_t run; /* system info */ + molfile_qm_basis_t basis; /* basis set info */ + molfile_qm_hessian_t hess; /* hessian info */ +} molfile_qm_t; + + + +/** + * Enumeration of all of the wavefunction types that can be read + * from QM file reader plugins. + * + * CANON = canonical (i.e diagonalized) wavefunction + * GEMINAL = GVB-ROHF geminal pairs + * MCSCFNAT = Multi-Configuration SCF natural orbitals + * MCSCFOPT = Multi-Configuration SCF optimized orbitals + * CINATUR = Configuration-Interaction natural orbitals + * BOYS = Boys localization + * RUEDEN = Ruedenberg localization + * PIPEK = Pipek-Mezey population localization + * + * NBO related localizations: + * -------------------------- + * NAO = Natural Atomic Orbitals + * PNAO = pre-orthogonal NAOs + * NBO = Natural Bond Orbitals + * PNBO = pre-orthogonal NBOs + * NHO = Natural Hybrid Orbitals + * PNHO = pre-orthogonal NHOs + * NLMO = Natural Localized Molecular Orbitals + * PNLMO = pre-orthogonal NLMOs + * + * UNKNOWN = Use this for any type not listed here + * You can use the string field for description + */ +enum molfile_qm_wavefunc_type { + MOLFILE_WAVE_CANON, MOLFILE_WAVE_GEMINAL, + MOLFILE_WAVE_MCSCFNAT, MOLFILE_WAVE_MCSCFOPT, + MOLFILE_WAVE_CINATUR, + MOLFILE_WAVE_PIPEK, MOLFILE_WAVE_BOYS, MOLFILE_WAVE_RUEDEN, + MOLFILE_WAVE_NAO, MOLFILE_WAVE_PNAO, MOLFILE_WAVE_NHO, + MOLFILE_WAVE_PNHO, MOLFILE_WAVE_NBO, MOLFILE_WAVE_PNBO, + MOLFILE_WAVE_PNLMO, MOLFILE_WAVE_NLMO, MOLFILE_WAVE_MOAO, + MOLFILE_WAVE_NATO, MOLFILE_WAVE_UNKNOWN +}; + + +/** + * Enumeration of all of the supported QM related charge + * types + */ +enum molfile_qm_charge_type { + MOLFILE_QMCHARGE_UNKNOWN, + MOLFILE_QMCHARGE_MULLIKEN, MOLFILE_QMCHARGE_LOWDIN, + MOLFILE_QMCHARGE_ESP, MOLFILE_QMCHARGE_NPA +}; + + + +/** + * Sizes of various QM-related, per-timestep data arrays + * which must be allocated by the caller (VMD) so that the plugin + * can fill in the arrays with data. + */ +typedef struct molfile_qm_timestep_metadata { + unsigned int count; /**< total # timesteps; -1 if unknown */ + unsigned int avg_bytes_per_timestep; /**< bytes per timestep */ + int has_gradient; /**< if timestep contains gradient */ + int num_scfiter; /**< # scf iterations for this ts */ + int num_orbitals_per_wavef[MOLFILE_MAXWAVEPERTS]; /**< # orbitals for each wavefunction */ + int has_orben_per_wavef[MOLFILE_MAXWAVEPERTS]; /**< orbital energy flags */ + int has_occup_per_wavef[MOLFILE_MAXWAVEPERTS]; /**< orbital occupancy flags */ + int num_wavef ; /**< # wavefunctions in this ts */ + int wavef_size; /**< size of one wavefunction + * (# of gaussian basis fctns) */ + int num_charge_sets; /**< # of charge values per atom */ +} molfile_qm_timestep_metadata_t; + + +/** + * QM wavefunction + */ +typedef struct { + int type; /**< MOLFILE_WAVE_CANON, MOLFILE_WAVE_MCSCFNAT, ... */ + int spin; /**< 1 for alpha, -1 for beta */ + int excitation; /**< 0 for ground state, 1,2,3,... for excited states */ + int multiplicity; /**< spin multiplicity of the state, zero if unknown */ + char info[MOLFILE_BUFSIZ]; /**< string for additional type info */ + + double energy; /**< energy of the electronic state. + * i.e. HF-SCF energy, CI state energy, + * MCSCF energy, etc. */ + + float *wave_coeffs; /**< expansion coefficients for wavefunction in the + * form {orbital1(c1),orbital1(c2),.....,orbitalM(cN)} */ + float *orbital_energies; /**< list of orbital energies for wavefunction */ + float *occupancies; /**< orbital occupancies */ + int *orbital_ids; /**< orbital ID numbers; If NULL then VMD will + * assume 1,2,3,...num_orbs. */ +} molfile_qm_wavefunction_t; + + +/** + * QM per trajectory timestep info + * Note that each timestep can contain multiple wavefunctions. + */ +typedef struct { + molfile_qm_wavefunction_t *wave; /**< array of wavefunction objects */ + float *gradient; /**< force on each atom (=gradient of energy) */ + + double *scfenergies; /**< energies from the SCF cycles */ + double *charges; /**< per-atom charges */ + int *charge_types; /**< type of each charge set */ +} molfile_qm_timestep_t; + + +/************************************************************** + **************************************************************/ + + + + +/** + * Enumeration of all of the supported graphics objects that can be read + * from graphics file reader plugins. + */ +enum molfile_graphics_type { + MOLFILE_POINT, MOLFILE_TRIANGLE, MOLFILE_TRINORM, MOLFILE_NORMS, + MOLFILE_LINE, MOLFILE_CYLINDER, MOLFILE_CAPCYL, MOLFILE_CONE, + MOLFILE_SPHERE, MOLFILE_TEXT, MOLFILE_COLOR, MOLFILE_TRICOLOR +}; + +/** + * Individual graphics object/element data + */ +typedef struct { + int type; /* One of molfile_graphics_type */ + int style; /* A general style parameter */ + float size; /* A general size parameter */ + float data[9]; /* All data for the element */ +} molfile_graphics_t; + + +/* + * Types for raw graphics elements stored in files. Data for each type + * should be stored by the plugin as follows: + +type data style size +---- ---- ----- ---- +point x, y, z pixel size +triangle x1,y1,z1,x2,y2,z2,x3,y3,z3 +trinorm x1,y1,z1,x2,y2,z2,x3,y3,z3 + the next array element must be NORMS +tricolor x1,y1,z1,x2,y2,z2,x3,y3,z3 + the next array elements must be NORMS + the following element must be COLOR, with three RGB triples +norms x1,y1,z1,x2,y2,z2,x3,y3,z3 +line x1,y1,z1,x2,y2,z2 0=solid pixel width + 1=stippled +cylinder x1,y1,z1,x2,y2,z2 resolution radius +capcyl x1,y1,z1,x2,y2,z2 resolution radius +sphere x1,y1,z1 resolution radius +text x, y, z, up to 24 bytes of text pixel size +color r, g, b +*/ + + +/** + * Main file reader API. Any function in this struct may be NULL + * if not implemented by the plugin; the application checks this to determine + * what functionality is present in the plugin. + */ +typedef struct { + /** + * Required header + */ + vmdplugin_HEAD + + /** + * Filename extension for this file type. May be NULL if no filename + * extension exists and/or is known. For file types that match several + * common extensions, list them in a comma separated list such as: + * "pdb,ent,foo,bar,baz,ban" + * The comma separated list will be expanded when filename extension matching + * is performed. If multiple plugins solicit the same filename extensions, + * the one that lists the extension earliest in its list is selected. In the + * case of a "tie", the first one tried/checked "wins". + */ + const char *filename_extension; + + /** + * Try to open the file for reading. Return an opaque handle, or NULL on + * failure. Set the number of atoms; if the number of atoms cannot be + * determined, set natoms to MOLFILE_NUMATOMS_UNKNOWN. + * Filetype should be the name under which this plugin was registered; + * this is provided so that plugins can provide the same function pointer + * to handle multiple file types. + */ + void *(* open_file_read)(const char *filepath, const char *filetype, + int *natoms); + + /** + * Read molecular structure from the given file handle. atoms is allocated + * by the caller and points to space for natoms. + * On success, place atom information in the passed-in pointer. + * optflags specifies which optional fields in the atoms will be set by + * the plugin. + */ + int (*read_structure)(void *, int *optflags, molfile_atom_t *atoms); + + /** + * Read bond information for the molecule. On success the arrays from + * and to should point to the (one-based) indices of bonded atoms. + * Each unique bond should be specified only once, so file formats that list + * bonds twice will need post-processing before the results are returned to + * the caller. + * If the plugin provides bond information, but the file loaded doesn't + * actually contain any bond info, the nbonds parameter should be + * set to 0 and from/to should be set to NULL to indicate that no bond + * information was actually present, and automatic bond search should be + * performed. + * + * If the plugin provides bond order information, the bondorder array + * will contain the bond order for each from/to pair. If not, the bondorder + * pointer should be set to NULL, in which case the caller will provide a + * default bond order value of 1.0. + * + * If the plugin provides bond type information, the bondtype array + * will contain the bond type index for each from/to pair. These numbers + * are consecutive integers starting from 0. + * the bondtypenames list, contains the corresponding names, if available, + * as a NULL string terminated list. nbondtypes is provided for convenience + * and consistency checking. + * + * These arrays must be freed by the plugin in the close_file_read function. + * This function can be called only after read_structure(). + * Return MOLFILE_SUCCESS if no errors occur. + */ + int (*read_bonds)(void *, int *nbonds, int **from, int **to, float **bondorder, + int **bondtype, int *nbondtypes, char ***bondtypename); + + /** + * XXX this function will be augmented and possibly superceded by a + * new QM-capable version named read_timestep(), when finished. + * + * Read the next timestep from the file. Return MOLFILE_SUCCESS, or + * MOLFILE_EOF on EOF. If the molfile_timestep_t argument is NULL, then + * the frame should be skipped. Otherwise, the application must prepare + * molfile_timestep_t by allocating space in coords for the corresponding + * number of coordinates. + * The natoms parameter exists because some coordinate file formats + * (like CRD) cannot determine for themselves how many atoms are in a + * timestep; the app must therefore obtain this information elsewhere + * and provide it to the plugin. + */ + int (* read_next_timestep)(void *, int natoms, molfile_timestep_t *); + + /** + * Close the file and release all data. The handle cannot be reused. + */ + void (* close_file_read)(void *); + + /** + * Open a coordinate file for writing using the given header information. + * Return an opaque handle, or NULL on failure. The application must + * specify the number of atoms to be written. + * filetype should be the name under which this plugin was registered. + */ + void *(* open_file_write)(const char *filepath, const char *filetype, + int natoms); + + /** + * Write structure information. Return success. + */ + int (* write_structure)(void *, int optflags, const molfile_atom_t *atoms); + + /** + * Write a timestep to the coordinate file. Return MOLFILE_SUCCESS if no + * errors occur. If the file contains structure information in each + * timestep (like a multi-entry PDB), it will have to cache the information + * from the initial calls from write_structure. + */ + int (* write_timestep)(void *, const molfile_timestep_t *); + + /** + * Close the file and release all data. The handle cannot be reused. + */ + void (* close_file_write)(void *); + + /** + * Retrieve metadata pertaining to volumetric datasets in this file. + * Set nsets to the number of volumetric data sets, and set *metadata + * to point to an array of molfile_volumetric_t. The array is owned by + * the plugin and should be freed by close_file_read(). The application + * may call this function any number of times. + */ + int (* read_volumetric_metadata)(void *, int *nsets, + molfile_volumetric_t **metadata); + + /** + * Read the specified volumetric data set into the space pointed to by + * datablock. The set is specified with a zero-based index. The space + * allocated for the datablock must be equal to + * xsize * ysize * zsize. No space will be allocated for colorblock + * unless has_color is nonzero; in that case, colorblock should be + * filled in with three RGB floats per datapoint. + */ + int (* read_volumetric_data)(void *, int set, float *datablock, + float *colorblock); +#if vmdplugin_ABIVERSION > 16 + int (* read_volumetric_data_ex)(void *, molfile_volumetric_readwrite_t *v); +#endif + + /** + * Read raw graphics data stored in this file. Return the number of data + * elements and the data itself as an array of molfile_graphics_t in the + * pointer provided by the application. The plugin is responsible for + * freeing the data when the file is closed. + */ + int (* read_rawgraphics)(void *, int *nelem, const molfile_graphics_t **data); + + /** + * Read molecule metadata such as what database (if any) this file/data + * came from, what the accession code for the database is, textual remarks + * and other notes pertaining to the contained structure/trajectory/volume + * and anything else that's informative at the whole file level. + */ + int (* read_molecule_metadata)(void *, molfile_metadata_t **metadata); + + /** + * Write bond information for the molecule. The arrays from + * and to point to the (one-based) indices of bonded atoms. + * Each unique bond will be specified only once by the caller. + * File formats that list bonds twice will need to emit both the + * from/to and to/from versions of each. + * This function must be called before write_structure(). + * + * Like the read_bonds() routine, the bondorder pointer is set to NULL + * if the caller doesn't have such information, in which case the + * plugin should assume a bond order of 1.0 if the file format requires + * bond order information. + * + * Support for bond types follows the bondorder rules. bondtype is + * an integer array of the size nbonds that contains the bond type + * index (consecutive integers starting from 0) and bondtypenames + * contain the corresponding strings, in case the naming/numbering + * scheme is different from the index numbers. + * if the pointers are set to NULL, then this information is not available. + * bondtypenames can only be used of bondtypes is also given. + * Return MOLFILE_SUCCESS if no errors occur. + */ + int (* write_bonds)(void *, int nbonds, int *from, int *to, float *bondorder, + int *bondtype, int nbondtypes, char **bondtypename); + + /** + * Write the specified volumetric data set into the space pointed to by + * datablock. The * allocated for the datablock must be equal to + * xsize * ysize * zsize. No space will be allocated for colorblock + * unless has_color is nonzero; in that case, colorblock should be + * filled in with three RGB floats per datapoint. + */ + int (* write_volumetric_data)(void *, molfile_volumetric_t *metadata, + float *datablock, float *colorblock); +#if vmdplugin_ABIVERSION > 16 + int (* write_volumetric_data_ex)(void *, molfile_volumetric_t *metadata, + molfile_volumetric_readwrite_t *v); +#endif + + /** + * Read in Angles, Dihedrals, Impropers, and Cross Terms and optionally types. + * (Cross terms pertain to the CHARMM/NAMD CMAP feature) + */ + int (* read_angles)(void *handle, int *numangles, int **angles, int **angletypes, + int *numangletypes, char ***angletypenames, int *numdihedrals, + int **dihedrals, int **dihedraltypes, int *numdihedraltypes, + char ***dihedraltypenames, int *numimpropers, int **impropers, + int **impropertypes, int *numimpropertypes, char ***impropertypenames, + int *numcterms, int **cterms, int *ctermcols, int *ctermrows); + + /** + * Write out Angles, Dihedrals, Impropers, and Cross Terms + * (Cross terms pertain to the CHARMM/NAMD CMAP feature) + */ + int (* write_angles)(void *handle, int numangles, const int *angles, const int *angletypes, + int numangletypes, const char **angletypenames, int numdihedrals, + const int *dihedrals, const int *dihedraltypes, int numdihedraltypes, + const char **dihedraltypenames, int numimpropers, + const int *impropers, const int *impropertypes, int numimpropertypes, + const char **impropertypenames, int numcterms, const int *cterms, + int ctermcols, int ctermrows); + + + /** + * Retrieve metadata pertaining to timestep independent + * QM datasets in this file. + * + * The metadata are the sizes of the QM related data structure + * arrays that will be populated by the plugin when + * read_qm_rundata() is called. Since the allocation of these + * arrays is done by VMD rather than the plugin, VMD needs to + * know the sizes beforehand. Consequently read_qm_metadata() + * has to be called before read_qm_rundata(). + */ + int (* read_qm_metadata)(void *, molfile_qm_metadata_t *metadata); + + + /** + * Read timestep independent QM data. + * + * Typical data that are defined only once per trajectory are + * general info about the calculation (such as the used method), + * the basis set and normal modes. + * The data structures to be populated must have been allocated + * before by VMD according to sizes obtained through + * read_qm_metadata(). + */ + int (* read_qm_rundata)(void *, molfile_qm_t *qmdata); + + + /** + * Query the molfile plugin to determine whether or not memory + * allocations used for atomic coordinates and PBC unit cell information + * need to be aligned to a particular virtual memory or filesystem + * page size boundary to facilitate kernel-bypass unbuffered I/O, + * e.g., as used by jsplugin. This API should be called prior to the + * first call to read a timestep. The required page alignment size + * (in bytes) is returned to the caller. If this API has not been + * called, then the molfile plugin should revert to standard + * kernel-buffered I/O and suffer the associated performance loss. + * The caller can be assured that the plugin will not request any + * page alignment size that is greater than the value of + * MOLFILE_DIRECTIO_MAX_BLOCK_SIZE, both as a runtime sanity check, + * and to ensure that a caller that is unable to perform the max + * aligned allocation doesn't call the API in the first place. + * If a page-aligned allocation is not required for the file being read, + * the plugin will return an alignment size of 1. + */ +#if vmdplugin_ABIVERSION > 17 + int (* read_timestep_pagealign_size)(void *, int *pagealignsize); +#endif + + + /** + * Read the next timestep from the file. Return MOLFILE_SUCCESS, or + * MOLFILE_EOF on EOF. If the molfile_timestep_t or molfile_qm_metadata_t + * arguments are NULL, then the coordinate or qm data should be skipped. + * Otherwise, the application must prepare molfile_timestep_t and + * molfile_qm_timestep_t by allocating space for the corresponding + * number of coordinates, orbital wavefunction coefficients, etc. + * Since it is common for users to want to load only the final timestep + * data from a QM run, the application may provide any combination of + * valid, or NULL pointers for the molfile_timestep_t and + * molfile_qm_timestep_t parameters, depending on what information the + * user is interested in. + * The natoms and qm metadata parameters exist because some file formats + * cannot determine for themselves how many atoms etc are in a + * timestep; the app must therefore obtain this information elsewhere + * and provide it to the plugin. + */ + int (* read_timestep)(void *, int natoms, molfile_timestep_t *, + molfile_qm_metadata_t *, molfile_qm_timestep_t *); + + int (* read_timestep_metadata)(void *, molfile_timestep_metadata_t *); + int (* read_qm_timestep_metadata)(void *, molfile_qm_timestep_metadata_t *); + + +#if defined(EXPERIMENTAL_DIRECTIO_APIS) + /** + * Calculate file offsets and I/O lengths for performing + * kernel-bypass direct I/O or using GPU-Direct Storage APIs, + * thereby enabling peak I/O rates to be achieved for analysis + * worksloads like clustering of trajectories. + */ + int (* calc_fileoffsets_timestep)(void *, + molfile_ssize_t frameindex, + int firstatom, + int lastatom, + int *firstatom, + int *pageoffset, + molfile_ssize_t *startoffset, + molfile_ssize_t *readlen); +#endif + + +#if defined(DESRES_READ_TIMESTEP2) + /** + * Read a specified timestep! + */ + int (* read_timestep2)(void *, molfile_ssize_t index, molfile_timestep_t *); + + /** + * write up to count times beginning at index start into the given + * space. Return the number read, or -1 on error. + */ + molfile_ssize_t (* read_times)( void *, + molfile_ssize_t start, + molfile_ssize_t count, + double * times ); +#endif + + /** + * Console output, READ-ONLY function pointer. + * Function pointer that plugins can use for printing to the host + * application's text console. This provides a clean way for plugins + * to send message strings back to the calling application, giving the + * caller the ability to prioritize, buffer, and redirect console messages + * to an appropriate output channel, window, etc. This enables the use of + * graphical consoles like TkCon without losing console output from plugins. + * If the function pointer is NULL, no console output service is provided + * by the calling application, and the output should default to stdout + * stream. If the function pointer is non-NULL, all output will be + * subsequently dealt with by the calling application. + * + * XXX this should really be put into a separate block of + * application-provided read-only function pointers for any + * application-provided services + */ + int (* cons_fputs)(const int, const char*); + +} molfile_plugin_t; + +#endif + diff --git a/src/vmdplugin/vmddir.h b/src/vmdplugin/vmddir.h new file mode 100644 index 0000000000..f6b10d2781 --- /dev/null +++ b/src/vmdplugin/vmddir.h @@ -0,0 +1,161 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2016 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ + +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: vmddir.h,v $ + * $Author: johns $ $Locker: $ $State: Exp $ + * $Revision: 1.12 $ $Date: 2016/11/28 05:01:54 $ + * + ***************************************************************************/ + +#include +#include + +#if defined(_MSC_VER) || defined(__MINGW32__) +#include + +typedef struct { + HANDLE h; + WIN32_FIND_DATA fd; +} VMDDIR; + +#else +#include + +typedef struct { + DIR * d; +} VMDDIR; +#endif + + + +static VMDDIR * vmd_opendir(const char *); +static char * vmd_readdir(VMDDIR *); +static void vmd_closedir(VMDDIR *); +static int vmd_file_is_executable(const char * filename); + + +#define VMD_FILENAME_MAX 1024 + +#if defined(_MSC_VER) || defined(__MINGW32__) + +/* Windows version */ + +static VMDDIR * vmd_opendir(const char * filename) { + VMDDIR * d; + char dirname[VMD_FILENAME_MAX]; + + strcpy(dirname, filename); + strcat(dirname, "\\*"); + d = (VMDDIR *) malloc(sizeof(VMDDIR)); + if (d != NULL) { +#if _MSC_VER < 1400 || !(defined(UNICODE) || defined(_UNICODE)) + d->h = FindFirstFile(dirname, &(d->fd)); +#else + wchar_t szBuff[VMD_FILENAME_MAX]; + swprintf(szBuff, L"%p", dirname); + d->h = FindFirstFile(szBuff, &(d->fd)); +#endif + if (d->h == ((HANDLE)(-1))) { + free(d); + return NULL; + } + } + return d; +} + +static char * vmd_readdir(VMDDIR * d) { + if (FindNextFile(d->h, &(d->fd))) { +#if _MSC_VER < 1400 || !(defined(UNICODE) || defined(_UNICODE)) + return d->fd.cFileName; +#else + int len = wcslen(d->fd.cFileName); + char* ascii = new char[len + 1]; + wcstombs(ascii, d->fd.cFileName, len); + return ascii; +#endif + } + return NULL; +} + +static void vmd_closedir(VMDDIR * d) { + if (d->h != NULL) { + FindClose(d->h); + } + free(d); +} + + +static int vmd_file_is_executable(const char * filename) { + FILE * fp; + if ((fp=fopen(filename, "rb")) != NULL) { + fclose(fp); + return 1; + } + + return 0; +} + +#else + +/* Unix version */ + +#include +#include + +static VMDDIR * vmd_opendir(const char * filename) { + VMDDIR * d; + + d = (VMDDIR *) malloc(sizeof(VMDDIR)); + if (d != NULL) { + d->d = opendir(filename); + if (d->d == NULL) { + free(d); + return NULL; + } + } + + return d; +} + +static char * vmd_readdir(VMDDIR * d) { + struct dirent * p; + if ((p = readdir(d->d)) != NULL) { + return p->d_name; + } + + return NULL; +} + +static void vmd_closedir(VMDDIR * d) { + if (d->d != NULL) { + closedir(d->d); + } + free(d); +} + + +static int vmd_file_is_executable(const char * filename) { + struct stat buf; + if (!stat(filename, &buf)) { + if (buf.st_mode & S_IXUSR || + buf.st_mode & S_IXGRP || + buf.st_mode & S_IXOTH) { + return 1; + } + } + return 0; +} + +#endif + + + + diff --git a/src/vmdplugin/vmdplugin.h b/src/vmdplugin/vmdplugin.h new file mode 100644 index 0000000000..1bd32cf4d2 --- /dev/null +++ b/src/vmdplugin/vmdplugin.h @@ -0,0 +1,191 @@ +/*************************************************************************** + *cr + *cr (C) Copyright 1995-2006 The Board of Trustees of the + *cr University of Illinois + *cr All Rights Reserved + *cr + ***************************************************************************/ + +/*************************************************************************** + * RCS INFORMATION: + * + * $RCSfile: vmdplugin.h,v $ + * $Author: johns $ $Locker: $ $State: Exp $ + * $Revision: 1.34 $ $Date: 2018/05/02 03:12:56 $ + * + ***************************************************************************/ + +/** @file + * This header must be included by every VMD plugin library. It defines the + * API for every plugin so that VMD can organize the plugins it finds. + */ + +#ifndef VMD_PLUGIN_H +#define VMD_PLUGIN_H + + +/* + * Preprocessor tricks to make it easier for us to redefine the names of + * functions when building static plugins. + */ +#if !defined(VMDPLUGIN) +/** + * macro defining VMDPLUGIN if it hasn't already been set to the name of + * a static plugin that is being compiled. This is the catch-all case. + */ +#define VMDPLUGIN vmdplugin +#endif +/** concatenation macro, joins args x and y together as a single string */ +#define xcat(x, y) cat(x, y) +/** concatenation macro, joins args x and y together as a single string */ +#define cat(x, y) x ## y + +/* + * macros to correctly define plugin function names depending on whether + * the plugin is being compiled for static linkage or dynamic loading. + * When compiled for static linkage, each plugin needs to have unique + * function names for all of its entry points. When compiled for dynamic + * loading, the plugins must name their entry points consistently so that + * the plugin loading mechanism can find the register, register_tcl, init, + * and fini routines via dlopen() or similar operating system interfaces. + */ +/*@{*/ +/** Macro names entry points correctly for static linkage or dynamic loading */ +#define VMDPLUGIN_register xcat(VMDPLUGIN, _register) +#define VMDPLUGIN_register_tcl xcat(VMDPLUGIN, _register_tcl) +#define VMDPLUGIN_init xcat(VMDPLUGIN, _init) +#define VMDPLUGIN_fini xcat(VMDPLUGIN, _fini) +/*@}*/ + + +/** "WIN32" is defined on both WIN32 and WIN64 platforms... */ +#if (defined(WIN32)) +#define WIN32_LEAN_AND_MEAN +#include + +#if !defined(STATIC_PLUGIN) +#if defined(VMDPLUGIN_EXPORTS) +/** + * Only define DllMain for plugins, not in VMD or in statically linked plugins + * VMDPLUGIN_EXPORTS is only defined when compiling dynamically loaded plugins + */ +BOOL APIENTRY DllMain( HANDLE hModule, + DWORD ul_reason_for_call, + LPVOID lpReserved + ) +{ + return TRUE; +} + +#define VMDPLUGIN_API __declspec(dllexport) +#else +#define VMDPLUGIN_API __declspec(dllimport) +#endif /* VMDPLUGIN_EXPORTS */ +#else /* ! STATIC_PLUGIN */ +#define VMDPLUGIN_API +#endif /* ! STATIC_PLUGIN */ +#else +/** If we're not compiling on Windows, then this macro is defined empty */ +#define VMDPLUGIN_API +#endif + +/** define plugin linkage correctly for both C and C++ based plugins */ +#ifdef __cplusplus +#define VMDPLUGIN_EXTERN extern "C" VMDPLUGIN_API +#else +#define VMDPLUGIN_EXTERN extern VMDPLUGIN_API +#endif /* __cplusplus */ + +/* + * Plugin API functions start here + */ + + +/** + * Init routine: called the first time the library is loaded by the + * application and before any other API functions are referenced. + * Return 0 on success. + */ +VMDPLUGIN_EXTERN int VMDPLUGIN_init(void); + +/** + * Macro for creating a struct header used in all plugin structures. + * + * This header should be placed at the top of every plugin API definition + * so that it can be treated as a subtype of the base plugin type. + * + * abiversion: Defines the ABI for the base plugin type (not for other plugins) + * type: A string descriptor of the plugin type. + * name: A name for the plugin. + * author: A string identifier, possibly including newlines. + * Major and minor version. + * is_reentrant: Whether this library can be run concurrently with itself. + */ +#define vmdplugin_HEAD \ + int abiversion; \ + const char *type; \ + const char *name; \ + const char *prettyname; \ + const char *author; \ + int majorv; \ + int minorv; \ + int is_reentrant; + +/** + * Typedef for generic plugin header, individual plugins can + * make their own structures as long as the header info remains + * the same as the generic plugin header, most easily done by + * using the vmdplugin_HEAD macro. + */ +typedef struct { + vmdplugin_HEAD +} vmdplugin_t; + +/** + * Use this macro to initialize the abiversion member of each plugin + */ +#define vmdplugin_ABIVERSION 18 + +/*@{*/ +/** Use this macro to indicate a plugin's thread-safety at registration time */ +#define VMDPLUGIN_THREADUNSAFE 0 +#define VMDPLUGIN_THREADSAFE 1 +/*@}*/ + +/*@{*/ +/** Error return code for use in the plugin registration and init functions */ +#define VMDPLUGIN_SUCCESS 0 +#define VMDPLUGIN_ERROR -1 +/*@}*/ + +/** + * Function pointer typedef for register callback functions + */ +typedef int (*vmdplugin_register_cb)(void *, vmdplugin_t *); + +/** + * Allow the library to register plugins with the application. + * The callback should be called using the passed-in void pointer, which + * should not be interpreted in any way by the library. Each vmdplugin_t + * pointer passed to the application should point to statically-allocated + * or heap-allocated memory and should never be later modified by the plugin. + * Applications must be permitted to retain only a copy of the the plugin + * pointer, without making any deep copy of the items in the struct. + */ +VMDPLUGIN_EXTERN int VMDPLUGIN_register(void *, vmdplugin_register_cb); + +/** + * Allow the library to register Tcl extensions. + * This API is optional; if found by dlopen, it will be called after first + * calling init and register. + */ +VMDPLUGIN_EXTERN int VMDPLUGIN_register_tcl(void *, void *tcl_interp, + vmdplugin_register_cb); + +/** + * The Fini method is called when the application will no longer use + * any plugins in the library. + */ +VMDPLUGIN_EXTERN int VMDPLUGIN_fini(void); + +#endif /* VMD_PLUGIN_H */