Add ability to extract box information to a data set

doc/cpptraj.lyx

@@ -24119,14 +24119,27 @@ box
 
 \begin_layout LyX-Code
 box {[x <xval>] [y <yval>] [z <zval>] {[alpha <a>] [beta <b>] [gamma <g>]
+ |
 \end_layout
 
 \begin_layout LyX-Code
-	 [truncoct]} | nobox | auto [offset <offset>] [radii {vdw|gb|parse|none}]}
-\begin_inset Separator latexpar
-\end_inset
+	 [truncoct] [x <length>] |
+\end_layout
 
+\begin_layout LyX-Code
+     nobox |
+\end_layout
 
+\begin_layout LyX-Code
+     auto [offset <offset>] [radii {vdw|gb|parse|none}] |
+\end_layout
+
+\begin_layout LyX-Code
+     getbox {ucell|frac|shape} [name <setname>] [out <file>]
+\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
+
+<length>] 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
+
+<setname> The name of the 3x3 matrix data set.
+\end_layout
+
+\begin_layout Description
+out
+\begin_inset space ~
+\end_inset
+
+<file> 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

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 <offset>] [radii {vdw|gb|parse|none}]}\n",
+          "\t auto [offset <offset>] [radii {vdw|gb|parse|none}] |\n"
+          "\t getbox {ucell|frac|shape} [name <setname>] [out <file>]\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() );

src/Action_Box.h

@@ -14,14 +14,17 @@ 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).
     ModeType mode_;        ///< How box info will be assigned.
     double offset_;        ///< Offset for AUTO
     RadiiType radiiMode_;  ///< Radii type to use for AUTO
     std::vector<double> Radii_; ///< Hold radius for each atom for AUTO
+    DataSet* set_;              ///< Hold output data set for GETBOX
+    GetType getmode_;           ///> Mode for GETBOX
 };
 #endif

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
 }
 
 // -----------------------------------------------------------------------------

src/Version.h

@@ -12,7 +12,7 @@
  * Whenever a number that precedes <revision> 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

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 <<EOF
+parm ../tz2.ortho.parm7
+trajin ../tz2.ortho.nc
+box getbox ucell name UC out ortho.ucell.dat
+box getbox frac  name FC out ortho.frac.dat
+box getbox shape name SP out ortho.shape.dat
+run
+clear trajin
+clear parm
+parm ../tz2.truncoct.parm7
+trajin ../tz2.truncoct.nc
+box getbox ucell name UC1 out truncoct.ucell.dat
+box getbox frac  name FC1 out truncoct.frac.dat
+box getbox shape name SP1 out truncoct.shape.dat
+EOF
+RunCpptraj "$UNITNAME"
+DoTest ortho.ucell.dat.save ortho.ucell.dat
+DoTest ortho.frac.dat.save ortho.frac.dat
+DoTest ortho.shape.dat.save ortho.shape.dat 
+DoTest truncoct.ucell.dat.save truncoct.ucell.dat
+DoTest truncoct.frac.dat.save truncoct.frac.dat
+DoTest truncoct.shape.dat.save truncoct.shape.dat
+
 EndTest
 
 exit 0

test/Test_Box/ortho.frac.dat.save

@@ -0,0 +1,11 @@
+#Frame                                                                                                                     FC
+       1  0.028358513  0.000000000  0.000000000  0.000000000  0.023897405  0.000000000  0.000000000  0.000000000  0.027648269
+       2  0.028363710  0.000000000  0.000000000  0.000000000  0.023901784  0.000000000  0.000000000  0.000000000  0.027653335
+       3  0.028365167  0.000000000  0.000000000  0.000000000  0.023903013  0.000000000  0.000000000  0.000000000  0.027654757
+       4  0.028360979  0.000000000  0.000000000  0.000000000  0.023899483  0.000000000  0.000000000  0.000000000  0.027650673
+       5  0.028375125  0.000000000  0.000000000  0.000000000  0.023911403  0.000000000  0.000000000  0.000000000  0.027664464
+       6  0.028375014  0.000000000  0.000000000  0.000000000  0.023911310  0.000000000  0.000000000  0.000000000  0.027664357
+       7  0.028364030  0.000000000  0.000000000  0.000000000  0.023902054  0.000000000  0.000000000  0.000000000  0.027653648
+       8  0.028369573  0.000000000  0.000000000  0.000000000  0.023906725  0.000000000  0.000000000  0.000000000  0.027659052
+       9  0.028359252  0.000000000  0.000000000  0.000000000  0.023898028  0.000000000  0.000000000  0.000000000  0.027648989
+      10  0.028349880  0.000000000  0.000000000  0.000000000  0.023890130  0.000000000  0.000000000  0.000000000  0.027639852

test/Test_Box/ortho.shape.dat.save

@@ -0,0 +1,11 @@
+#Frame                                                                                                                     SP
+       1 35.262779662  0.000000000  0.000000000  0.000000000 41.845547680  0.000000000  0.000000000  0.000000000 36.168629529
+       2 35.256319297  0.000000000  0.000000000  0.000000000 41.837881310  0.000000000  0.000000000  0.000000000 36.162003206
+       3 35.254507178  0.000000000  0.000000000  0.000000000 41.835730909  0.000000000  0.000000000  0.000000000 36.160144536
+       4 35.259713496  0.000000000  0.000000000  0.000000000 41.841909130  0.000000000  0.000000000  0.000000000 36.165484598
+       5 35.242136053  0.000000000  0.000000000  0.000000000 41.821050373  0.000000000  0.000000000  0.000000000 36.147455615
+       6 35.242273127  0.000000000  0.000000000  0.000000000 41.821213036  0.000000000  0.000000000  0.000000000 36.147596211
+       7 35.255920399  0.000000000  0.000000000  0.000000000 41.837407946  0.000000000  0.000000000  0.000000000 36.161594061
+       8 35.249032888  0.000000000  0.000000000  0.000000000 41.829234692  0.000000000  0.000000000  0.000000000 36.154529620
+       9 35.261861086  0.000000000  0.000000000  0.000000000 41.844457626  0.000000000  0.000000000  0.000000000 36.167687356
+      10 35.273518528  0.000000000  0.000000000  0.000000000 41.858291250  0.000000000  0.000000000  0.000000000 36.179644260

test/Test_Box/ortho.ucell.dat.save

@@ -0,0 +1,11 @@
+#Frame                                                                                                                     UC
+       1 35.262779662  0.000000000  0.000000000  0.000000000 41.845547680  0.000000000  0.000000000  0.000000000 36.168629529
+       2 35.256319297  0.000000000  0.000000000  0.000000000 41.837881310  0.000000000  0.000000000  0.000000000 36.162003206
+       3 35.254507178  0.000000000  0.000000000  0.000000000 41.835730909  0.000000000  0.000000000  0.000000000 36.160144536
+       4 35.259713496  0.000000000  0.000000000  0.000000000 41.841909130  0.000000000  0.000000000  0.000000000 36.165484598
+       5 35.242136053  0.000000000  0.000000000  0.000000000 41.821050373  0.000000000  0.000000000  0.000000000 36.147455615
+       6 35.242273127  0.000000000  0.000000000  0.000000000 41.821213036  0.000000000  0.000000000  0.000000000 36.147596211
+       7 35.255920399  0.000000000  0.000000000  0.000000000 41.837407946  0.000000000  0.000000000  0.000000000 36.161594061
+       8 35.249032888  0.000000000  0.000000000  0.000000000 41.829234692  0.000000000  0.000000000  0.000000000 36.154529620
+       9 35.261861086  0.000000000  0.000000000  0.000000000 41.844457626  0.000000000  0.000000000  0.000000000 36.167687356
+      10 35.273518528  0.000000000  0.000000000  0.000000000 41.858291250  0.000000000  0.000000000  0.000000000 36.179644260

test/Test_Box/truncoct.frac.dat.save

@@ -0,0 +1,11 @@
+#Frame                                                                                                                    FC1
+       1  0.023563316  0.008330890  0.014429523 -0.000000000  0.024992671  0.014429523  0.000000000 -0.000000000  0.028859049
+       2  0.023569774  0.008333173  0.014433478 -0.000000000  0.024999521  0.014433478  0.000000000 -0.000000000  0.028866959
+       3  0.023564506  0.008331310  0.014430252 -0.000000000  0.024993933  0.014430252  0.000000000 -0.000000000  0.028860507
+       4  0.023569960  0.008333239  0.014433592 -0.000000000  0.024999718  0.014433592  0.000000000 -0.000000000  0.028867187
+       5  0.023568827  0.008332838  0.014432898 -0.000000000  0.024998516  0.014432898  0.000000000 -0.000000000  0.028865799
+       6  0.023569149  0.008332952  0.014433095 -0.000000000  0.024998857  0.014433095  0.000000000 -0.000000000  0.028866193
+       7  0.023568591  0.008332755  0.014432753 -0.000000000  0.024998266  0.014432753  0.000000000 -0.000000000  0.028865510
+       8  0.023563974  0.008331122  0.014429926 -0.000000000  0.024993368  0.014429926  0.000000000 -0.000000000  0.028859855
+       9  0.023558188  0.008329077  0.014426383 -0.000000000  0.024987232  0.014426383  0.000000000 -0.000000000  0.028852769
+      10  0.023569101  0.008332935  0.014433065 -0.000000000  0.024998806  0.014433065  0.000000000 -0.000000000  0.028866134

test/Test_Box/truncoct.shape.dat.save

@@ -0,0 +1,11 @@
+#Frame                                                                                                                    SP1
+       1 40.836801372 -8.167359385 -8.167359385 -8.167359385 40.836801372 -8.167359385 -8.167359385 -8.167359385 40.836801372
+       2 40.825611756 -8.165121462 -8.165121462 -8.165121462 40.825611756 -8.165121462 -8.165121462 -8.165121462 40.825611756
+       3 40.834738327 -8.166946776 -8.166946776 -8.166946776 40.834738327 -8.166946776 -8.166946776 -8.166946776 40.834738327
+       4 40.825289589 -8.165057028 -8.165057028 -8.165057028 40.825289589 -8.165057028 -8.165057028 -8.165057028 40.825289589
+       5 40.827252219 -8.165449554 -8.165449554 -8.165449554 40.827252219 -8.165449554 -8.165449554 -8.165449554 40.827252219
+       6 40.826695354 -8.165338181 -8.165338181 -8.165338181 40.826695354 -8.165338181 -8.165338181 -8.165338181 40.826695354
+       7 40.827661528 -8.165531416 -8.165531416 -8.165531416 40.827661528 -8.165531416 -8.165531416 -8.165531416 40.827661528
+       8 40.835661196 -8.167131350 -8.167131350 -8.167131350 40.835661196 -8.167131350 -8.167131350 -8.167131350 40.835661196
+       9 40.845690102 -8.169137131 -8.169137131 -8.169137131 40.845690102 -8.169137131 -8.169137131 -8.169137131 40.845690102
+      10 40.826778085 -8.165354728 -8.165354728 -8.165354728 40.826778085 -8.165354728 -8.165354728 -8.165354728 40.826778085

test/Test_Box/truncoct.ucell.dat.save

@@ -0,0 +1,11 @@
+#Frame                                                                                                                    UC1
+       1 42.438848535  0.000000000  0.000000000 -14.146281704 40.011730517  0.000000000 -14.146281704 -20.005862837 34.651176475
+       2 42.427219944  0.000000000  0.000000000 -14.142405507 40.000766976  0.000000000 -14.142405507 -20.000381067 34.641681770
+       3 42.436704555  0.000000000  0.000000000 -14.145567044 40.009709153  0.000000000 -14.145567044 -20.004852155 34.649425923
+       4 42.426885139  0.000000000  0.000000000 -14.142293905 40.000451319  0.000000000 -14.142293905 -20.000223239 34.641408402
+       5 42.428924764  0.000000000  0.000000000 -14.142973780 40.002374295  0.000000000 -14.142973780 -20.001184727 34.643073749
+       6 42.428346053  0.000000000  0.000000000 -14.142780876 40.001828681  0.000000000 -14.142780876 -20.000911920 34.642601233
+       7 42.429350131  0.000000000  0.000000000 -14.143115569 40.002775335  0.000000000 -14.143115569 -20.001385247 34.643421060
+       8 42.437663629  0.000000000  0.000000000 -14.145886735 40.010613376  0.000000000 -14.145886735 -20.005304267 34.650209003
+       9 42.448085974  0.000000000  0.000000000 -14.149360850 40.020439657  0.000000000 -14.149360850 -20.010217407 34.658718812
+      10 42.428432030  0.000000000  0.000000000 -14.142809535 40.001909741  0.000000000 -14.142809535 -20.000952450 34.642671433