Add ability to move and rotate (via RMS fitting on coordinates) grids

src/Action_Dipole.cpp

@@ -59,7 +59,7 @@ Action::RetType Action_Dipole::Init(ArgList& actionArgs, ActionInit& init, int d
   mprintf("\tGrid will be printed to file %s\n",outfile_->Filename().full());
   mprintf("\tMask expression: [%s]\n",mask_.MaskString());
   if (max_ > 0)
-    mprintf("\tOnly keeping density >= to %.0lf%% of the maximum density\n", max_);
+    mprintf("\tOnly keeping density >= to %.0f%% of the maximum density\n", max_);
 
   return Action::OK;
 }
@@ -104,12 +104,15 @@ Action::RetType Action_Dipole::Setup(ActionSetup& setup) {
 Action::RetType Action_Dipole::DoAction(int frameNum, ActionFrame& frm) {
   Vec3 cXYZ, dipolar_vector, COM;
 
+  // Move grid if necessary
+  MoveGrid( frm.Frm(), *grid_ );
+
   // Set up center to origin or box center
-  if (GridMode() == GridAction::BOX) 
+  if (GridOffsetType() == GridAction::BOX_CENTER) 
     cXYZ = frm.Frm().BoxCrd().Center();
-  else if (GridMode() == GridAction::MASKCENTER)
+  else if (GridOffsetType() == GridAction::MASK_CENTER)
     cXYZ = frm.Frm().VGeometricCenter( CenterMask() );
-  else // GridAction::ORIGIN/SPECIFIEDCENTER
+  else // GridAction::NO_OFFSET
     cXYZ.Zero();
 
   // Traverse over solvent molecules.
@@ -187,6 +190,8 @@ int Action_Dipole::SyncAction() {
   * comes with Midas/Plus.
   */
 void Action_Dipole::Print() {
+  if (grid_ == 0) return
+  FinishGrid( *grid_ );
   double max_density;
 
   // Write header

src/Action_Grid.cpp

@@ -128,11 +128,13 @@ Action::RetType Action_Grid::Setup(ActionSetup& setup) {
 
 // Action_Grid::DoAction()
 Action::RetType Action_Grid::DoAction(int frameNum, ActionFrame& frm) {
+  // Move grid if necessary
+  MoveGrid( frm.Frm(), *grid_ );
   if (useMaskArray_) {
     Vec3 offset(0.0);
-    if (GridMode() == BOX)
+    if (GridOffsetType() == BOX_CENTER)
       offset = frm.Frm().BoxCrd().Center();
-    else if (GridMode() == MASKCENTER)
+    else if (GridOffsetType() == MASK_CENTER)
       offset = frm.Frm().VGeometricCenter( CenterMask() );
     for (Cpptraj::MaskArray::const_iterator mask = mArray_.begin();
                                             mask != mArray_.end(); ++mask)
@@ -146,7 +148,8 @@ Action::RetType Action_Grid::DoAction(int frameNum, ActionFrame& frm) {
 
 // Action_Grid::print()
 void Action_Grid::Print() {
-  if (nframes_ < 1) return;
+  if (nframes_ < 1 || grid_ == 0) return;
+  FinishGrid( *grid_ );
   // Perform normalization and find max.
   double gridMax = 0.0;
   if (normalize_ == NONE) {

src/Action_GridFreeEnergy.cpp

@@ -82,12 +82,16 @@ Action::RetType Action_GridFreeEnergy::Setup(ActionSetup& setup) {
 
 // Action_GridFreeEnergy::action()
 Action::RetType Action_GridFreeEnergy::DoAction(int frameNum, ActionFrame& frm) {
+  // Move grid if necessary
+  MoveGrid( frm.Frm(), *grid_ );
   GridFrame( frm.Frm(), mask_, *grid_ );
   return Action::OK;
 }
 
 // Action_GridFreeEnergy::print()
 void Action_GridFreeEnergy::Print() {
+  if (grid_ == 0) return;
+  FinishGrid( *grid_ );
   /* How times does this occupancy count value arise?
    *    i.e. if  
    *                 voxelOccupancyCount[50] = 10 

src/Action_Vector.cpp

@@ -5,11 +5,13 @@
 #include "CpptrajStdio.h"
 #include "DistRoutines.h" // MinImagedVec, includes Matrix_3x3 for principal
 #include "DataSet_Vector.h"
+#include "DataSet_3D.h"
 
 // CONSTRUCTOR
 Action_Vector::Action_Vector() :
   Vec_(0),
   Magnitude_(0),
+  gridSet_(0),
   vcorr_(0),
   ptrajoutput_(false),
   needBoxInfo_(false),
@@ -20,12 +22,12 @@ Action_Vector::Action_Vector() :
 // Action_Vector::Help()
 void Action_Vector::Help() const {
   mprintf("\t[<name>] <Type> [out <filename> [ptrajoutput]] [<mask1>] [<mask2>]\n"
-          "\t[magnitude] [ired]\n"
+          "\t[magnitude] [ired] [gridset <grid>]\n"
           "\t<Type> = { mask     | minimage  | dipole | center   | corrplane | \n"
           "\t           box      | boxcenter | ucellx | ucelly   | ucellz    | \n"
           "\t           momentum | principal [x|y|z]  | velocity | force       }\n" 
-          "  Calculate the specified coordinate vector.\n"
-          "    mask: (Default) Vector from <mask1> to <mask2>.\n"
+          "  Calculate the vector of specified <Type>:\n"
+          "    mask             : (Default) Vector from <mask1> to <mask2>.\n"
           "    minimage         : Store the minimum image vector between atoms in <mask1> and <mask2>.\n"
           "    dipole           : Dipole and center of mass of the atoms specified in <mask1>\n"
           "    center           : Store the center of mass of atoms in <mask1>.\n"
@@ -137,6 +139,21 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d
   if (mode_ == BOX || mode_ == BOX_X || mode_ == BOX_Y || mode_ == BOX_Z ||
       mode_ == BOX_CTR || mode_ == MINIMAGE)
     needBoxInfo_ = true;
+  gridSet_ = 0;
+  if (needBoxInfo_) {
+    std::string gridSetArg = actionArgs.GetStringKey("gridset");
+    if (!gridSetArg.empty()) {
+      DataSetList gridSetList = init.DSL().SelectGroupSets( gridSetArg, DataSet::GRID_3D );
+      if (gridSetList.empty()) {
+        mprinterr("Error: %s does not select any grid data set.\n", gridSetArg.c_str());
+        return Action::ERR;
+      }
+      if (gridSetList.size() > 1) {
+        mprintf("Warning: %s selects more than 1 grid data set. Only using the first set.\n", gridSetArg.c_str());
+      }
+      gridSet_ = (DataSet_3D*)gridSetList[0];
+    }
+  }
   // Check if IRED vector
   bool isIred = actionArgs.hasKey("ired"); 
   // Vector Mask
@@ -190,6 +207,8 @@ Action::RetType Action_Vector::Init(ArgList& actionArgs, ActionInit& init, int d
     mprintf(" %s", filename.c_str());
   }
   mprintf("\n");
+  if (gridSet_ != 0)
+    mprintf("\tExtracting box vectors from grid set '%s'\n", gridSet_->legend());
 
   return Action::OK;
 }
@@ -418,16 +437,21 @@ void Action_Vector::CorrPlane(Frame const& currentFrame) {
 }
 
 //  Action_Vector::UnitCell()
-void Action_Vector::UnitCell(Box const& box) {
+void Action_Vector::UnitCell(Box const& box, Vec3 const& oxyz) {
   switch ( mode_ ) {
-    case BOX_X   : Vec_->AddVxyzo( box.UnitCell().Row1(), Vec3(0.0) ); break;
-    case BOX_Y   : Vec_->AddVxyzo( box.UnitCell().Row2(), Vec3(0.0) ); break;
-    case BOX_Z   : Vec_->AddVxyzo( box.UnitCell().Row3(), Vec3(0.0) ); break;
+    case BOX_X   : Vec_->AddVxyzo( box.UnitCell().Row1(), oxyz ); break;
+    case BOX_Y   : Vec_->AddVxyzo( box.UnitCell().Row2(), oxyz ); break;
+    case BOX_Z   : Vec_->AddVxyzo( box.UnitCell().Row3(), oxyz ); break;
     case BOX_CTR : Vec_->AddVxyz( box.UnitCell().TransposeMult(Vec3(0.5)) ); break;
     default: return;
   }
 }
 
+/** Store box vector (A, B, C) lengths as a single vector (|A|, |B|, |C|). */
+void Action_Vector::BoxLengths(Box const& box) {
+  Vec_->AddVxyz( box.Lengths() );
+}
+
 // Action_Vector::MinImage()
 void Action_Vector::MinImage(Frame const& frm) {
   Vec3 com1 = frm.VCenterOfMass(mask_);
@@ -461,11 +485,21 @@ Action::RetType Action_Vector::DoAction(int frameNum, ActionFrame& frm) {
     case PRINCIPAL_Y :
     case PRINCIPAL_Z : Principal(frm.Frm()); break;
     case CORRPLANE   : CorrPlane(frm.Frm()); break;
-    case BOX         : Vec_->AddVxyz( frm.Frm().BoxCrd().Lengths() ); break;
+    case BOX         :
+      if (gridSet_ != 0)
+        BoxLengths( gridSet_->Bin().GridBox() );
+      else
+        BoxLengths( frm.Frm().BoxCrd() );
+      break;
     case BOX_X       : 
     case BOX_Y       : 
     case BOX_Z       : 
-    case BOX_CTR     : UnitCell( frm.Frm().BoxCrd() ); break;
+    case BOX_CTR     :
+      if (gridSet_ != 0)
+        UnitCell( gridSet_->Bin().GridBox(), gridSet_->Bin().GridOrigin() );
+      else
+        UnitCell( frm.Frm().BoxCrd(), Vec3(0.0) );
+      break;
     case MINIMAGE    : MinImage( frm.Frm() ); break; 
     default          : return Action::ERR; // NO_OP
   } // END switch over vectorMode

src/Action_Vector.h

@@ -2,6 +2,7 @@
 #define INC_ACTION_VECTOR_H
 #include "Action.h"
 class DataSet_Vector;
+class DataSet_3D;
 class Action_Vector : public Action {
   public:
     Action_Vector();
@@ -29,11 +30,13 @@ class Action_Vector : public Action {
     void Dipole(Frame const&);
     void Principal(Frame const&);
     void CorrPlane(Frame const&);
-    void UnitCell(Box const&);
+    void UnitCell(Box const&, Vec3 const&);
+    void BoxLengths(Box const&);
     void MinImage(Frame const&);
 
     DataSet_Vector* Vec_;   ///< Hold vector values
     DataSet* Magnitude_;    ///< Hold vector magnitudes if requested
+    DataSet_3D* gridSet_;   ///< Hold grid set for getting box vectors from grid.
     double* vcorr_;         ///< Temp. space for calculating CorrPlane
     vectorMode mode_;       ///< Vector calculation mode
     bool ptrajoutput_;      ///< If true output in ptraj format

src/Action_Volmap.cpp

@@ -214,11 +214,10 @@ Action::RetType Action_Volmap::Init(ArgList& actionArgs, ActionInit& init, int d
       mprinterr("Error: Currently only works with orthogonal grids.\n");
       return Action::ERR;
     }
-    GridBin_Ortho const& gbo = static_cast<GridBin_Ortho const&>( grid_->Bin() );
-    dx_ = gbo.DX();
-    dy_ = gbo.DY();
-    dz_ = gbo.DZ();
-    Vec3 const& oxyz = gbo.GridOrigin();
+    dx_ = grid_->Bin().DX();
+    dy_ = grid_->Bin().DY();
+    dz_ = grid_->Bin().DZ();
+    Vec3 const& oxyz = grid_->Bin().GridOrigin();
     xmin_ = oxyz[0];
     ymin_ = oxyz[1];
     zmin_ = oxyz[2];

src/Box.cpp

@@ -68,11 +68,11 @@ void Box::swap(Box& rhs) {
 }
 
 #ifdef MPI
-/** Sync box info to Comm master. */
-int Box::SyncBox(Parallel::Comm const& commIn) {
+/** Broadcast box info from Comm master. */
+int Box::BroadcastBox(Parallel::Comm const& commIn) {
   commIn.MasterBcast( box_,    6, MPI_DOUBLE );
-  unitCell_.SyncMatrix( commIn );
-  fracCell_.SyncMatrix( commIn );
+  unitCell_.BroadcastMatrix( commIn );
+  fracCell_.BroadcastMatrix( commIn );
   commIn.MasterBcast( &cellVolume_, 1, MPI_DOUBLE );
   return 0;
 }

src/Box.h

@@ -24,7 +24,7 @@ class Box {
     /// SWAP
     void swap(Box&);
 #   ifdef MPI
-    int SyncBox(Parallel::Comm const&);
+    int BroadcastBox(Parallel::Comm const&);
     int SendBox(int, Parallel::Comm const&) const;
     int RecvBox(int, Parallel::Comm const&);
 #   endif

src/CoordinateInfo.cpp

@@ -136,7 +136,11 @@ std::string CoordinateInfo::InfoString() const {
 
 #ifdef MPI
 #define CINFOMPISIZE 12
-int CoordinateInfo::SyncCoordInfo(Parallel::Comm const& commIn) {
+/** Broadcast coordinate info from Master. The number of things to be 
+  * broadcast is controlled by CINFOMPISIZE; this must be increased
+  * if more items added to CoordinateInfo.
+  */
+int CoordinateInfo::BroadcastCoordInfo(Parallel::Comm const& commIn) {
   // ensSize, hasvel, hastemp, hastime, hasfrc, NrepDims, Dim1, ..., DimN, 
   int* iArray;
   int iSize;
@@ -181,7 +185,7 @@ int CoordinateInfo::SyncCoordInfo(Parallel::Comm const& commIn) {
       remdDim_.AddRemdDimension( iArray[ii] );
   }
   delete[] iArray;
-  box_.SyncBox( commIn );
+  box_.BroadcastBox( commIn );
   return 0;
 }
 #undef CINFOMPISIZE

src/CoordinateInfo.h

@@ -51,7 +51,7 @@ class CoordinateInfo {
     /// \return string containing info on present metadata
     std::string InfoString() const;
 #   ifdef MPI
-    int SyncCoordInfo(Parallel::Comm const&);
+    int BroadcastCoordInfo(Parallel::Comm const&);
 #   endif
     /// \return True if Frame would need to be re-setup based on CoordinateInfo
     bool operator !=(CoordinateInfo const& rhs) const {

src/DataIO_Std.cpp

@@ -795,13 +795,11 @@ int DataIO_Std::Read_3D(std::string const& fname,
     // Assume XYZ coords are of bin corners. Need to offset coords by half
     // the voxel size.
     if (!ds->Bin().IsOrthoGrid()) {
-      GridBin_Nonortho const& b = static_cast<GridBin_Nonortho const&>( ds->Bin() );
-      offset = b.Ucell().TransposeMult(Vec3( 1/(2*(double)ds->NX()),
-                                             1/(2*(double)ds->NY()),
-                                             1/(2*(double)ds->NZ()) ));
+      offset = ds->Bin().Ucell().TransposeMult(Vec3( 1/(2*(double)ds->NX()),
+                                                     1/(2*(double)ds->NY()),
+                                                     1/(2*(double)ds->NZ()) ));
     } else {
-      GridBin_Ortho const& b = static_cast<GridBin_Ortho const&>( ds->Bin() );
-      offset = Vec3(b.DX()/2, b.DY()/2, b.DZ()/2);
+      offset = Vec3(ds->Bin().DX()/2, ds->Bin().DY()/2, ds->Bin().DZ()/2);
     }
   }
   if (debug_ > 0)
@@ -1456,20 +1454,18 @@ int DataIO_Std::WriteSet3D( DataSet const& setIn, CpptrajFile& file ) {
                 set.Bin().GridOrigin()[1],
                 set.Bin().GridOrigin()[2]);
     if (set.Bin().IsOrthoGrid()) {
-      GridBin_Ortho const& b = static_cast<GridBin_Ortho const&>( set.Bin() );
-      file.Printf("#delta %12.7f %12.7f %12.7f\n", b.DX(), b.DY(), b.DZ());
+      file.Printf("#delta %12.7f %12.7f %12.7f\n", set.Bin().DX(), set.Bin().DY(), set.Bin().DZ());
     } else {
-      GridBin_Nonortho const& b = static_cast<GridBin_Nonortho const&>( set.Bin() );
       file.Printf("#delta %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f %12.7f\n",
-                  b.Ucell()[0]/set.NX(),
-                  b.Ucell()[1]/set.NX(),
-                  b.Ucell()[2]/set.NX(),
-                  b.Ucell()[3]/set.NY(),
-                  b.Ucell()[4]/set.NY(),
-                  b.Ucell()[5]/set.NY(),
-                  b.Ucell()[6]/set.NZ(),
-                  b.Ucell()[7]/set.NZ(),
-                  b.Ucell()[8]/set.NZ());
+                  set.Bin().Ucell()[0]/set.NX(),
+                  set.Bin().Ucell()[1]/set.NX(),
+                  set.Bin().Ucell()[2]/set.NX(),
+                  set.Bin().Ucell()[3]/set.NY(),
+                  set.Bin().Ucell()[4]/set.NY(),
+                  set.Bin().Ucell()[5]/set.NY(),
+                  set.Bin().Ucell()[6]/set.NZ(),
+                  set.Bin().Ucell()[7]/set.NZ(),
+                  set.Bin().Ucell()[8]/set.NZ());
     }
     file.Printf("#%s %s %s %s\n", Xdim.Label().c_str(), 
                 Ydim.Label().c_str(), Zdim.Label().c_str(), set.legend());