diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index b7f63c42df..fbcbeccff2 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -24119,14 +24119,27 @@ box \begin_layout LyX-Code box {[x ] [y ] [z ] {[alpha ] [beta ] [gamma ] + | \end_layout \begin_layout LyX-Code - [truncoct]} | nobox | auto [offset ] [radii {vdw|gb|parse|none}]} -\begin_inset Separator latexpar -\end_inset + [truncoct] [x ] | +\end_layout +\begin_layout LyX-Code + nobox | +\end_layout +\begin_layout LyX-Code + auto [offset ] [radii {vdw|gb|parse|none}] | +\end_layout + +\begin_layout LyX-Code + getbox {ucell|frac|shape} [name ] [out ] +\end_layout + +\begin_layout LyX-Code + } \end_layout \begin_deeper @@ -24183,7 +24196,15 @@ Change box length(s) to specified value(s). \end_layout \begin_layout Description -[truncoct] Set box angles to truncated octahedron. +[truncoct +\begin_inset space ~ +\end_inset + +[x +\begin_inset space ~ +\end_inset + +] Set box angles (and optionally a length) to truncated octahedron. \end_layout \begin_layout Description @@ -24213,6 +24234,44 @@ radii Born, PARSE, or no radii. \end_layout +\end_deeper +\begin_layout Description +getbox Save existing box information to a 3x3 matrix data set. +\end_layout + +\begin_deeper +\begin_layout Description +{ucell|frac|shape} Specify the kind of box information to save: +\series bold +ucell +\series default + saves the unit cell vectors, +\series bold +frac +\series default + saves the fractional unit cell vectors, and +\series bold +shape +\series default + saves the symmetric shape matrix unit cell vectors. +\end_layout + +\begin_layout Description +name +\begin_inset space ~ +\end_inset + + The name of the 3x3 matrix data set. +\end_layout + +\begin_layout Description +out +\begin_inset space ~ +\end_inset + + File to write the 3x3 matrix data to. +\end_layout + \end_deeper \end_deeper \begin_layout Standard @@ -24228,6 +24287,8 @@ Modify box information during trajectory processing. auto \series default keyword. + If 'getbox' is specified, the existing box information will be saved to + a data set. \end_layout \begin_layout Subsection diff --git a/src/Action_Box.cpp b/src/Action_Box.cpp index 50780b1483..cc2c03e5c3 100644 --- a/src/Action_Box.cpp +++ b/src/Action_Box.cpp @@ -2,27 +2,36 @@ #include "CpptrajStdio.h" Action_Box::Action_Box() : - mode_(SET) + mode_(SET), + set_(0), + getmode_(GET_UNITCELL) {} void Action_Box::Help() const { mprintf("\t{%s |\n" "\t %s |\n" "\t nobox |\n" - "\t auto [offset ] [radii {vdw|gb|parse|none}]}\n", + "\t auto [offset ] [radii {vdw|gb|parse|none}] |\n" + "\t getbox {ucell|frac|shape} [name ] [out ]\n" + "\t}\n", BoxArgs::Keywords_XyzAbg(), BoxArgs::Keywords_TruncOct()); mprintf(" For each input frame, replace any box information with the information given.\n" " If 'truncoct' is specified, alpha, beta, and gamma will be set to the\n" - " appropriate angle for a truncated octahedral box. If 'nobox' is specified,\n" - " all existing box information will be removed. If 'auto' is specified, an\n" - " orthogonal box will be set for existing atoms using the specified distance\n" - " offset value, ensuring specified radii (default vdw) are enclosed.\n"); + " appropriate angle for a truncated octahedral box.\n" + " If 'nobox' is specified, all existing box information will be removed.\n" + " If 'auto' is specified, an orthogonal box will be set for existing atoms\n" + " using the specified distance offset value, ensuring specified radii\n" + " (default vdw) are enclosed.\n" + " If 'getbox' is specified, the existing box information will be saved to\n" + " a data set.\n"); } // Action_Box::Init() Action::RetType Action_Box::Init(ArgList& actionArgs, ActionInit& init, int debugIn) { // Get keywords + std::string dsname; + DataFile* outfile = 0; if ( actionArgs.hasKey("nobox") ) mode_ = REMOVE; else if (actionArgs.hasKey("auto")) { @@ -52,12 +61,37 @@ Action::RetType Action_Box::Init(ArgList& actionArgs, ActionInit& init, int debu return Action::ERR; } } + } else if (actionArgs.Contains("getbox")) { + mode_ = GETBOX; + std::string getbox = actionArgs.GetStringKey("getbox"); + if (getbox == "ucell") + getmode_ = GET_UNITCELL; + else if (getbox == "frac") + getmode_ = GET_FRACCELL; + else if (getbox == "shape") + getmode_ = GET_SHAPE; + else { + mprinterr("Error: Expected getbox {ucell|frac|shape}, got %s\n", getbox.c_str()); + return Action::ERR; + } + dsname = actionArgs.GetStringKey("name"); + outfile = init.DFL().AddDataFile( actionArgs.GetStringKey("out"), actionArgs ); } else { mode_ = SET; // TODO check for bad args? if (boxArgs_.SetBoxArgs( actionArgs )) return Action::ERR; } + // Create set + if (mode_ == GETBOX) { + set_ = init.DSL().AddSet( DataSet::MAT3X3, dsname, "BOX" ); + if (set_ == 0) { + mprinterr("Error: Could not create box set.\n"); + return Action::ERR; + } + if (outfile != 0) outfile->AddDataSet( set_ ); + } + mprintf(" BOX:"); if (mode_ == REMOVE) mprintf(" Removing box information.\n"); @@ -72,6 +106,11 @@ Action::RetType Action_Box::Init(ArgList& actionArgs, ActionInit& init, int debu mprintf("\tWill use VDW, GB, or PARSE radii if available (with that priority).\n"); break; } + } else if (mode_ == GETBOX) { + static const char* getstr[] = { "unit cell", "fractional cell", "shape matrix" }; + mprintf(" Getting %s information from box.\n", getstr[getmode_]); + mprintf("\tOutput set: %s\n", set_->legend()); + if (outfile != 0) mprintf("\tOutput file: %s\n", outfile->DataFilename().full()); } else { boxArgs_.PrintXyzAbg(); } @@ -84,6 +123,13 @@ Action::RetType Action_Box::Setup(ActionSetup& setup) { if (mode_ == REMOVE) { mprintf("\tRemoving box info.\n"); cInfo_.SetBox( Box() ); + } else if (mode_ == GETBOX) { + // Check box type + if (!setup.CoordInfo().TrajBox().HasBox()) { + mprintf("Warning: Topology %s does not contain box information.\n", + setup.Top().c_str()); + return Action::SKIP; + } } else { // SET, AUTO // Fill in missing box information from current box @@ -143,6 +189,22 @@ Action::RetType Action_Box::DoAction(int frameNum, ActionFrame& frm) { frm.ModifyFrm().ModifyBox().SetNoBox(); } else if (mode_ == AUTO) { frm.ModifyFrm().SetOrthoBoundingBox(Radii_, offset_); + } else if (mode_ == GETBOX) { + if (getmode_ == GET_UNITCELL) { + Matrix_3x3 const& ucell = frm.Frm().BoxCrd().UnitCell(); + set_->Add( frameNum, ucell.Dptr() ); + } else if (getmode_ == GET_FRACCELL) { + Matrix_3x3 const& frac = frm.Frm().BoxCrd().FracCell(); + set_->Add( frameNum, frac.Dptr() ); + } else if (getmode_ == GET_SHAPE) { + double shape[6]; + frm.Frm().BoxCrd().GetSymmetricShapeMatrix(shape); + double ucell[9]; + ucell[0] = shape[0]; ucell[1] = shape[1]; ucell[2] = shape[3]; + ucell[3] = shape[1]; ucell[4] = shape[2]; ucell[5] = shape[4]; + ucell[6] = shape[3]; ucell[7] = shape[4]; ucell[8] = shape[5]; + set_->Add( frameNum, ucell ); + } } else { boxArgs_.SetMissingInfo( frm.Frm().BoxCrd() ); frm.ModifyFrm().ModifyBox().AssignFromXyzAbg( boxArgs_.XyzAbg() ); diff --git a/src/Action_Box.h b/src/Action_Box.h index 31ba075489..155d596c81 100644 --- a/src/Action_Box.h +++ b/src/Action_Box.h @@ -14,8 +14,9 @@ class Action_Box : public Action { Action::RetType DoAction(int, ActionFrame&); void Print() {} - enum ModeType { SET = 0, REMOVE, AUTO }; + enum ModeType { SET = 0, REMOVE, AUTO, GETBOX }; enum RadiiType { UNSPECIFIED = 0, GB, PARSE, VDW, NONE }; + enum GetType { GET_UNITCELL = 0, GET_FRACCELL, GET_SHAPE }; CoordinateInfo cInfo_; ///< For holding modified coordinate info. BoxArgs boxArgs_; ///< Hold arguments for setting box (SET). @@ -23,5 +24,7 @@ class Action_Box : public Action { double offset_; ///< Offset for AUTO RadiiType radiiMode_; ///< Radii type to use for AUTO std::vector Radii_; ///< Hold radius for each atom for AUTO + DataSet* set_; ///< Hold output data set for GETBOX + GetType getmode_; ///> Mode for GETBOX }; #endif diff --git a/src/Box.cpp b/src/Box.cpp index bad5d9c7b3..d7d67793eb 100644 --- a/src/Box.cpp +++ b/src/Box.cpp @@ -490,12 +490,12 @@ void Box::CalcShapeFromXyzAbg(double* shape, const double* box) double B = sqrt( Evals[1] ); double C = sqrt( Evals[2] ); - shape[0] = A*HtH[0]*HtH[0] + B*HtH[1]*HtH[1] + C*HtH[2]*HtH[2]; - shape[2] = A*HtH[3]*HtH[3] + B*HtH[4]*HtH[4] + C*HtH[5]*HtH[5]; - shape[5] = A*HtH[6]*HtH[6] + B*HtH[7]*HtH[7] + C*HtH[8]*HtH[8]; - shape[1] = A*HtH[0]*HtH[3] + B*HtH[1]*HtH[4] + C*HtH[2]*HtH[5]; - shape[3] = A*HtH[0]*HtH[6] + B*HtH[1]*HtH[7] + C*HtH[2]*HtH[8]; - shape[4] = A*HtH[3]*HtH[6] + B*HtH[4]*HtH[7] + C*HtH[5]*HtH[8]; + shape[0] = A*HtH[0]*HtH[0] + B*HtH[1]*HtH[1] + C*HtH[2]*HtH[2]; // XX + shape[2] = A*HtH[3]*HtH[3] + B*HtH[4]*HtH[4] + C*HtH[5]*HtH[5]; // YY + shape[5] = A*HtH[6]*HtH[6] + B*HtH[7]*HtH[7] + C*HtH[8]*HtH[8]; // ZZ + shape[1] = A*HtH[0]*HtH[3] + B*HtH[1]*HtH[4] + C*HtH[2]*HtH[5]; // XY + shape[3] = A*HtH[0]*HtH[6] + B*HtH[1]*HtH[7] + C*HtH[2]*HtH[8]; // XZ + shape[4] = A*HtH[3]*HtH[6] + B*HtH[4]*HtH[7] + C*HtH[5]*HtH[8]; // YZ } // ----------------------------------------------------------------------------- diff --git a/src/Version.h b/src/Version.h index 0e0d25aa2f..f4faf6227f 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.23.0" +#define CPPTRAJ_INTERNAL_VERSION "V6.23.1" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/test/Test_Box/RunTest.sh b/test/Test_Box/RunTest.sh index cc2524bf42..8b39745135 100755 --- a/test/Test_Box/RunTest.sh +++ b/test/Test_Box/RunTest.sh @@ -4,7 +4,8 @@ CleanFiles box.in addbox.rst7 addbox.rst7.? addbox.rst7.10 \ modX.rst7 modX.rst7.? modX.rst7.10 \ - frame1.rst7 tz2.box.rst7 tz2.vdw.rst7 + frame1.rst7 tz2.box.rst7 tz2.vdw.rst7 \ + *.ucell.dat *.frac.dat *.shape.dat TESTNAME='Box tests' Requires netcdf maxthreads 10 @@ -79,6 +80,30 @@ EOF DoTest tz2.vdw.rst7.save tz2.vdw.rst7 fi +UNITNAME='Box test (get box info)' +cat > box.in <