diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index fe0f909900..c3883e2633 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -17014,7 +17014,8 @@ Options for XYZ format: \end_layout \begin_layout LyX-Code -ftype {atomxyz|xyz} titletype {none|single|perframe} width <#> prec <#> +[ftype {namexyz|atomxyz|xyz}] [titletype {none|single|perframe}] [width + <#>] [prec <#>] \end_layout \begin_deeper @@ -17023,7 +17024,16 @@ ftype \begin_inset space ~ \end_inset -{atomxyz|xyz} Choose either 'ATOM X Y Z' (default) or 'X Y Z' output format. +{atomxyz|xyz} Choose either 'NAME X Y Z' (default), 'ATOM X Y Z', or 'X + Y Z' output format. + +\series bold +'namexyz' +\series default + format is the standard XYZ format, where each frame is preceded by the + number of atoms and a comment. + The comment written by CPPTRAJ will include the set number and box information + (if present). \end_layout \begin_layout Description @@ -17033,6 +17043,11 @@ titletype {none|single|perframe} No title, one title (default), or title before every frame. + Only applies if not +\series bold +'namexyz' +\series default +. \end_layout \begin_layout Description diff --git a/src/Traj_XYZ.cpp b/src/Traj_XYZ.cpp index d41fb387c2..6d0808274c 100644 --- a/src/Traj_XYZ.cpp +++ b/src/Traj_XYZ.cpp @@ -15,45 +15,67 @@ Traj_XYZ::Traj_XYZ() : set_(0), width_(0), prec_(-1), - fmt_(0) + fmt_(0), + hasBox_(false) {} -/** \param line1 First line. Will be set to title line if title is present, cleared otherwise. +/** Parse the given line, try to determine what format it corresponds to. */ +Traj_XYZ::LineFmtType Traj_XYZ::DetermineLineFormat(std::string const& lineIn) +{ + if (lineIn.empty()) return UNKNOWN_LINE_FORMAT; + ArgList line(lineIn, " \t\r\n"); + if (line.Nargs() == 1 && validInteger(line[0])) + return SINGLE_INTEGER; + else if (line.Nargs() == 3 && validDouble(line[0]) && + validDouble(line[1]) && + validDouble(line[2])) + return THREE_DOUBLES; + else if (line.Nargs() == 4 && validInteger(line[0]) && + validDouble(line[1]) && + validDouble(line[2]) && + validDouble(line[3])) + return INTEGER_AND_THREE_DOUBLES; + else if (line.Nargs() == 4 && !validDouble(line[0]) && // TODO use isalpha(line[0][0])? + validDouble(line[1]) && + validDouble(line[2]) && + validDouble(line[3])) + return STRING_AND_THREE_DOUBLES; + else + return UNKNOWN_LINE_FORMAT; +} + +/** \param title Will be set to title line if title is present, cleared otherwise. + * \param line1 First line. * \param line2 Second line. + * \param line3 Third line. */ -Traj_XYZ::Type Traj_XYZ::DetermineFormat(std::string& line1, - std::string const& line2) -const +Traj_XYZ::Type Traj_XYZ::DetermineFormat(std::string& title, + std::string const& line1, + std::string const& line2, + std::string const& line3) { - std::string line = line1; - // This line will either be a title line, atom XYZ, or XYZ - RemoveLeadingWhitespace( line ); - RemoveTrailingWhitespace( line ); - if (!line.empty() && line[0] == '#') { - line1 = line; - // Ignore the title. Get the next line, which must be atom XYZ or XYZ - line = line2; - } else - line1.clear(); - static const unsigned int tkSize = 64; - char tk0[tkSize]; - char tk1[tkSize]; - char tk2[tkSize]; - char tk3[tkSize]; - int nscan = sscanf(line.c_str(), "%s %s %s %s", tk0, tk1, tk2, tk3); - //mprintf("DEBUG: '%s' '%s' '%s' '%s'\n", tk0, tk1, tk2, tk3); - if (nscan == 4) { - if ( validInteger( std::string(tk0) ) && - validDouble( std::string(tk1) ) && - validDouble( std::string(tk2) ) && - validDouble( std::string(tk3) ) ) + title.clear(); + LineFmtType l1fmt = DetermineLineFormat(line1); + LineFmtType l2fmt = DetermineLineFormat(line2); + LineFmtType l3fmt = DetermineLineFormat(line3); + + if (l1fmt == SINGLE_INTEGER && l2fmt == UNKNOWN_LINE_FORMAT && l3fmt == STRING_AND_THREE_DOUBLES) + { + title = line2; + return NAME_XYZ; + } else if (l1fmt == UNKNOWN_LINE_FORMAT) { + if (l2fmt == INTEGER_AND_THREE_DOUBLES) { + title = line1; return ATOM_XYZ; - } else if (nscan == 3) { - if ( validDouble( std::string(tk0) ) && - validDouble( std::string(tk1) ) && - validDouble( std::string(tk2) ) ) + } else if (l2fmt == THREE_DOUBLES) { + title = line1; return XYZ; - } + } + } else if (l1fmt == INTEGER_AND_THREE_DOUBLES) { + return ATOM_XYZ; + } else if (l1fmt == THREE_DOUBLES) + return XYZ; + return UNKNOWN; } @@ -63,7 +85,9 @@ bool Traj_XYZ::ID_TrajFormat(CpptrajFile& fileIn) { if (fileIn.OpenFile()) return false; std::string line1 = fileIn.GetLine(); std::string line2 = fileIn.GetLine(); - if ( DetermineFormat(line1, line2) == UNKNOWN ) return false; + std::string line3 = fileIn.GetLine(); + std::string title; + if ( DetermineFormat(title, line1, line2, line3) == UNKNOWN ) return false; return true; } @@ -71,9 +95,10 @@ bool Traj_XYZ::ID_TrajFormat(CpptrajFile& fileIn) { /** Print trajectory info to stdout. */ void Traj_XYZ::Info() { switch (ftype_) { - case UNKNOWN : - case XYZ : mprintf("is an XYZ trajectory"); break; + case UNKNOWN : mprintf("is an XYZ trajectory"); break; + case XYZ : mprintf("is an XYZ-only trajectory"); break; case ATOM_XYZ : mprintf("is an Atom-XYZ trajectory"); break; + case NAME_XYZ : mprintf("is a standard XYZ trajectory"); break; } } @@ -105,6 +130,37 @@ const char* Traj_XYZ::FMT_XYZ_ = "%lf %lf %lf"; const char* Traj_XYZ::FMT_ATOM_XYZ_ = "%*i %lf %lf %lf"; +const char* Traj_XYZ::FMT_NAME_XYZ_ = "%*s %lf %lf %lf"; + +/** Read box information from given line. */ +int Traj_XYZ::parseBoxLine(Box& xyzbox, ArgList const& boxline) +{ + double boxcrd[9]; + int iarg = 0; + while (iarg < boxline.Nargs()) { + int incr = 1; + if (boxline[iarg] == "X:") { + boxcrd[0] = convertToDouble(boxline[iarg+1]); + boxcrd[1] = convertToDouble(boxline[iarg+2]); + boxcrd[2] = convertToDouble(boxline[iarg+3]); + incr = 4; + } else if (boxline[iarg] == "Y:") { + boxcrd[3] = convertToDouble(boxline[iarg+1]); + boxcrd[4] = convertToDouble(boxline[iarg+2]); + boxcrd[5] = convertToDouble(boxline[iarg+3]); + incr = 4; + } else if (boxline[iarg] == "Z:") { + boxcrd[6] = convertToDouble(boxline[iarg+1]); + boxcrd[7] = convertToDouble(boxline[iarg+2]); + boxcrd[8] = convertToDouble(boxline[iarg+3]); + incr = 4; + } + iarg += incr; + } + xyzbox.SetupFromUcell(boxcrd); + return 0; +} + /** Set up trajectory for reading. * \return Number of frames in trajectory. */ @@ -112,50 +168,87 @@ int Traj_XYZ::setupTrajin(FileName const& fname, Topology* trajParm) { if (file_.OpenFileRead( fname )) return TRAJIN_ERR; // Initial format determiniation. + std::string title; std::string line1 = file_.GetLine(); std::string line2 = file_.GetLine(); - ftype_ = DetermineFormat(line1, line2); + std::string line3 = file_.GetLine(); + ftype_ = DetermineFormat(title, line1, line2, line3); switch (ftype_) { - case XYZ : fmt_ = FMT_XYZ_; break; - case ATOM_XYZ : fmt_ = FMT_ATOM_XYZ_; break; + case XYZ : + fmt_ = FMT_XYZ_; + mprintf("\tCoordinate lines are \n"); + break; + case ATOM_XYZ : + fmt_ = FMT_ATOM_XYZ_; + mprintf("\tCoordinate lines are <#> \n"); + break; + case NAME_XYZ : + fmt_ = FMT_NAME_XYZ_; + mprintf("\tCoordinate lines are \n"); + break; case UNKNOWN : mprinterr("Internal Error: '%s' does not appear to be XYZ format anymore.\n", file_.Filename().full()); return TRAJIN_ERR; } - if (line1.empty()) + // Initial header type determination + Box xyzbox; + if (title.empty()) titleType_ = NO_TITLE; - else + else if (ftype_ == NAME_XYZ) { + titleType_ = NATOM_COMMENT; + int nat = convertToInteger(line1); + if (nat != trajParm->Natom()) { + mprinterr("Error: # atoms in XYZ (%i) does not match # atoms in topology '%s' (%i)\n", + nat, trajParm->c_str(), trajParm->Natom()); + return TRAJIN_ERR; + } + // Determine if the comment line contains box information + // 0 1 2 3 4 5 6 + // Conf 1. Box X: 19.660 0.000 0.000 Y: 0.000 19.660 0.000 Z: 0.000 0.000 19.660 + ArgList boxline(line2, " \r\n\t"); + if (boxline.Nargs() > 12) { + if (boxline.hasKey("Box")) { + mprintf("\tComment line appears to have box information.\n"); + hasBox_ = true; + parseBoxLine(xyzbox, boxline); + } + } + } else titleType_ = SINGLE; - // Read one frame, see if the header is repeated. closeTraj(); int nframes = TRAJIN_UNK; - if (openTrajin()) return TRAJIN_ERR; - if (titleType_ != NO_TITLE) file_.Line(); - for (int at = 0; at != trajParm->Natom(); at++) { - const char* atline = file_.Line(); - if (atline == 0) { - mprinterr("Error: Unexpected EOF when reading first frame of '%s'\n", atline); - return TRAJIN_ERR; + + // Read one frame, see if the header is repeated. + if (titleType_ == SINGLE) { + if (openTrajin()) return TRAJIN_ERR; + if (titleType_ == SINGLE) + file_.Line(); + for (int at = 0; at != trajParm->Natom(); at++) { + const char* atline = file_.Line(); + if (atline == 0) { + mprinterr("Error: Unexpected EOF when reading first frame of '%s'\n", atline); + return TRAJIN_ERR; + } } + line2 = file_.GetLine(); + if (!line2.empty()) { + if (DetermineLineFormat(line2) == UNKNOWN_LINE_FORMAT) + titleType_ = MULTIPLE; + } else + nframes = 1; } - line2 = file_.GetLine(); - if (!line2.empty()) { - RemoveLeadingWhitespace( line2 ); - RemoveTrailingWhitespace( line2 ); - if (line2[0] == '#') - titleType_ = MULTIPLE; - } else - nframes = 1; + switch (titleType_) { case UNKNOWN_TITLE : - case NO_TITLE : mprintf("\tNo title detected.\n"); break; - case SINGLE : mprintf("\tSingle title detected.\n"); break; - case MULTIPLE : mprintf("\tTitle before each frame detected.\n"); break; + case NO_TITLE : mprintf("\tNo title detected.\n"); break; + case SINGLE : mprintf("\tSingle title detected.\n"); break; + case MULTIPLE : mprintf("\tTitle before each frame detected.\n"); break; + case NATOM_COMMENT : mprintf("\tStandard # atoms followed by comment detected.\n"); break; } - if (!line1.empty()) SetTitle( line1 ); - SetCoordInfo( CoordinateInfo(Box(), false, false, false) ); + if (!title.empty()) SetTitle( title ); + SetCoordInfo( CoordinateInfo(xyzbox, false, false, false) ); set_ = 0; @@ -163,16 +256,38 @@ int Traj_XYZ::setupTrajin(FileName const& fname, Topology* trajParm) } /** Read title. */ +void Traj_XYZ::ReadTitle(Box& xyzbox) { + if (titleType_ == SINGLE) { + file_.Line(); + titleType_ = NO_TITLE; + } else if (titleType_ == MULTIPLE) { + file_.Line(); + } else if (titleType_ == NATOM_COMMENT) { + const char* ptr = file_.Line(); // #atoms + if (ptr && hasBox_) { + std::string line2 = file_.GetLine(); // comment + ArgList boxline(line2, " \r\n\t"); + parseBoxLine(xyzbox, boxline); + } else + file_.Line(); // comment + } +} + +/** Read title. Skip box info */ void Traj_XYZ::ReadTitle() { if (titleType_ == SINGLE) { file_.Line(); titleType_ = NO_TITLE; - } else if (titleType_ == MULTIPLE) + } else if (titleType_ == MULTIPLE) { file_.Line(); + } else if (titleType_ == NATOM_COMMENT) { + file_.Line(); // #atoms + file_.Line(); // comment + } } /** Read specified frame into given buffer. */ -int Traj_XYZ::readXYZ(int set, int natom, double* xAddress) { +int Traj_XYZ::readXYZ(int set, int natom, double* xAddress, Box& xyzbox) { // If an earlier set is being requested, reopen the file. if (set < set_) { closeTraj(); @@ -185,7 +300,7 @@ int Traj_XYZ::readXYZ(int set, int natom, double* xAddress) { file_.Line(); set_++; } - ReadTitle(); + ReadTitle(xyzbox); // Read coordinates into frame double* xyz = xAddress; for (int at = 0; at != natom; at++) { @@ -200,24 +315,24 @@ int Traj_XYZ::readXYZ(int set, int natom, double* xAddress) { /** Read specified trajectory frame. */ int Traj_XYZ::readFrame(int set, Frame& frameIn) { - return readXYZ(set, frameIn.Natom(), frameIn.xAddress()); + return readXYZ(set, frameIn.Natom(), frameIn.xAddress(), frameIn.ModifyBox()); } /** Read velocities from specified frame. */ int Traj_XYZ::readVelocity(int set, Frame& frameIn) { - return readXYZ(set, frameIn.Natom(), frameIn.vAddress()); + return readXYZ(set, frameIn.Natom(), frameIn.vAddress(), frameIn.ModifyBox()); } /** Read forces from specified frame. */ int Traj_XYZ::readForce(int set, Frame& frameIn) { - return readXYZ(set, frameIn.Natom(), frameIn.fAddress()); + return readXYZ(set, frameIn.Natom(), frameIn.fAddress(), frameIn.ModifyBox()); } // ----------------------------------------------------------------------------- /** Write help. */ void Traj_XYZ::WriteHelp() { - mprintf("\tftype {atomxyz|xyz} : Choose either 'ATOM X Y Z' (default) or 'X Y Z' output format.\n" - "\ttitletype {none|single|perframe} : No title, one title (default), or title before every frame.\n" + mprintf("\tftype {namexyz|atomxyz|xyz} : Choose either 'NAME X Y Z' (default), 'ATOM X Y Z', or 'X Y Z' output format.\n" + "\ttitletype {none|single|perframe} : No title, one title (default), or title before every frame. Only applies if not 'namexyz'.\n" "\twidth <#> : Output format width.\n" "\tprec <#> : Output format precision.\n"); } @@ -230,6 +345,8 @@ int Traj_XYZ::processWriteArgs(ArgList& argIn, DataSetList const& DSLin) { ftype_ = XYZ; else if (arg == "atomxyz") ftype_ = ATOM_XYZ; + else if (arg == "namexyz") + ftype_ = NAME_XYZ; else { mprinterr("Error: Unrecognized key for 'ftype' = '%s'\n", arg.c_str()); return 1; @@ -260,13 +377,36 @@ int Traj_XYZ::setupTrajout(FileName const& fname, Topology* trajParm, { if (titleType_ == UNKNOWN_TITLE) titleType_ = SINGLE; + if (ftype_ == UNKNOWN) - ftype_ = ATOM_XYZ; + ftype_ = NAME_XYZ; + + if (ftype_ == NAME_XYZ) { + if (width_ == 0) + width_ = 14; + if (prec_ == -1) + prec_ = 8; + } TextFormat ffmt(TextFormat::DOUBLE, width_, prec_, 3); + + hasBox_ = cInfoIn.HasBox(); + if (hasBox_ && ftype_ != NAME_XYZ) { + mprintf("Warning: Box coordinates present but not using 'namexyz' format.\n" + "Warning: Box coordinates will not be written.\n"); + } + if (ftype_ == ATOM_XYZ) ofmt_ = "%i " + ffmt.Fmt(); else if (ftype_ == XYZ) ofmt_ = ffmt.Fmt(); + else if (ftype_ == NAME_XYZ) { + titleType_ = NATOM_COMMENT; + ofmt_ = "%-7s " + ffmt.Fmt(); + names_.clear(); + names_.reserve( trajParm->Natom() ); + for (Topology::atom_iterator at = trajParm->begin(); at != trajParm->end(); ++at) + names_.push_back( std::string(at->ElementName()) ); + } ofmt_.append("\n"); //mprintf("DEBUG: output format string: '%s'\n", ofmt_.c_str()); if (titleType_ != NO_TITLE && Title().empty()) @@ -280,8 +420,20 @@ int Traj_XYZ::writeFrame(int set, Frame const& frameOut) { if (titleType_ == SINGLE) { file_.Printf("#%s\n", Title().c_str()); titleType_ = NO_TITLE; - } else if (titleType_ == MULTIPLE) + } else if (titleType_ == MULTIPLE) { file_.Printf("#%s\n", Title().c_str()); + } else if (titleType_ == NATOM_COMMENT) { + file_.Printf("%i\n", frameOut.Natom()); + file_.Printf("Conf %i.", set+1); + if (hasBox_) { + Matrix_3x3 const& ucell = frameOut.BoxCrd().UnitCell(); + file_.Printf(" Box X: %.3f %.3f %.3f Y: %.3f %.3f %.3f Z: %.3f %.3f %.3f", + ucell[0], ucell[1], ucell[2], + ucell[3], ucell[4], ucell[5], + ucell[6], ucell[7], ucell[8]); + } + file_.Printf("\n"); + } const double* xyz = frameOut.xAddress(); if (ftype_ == ATOM_XYZ) { @@ -290,6 +442,9 @@ int Traj_XYZ::writeFrame(int set, Frame const& frameOut) { } else if (ftype_ == XYZ) { for (int at = 0; at != frameOut.Natom(); at++, xyz += 3) file_.Printf(ofmt_.c_str(), xyz[0], xyz[1], xyz[2]); + } else if (ftype_ == NAME_XYZ) { + for (int at = 0; at != frameOut.Natom(); at++, xyz += 3) + file_.Printf(ofmt_.c_str(), names_[at].c_str(), xyz[0], xyz[1], xyz[2]); } return 0; } diff --git a/src/Traj_XYZ.h b/src/Traj_XYZ.h index 5524e6bd93..b11335951d 100644 --- a/src/Traj_XYZ.h +++ b/src/Traj_XYZ.h @@ -35,24 +35,52 @@ class Traj_XYZ : public TrajectoryIO { void parallelCloseTraj(); // ------------------------------------------- # endif - - enum Type { UNKNOWN=0, XYZ, ATOM_XYZ }; - enum TitleType { NO_TITLE = 0, SINGLE, MULTIPLE, UNKNOWN_TITLE }; - - Type DetermineFormat(std::string&, std::string const&) const; + /// Data line types + enum Type { UNKNOWN=0, + XYZ, ///< + ATOM_XYZ, ///< <#> + NAME_XYZ, ///< + }; + /// Specify how/when headers should appear + enum TitleType { NO_TITLE = 0, ///< Never + SINGLE, ///< First frame only + MULTIPLE, ///< Every frame + NATOM_COMMENT, ///< Every frame has # atoms followed by a comment + UNKNOWN_TITLE + }; + /// Line format types + enum LineFmtType { SINGLE_INTEGER, + THREE_DOUBLES, + INTEGER_AND_THREE_DOUBLES, + STRING_AND_THREE_DOUBLES, + UNKNOWN_LINE_FORMAT }; + /// \return Format of given line + static LineFmtType DetermineLineFormat(std::string const&); + /// \return File format based on first three lines + static Type DetermineFormat(std::string&, std::string const&, std::string const&, + std::string const&); + /// Set given box from box line + static int parseBoxLine(Box&, ArgList const&); + /// Read past header, process box info + inline void ReadTitle(Box&); + /// Read past header inline void ReadTitle(); - int readXYZ(int, int, double*); + /// Read xyz coords + int readXYZ(int, int, double*, Box&); static const char* FMT_XYZ_; static const char* FMT_ATOM_XYZ_; + static const char* FMT_NAME_XYZ_; BufferedLine file_; - std::string ofmt_; - TitleType titleType_; - Type ftype_; - int set_; - int width_; - int prec_; - const char* fmt_; ///< Format for reading + std::string ofmt_; ///< Format for writing + TitleType titleType_; ///< Title type + Type ftype_; ///< Format type + int set_; ///< Current set (reading) + int width_; ///< Floating point width (writing) + int prec_; ///< Floating point precision (writing) + const char* fmt_; ///< Format for reading + bool hasBox_; ///< True if comment line contains box information + std::vector names_; ///< Names for writing NAME_XYZ format }; #endif diff --git a/src/Version.h b/src/Version.h index 3d3bbc5bea..91af6fcd4f 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 "V6.0.3" +#define CPPTRAJ_INTERNAL_VERSION "V6.1.0" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/test/Test_XYZformat/RunTest.sh b/test/Test_XYZformat/RunTest.sh index a3702c236e..b9ee5613a3 100755 --- a/test/Test_XYZformat/RunTest.sh +++ b/test/Test_XYZformat/RunTest.sh @@ -3,7 +3,7 @@ . ../MasterTest.sh CleanFiles xyz.in tz2.xyz test1.crd.save test?.crd tz2.st.xyz tz2.nt.at.xyz \ - tz2.nt.xyz tz2.mt.at.xyz tz2.mt.xyz + tz2.nt.xyz tz2.mt.at.xyz tz2.mt.xyz traj.crd cpptraj.traj.xyz TESTNAME='XYZ format tests' Requires notparallel @@ -16,7 +16,7 @@ parm ../tz2.parm7 trajin ../tz2.crd 1 10 trajout test1.crd.save crd # Single title, atom-xyz -trajout tz2.xyz +trajout tz2.xyz ftype atomxyz # Single title, xyz trajout tz2.st.xyz titletype single ftype xyz # No title, atom-xyz @@ -48,5 +48,24 @@ EOF ((N++)) done +# Read Standard format, With box coords in comments +cat > xyz.in < xyz.in <