Add 'restrict' keyword to density action

doc/cpptraj.lyx

@@ -21500,29 +21500,33 @@ density
 \end_layout
 
 \begin_layout LyX-Code
-density [out <filename>] 
+density [out <filename>] [name <set name>]
 \end_layout
 
 \begin_layout LyX-Code
         [ <mask1> ...
- <maskN> [name <set name>] [delta <resolution>] [{x|y|z}]
+ <maskN> [delta <resolution>] [{x|y|z}]
 \end_layout
 
 \begin_layout LyX-Code
-          [{number|mass|charge|electron}] [{bincenter|binedge}] ]
+          [{number|mass|charge|electron}] [{bincenter|binedge}]
 \begin_inset Separator latexpar
 \end_inset
 
 
 \end_layout
 
+\begin_layout LyX-Code
+          [restrict {cylinder|square} cutoff <cut>] ]
+\end_layout
+
 \begin_deeper
 \begin_layout Description
 [out
 \begin_inset space ~
 \end_inset
 
-<filename>] Output file for histogram (relative distances vs.
+<filename>] Output file for histogram(s) (relative distances vs.
  densities for each mask) or total density.
 \end_layout
 
@@ -21564,7 +21568,7 @@ name>] Output data set name.
 \end_layout
 
 \begin_layout Description
-[{x|y|z}] Coordinate for density calculation.
+[{x|y|z}] Coordinate (axis) for density calculation.
  (default z)
 \end_layout
 
@@ -21581,22 +21585,70 @@ name>] Output data set name.
  based on bin center (default) or bin edges.
 \end_layout
 
+\begin_layout Description
+[restrict
+\begin_inset space ~
+\end_inset
+
+{cylinder|square}] If 
+\series bold
+'restrict'
+\series default
+ is specified, only calculate the density that is within a cylinder or square
+ shape from the specified axis as defined by a distance cutoff.
+\end_layout
+
+\begin_deeper
+\begin_layout Description
+cutoff
+\begin_inset space ~
+\end_inset
+
+<cut> The distance cutoff for 
+\series bold
+'restrict'
+\series default
+.
+\end_layout
+
+\end_deeper
 \end_deeper
 \begin_layout Standard
-DataSet Aspects:
+DataSets Created if masks specified:
 \end_layout
 
 \begin_layout Description
-[avg] Average density over coordinate.
+<set
+\begin_inset space ~
+\end_inset
+
+name>[avg]:<idx> Average density over coordinate for mask number <idx>.
+\end_layout
+
+\begin_layout Description
+<set
+\begin_inset space ~
+\end_inset
+
+name>[sd]:<idx> Standard deviation of density over coordinate for mask number
+ <idx>.
+\end_layout
+
+\begin_layout Standard
+DataSets Created if no masks specified:
 \end_layout
 
 \begin_layout Description
-[sd] Standard deviation of density over coordinate.
+<set
+\begin_inset space ~
+\end_inset
+
+name> Total system density each frame.
 \end_layout
 
 \end_deeper
 \begin_layout Standard
-If no arguments are specified, calculate the total system density.
+If no atom masks are specified, calculate the total system density.
  Otherwise, calculate specified density along the given axis for atoms in
  specified mask(s).
  Defaults are shown in parentheses above.
@@ -21662,6 +21714,15 @@ density out ion_density.dat number delta 0.25 ":SOD" ":CLA"
 See also $AMBERHOME/AmberTools/test/cpptraj/Test_Density.
 \end_layout
 
+\begin_layout Standard
+It can be useful to write out the average and standard deviation as an XYDY
+ set to a Grace data file, e.g.
+\end_layout
+
+\begin_layout LyX-Code
+density :WAT@O out wato.agr xydy
+\end_layout
+
 \begin_layout Subsection
 diffusion
 \end_layout

src/Action_Density.cpp

@@ -14,15 +14,18 @@ Action_Density::Action_Density() :
   property_(NUMBER),
   binType_(CENTER),
   delta_(0.0),
-  density_(0)
+  density_(0),
+  restrictType_(NONE),
+  cutVal_(-1)
 {
   area_coord_[0] = DX; area_coord_[1] = DY;
 }
 
 void Action_Density::Help() const {
   mprintf("\t[out <filename>] [name <set name>]\n"
           "\t[ <mask1> ... <maskN> [delta <resolution>] [{x|y|z}]\n"
-          "\t  [{number|mass|charge|electron}] [{bincenter|binedge}] ]\n"
+          "\t  [{number|mass|charge|electron}] [{bincenter|binedge}]\n"
+          "\t  [restrict {cylinder|square} cutoff <cut>]  ]\n"
           "  If one or more masks are specified, calculate specified density of selected\n"
           "  atoms along a coordinate. Otherwise calculate the total system density.\n");
 }
@@ -57,6 +60,26 @@ Action::RetType Action_Density::Init(ArgList& actionArgs, ActionInit& init, int
     area_coord_[1] = DY;
   }
 
+  restrictType_ = NONE;
+  std::string restrictStr = actionArgs.GetStringKey("restrict");
+  if (!restrictStr.empty()) {
+    if (restrictStr == "cylinder")
+      restrictType_ = CYLINDER;
+    else if (restrictStr == "square")
+      restrictType_ = SQUARE;
+    else {
+      mprinterr("Error: Unrecognized option for 'restrict': %s\n", restrictStr.c_str());
+      return Action::ERR;
+    }
+  }
+  if (restrictType_ != NONE) {
+    cutVal_ = actionArgs.getKeyDouble("cutoff", -1);
+    if (cutVal_ <= 0) {
+      mprinterr("Error: If 'restrict' is specified, 'cutoff' must be specified and > 0\n");
+      return Action::ERR;
+    }
+  }
+
   property_ = NUMBER;
   if      (actionArgs.hasKey("number") )   property_ = NUMBER;
   else if (actionArgs.hasKey("mass") )     property_ = MASS;
@@ -133,6 +156,18 @@ Action::RetType Action_Density::Init(ArgList& actionArgs, ActionInit& init, int
     mprintf("\tData set name is '%s'\n", dsname.c_str());
     mprintf("\tData set aspect [avg] holds the mean, aspect [sd] holds standard deviation.\n");
     mprintf("\tBin coordinates will be to bin %s.\n", binStr[binType_]);
+    if (restrictType_ != NONE) {
+      if (restrictType_ == CYLINDER) {
+        mprintf("\tOnly atoms within a radius of %f Ang. from the axis will be binned.\n", cutVal_);
+        // Square the cutoff
+        cutVal_ *= cutVal_;
+      } else if (restrictType_ == SQUARE) {
+        mprintf("\tOnly atoms within a square of %f Ang. from the axis will be binned.\n", cutVal_);
+      } else {
+        mprinterr("Internal Error: Restrict type not handled.\n");
+        return Action::ERR;
+      }
+    }
   } else {
     mprintf(" No masks specified, calculating total system density in g/cm^3.\n");
     mprintf("\tData set name is '%s'\n", density_->legend());
@@ -149,11 +184,16 @@ Action::RetType Action_Density::HistSetup(ActionSetup& setup) {
   properties_.clear();
   properties_.reserve( masks_.size() );
 
+  int total_nselected = 0;
   for (std::vector<AtomMask>::iterator mask = masks_.begin();
                                        mask != masks_.end();
                                        mask++)
   {
     if (setup.Top().SetupIntegerMask(*mask) ) return Action::ERR;
+    mprintf("\t");
+    mask->BriefMaskInfo();
+    mprintf("\n");
+    total_nselected += mask->Nselected();
 
     std::vector<double> property;
 
@@ -172,9 +212,10 @@ Action::RetType Action_Density::HistSetup(ActionSetup& setup) {
 
     properties_.push_back(property);
 
-    mprintf("\t");
-    mask->BriefMaskInfo();
-    mprintf("\n");
+  }
+  if (total_nselected < 1) {
+    mprintf("Warning: Nothing selected by masks, skipping.\n");
+    return Action::SKIP;
   }
 
   return Action::OK;  
@@ -220,19 +261,39 @@ Action::RetType Action_Density::HistAction(int frameNum, ActionFrame& frm) {
     for (AtomMask::const_iterator atm = mask.begin(); atm != mask.end(); ++atm, midx++)
     {
       const double* XYZ = frm.Frm().XYZ( *atm );
-      if (XYZ[axis_] < 0) {
-        // Coordinate is negative. Need to subtract off delta so that proper bin
-        // is populated (-delta to 0.0).
-        bin = (XYZ[axis_] - delta_) / delta_;
-      } else {
-        // Coordinate is positive.
-        bin = XYZ[axis_] / delta_;
+      // Check if we are doing a geometric restriction
+      bool binValues = false;
+      if (restrictType_ == CYLINDER) {
+        //mprintf("DEBUG: CYLINDER RESTRICT: atom %i ordinates %i %i\n", *atm+1, area_coord_[0], area_coord_[1]);
+        double x = XYZ[area_coord_[0]];
+        double y = XYZ[area_coord_[1]];
+        //mprintf("DEBUG:    coords: x= %f  y= %f\n", x, y);
+        double d2 = x*x + y*y;
+        if (d2 < cutVal_)
+          binValues = true;
+        //mprintf("DEBUG:    binValues= %i\n", (int)binValues);
+      } else if (restrictType_ == SQUARE) {
+        if (fabs(XYZ[area_coord_[0]]) < cutVal_ && 
+            fabs(XYZ[area_coord_[1]]) < cutVal_)
+          binValues = true;
+      } else
+        binValues = true;
+      if (binValues) {
+        if (XYZ[axis_] < 0) {
+          // Coordinate is negative. Need to subtract off delta so that proper bin
+          // is populated (-delta to 0.0).
+          bin = (XYZ[axis_] - delta_) / delta_;
+        } else {
+          // Coordinate is positive.
+          bin = XYZ[axis_] / delta_;
+        }
+        //mprintf("DEBUG: frm=%6i mask=%3u atm=%8i crd=%8.3f bin=%li\n", frameNum+1, idx, *atm+1, XYZ[axis_], bin);
+        Sum[bin] += properties_[idx][midx];
       }
-      //mprintf("DEBUG: frm=%6i mask=%3u atm=%8i crd=%8.3f bin=%li\n", frameNum+1, idx, *atm+1, XYZ[axis_], bin);
-      Sum[bin] += properties_[idx][midx];
-    }
+    } // END loop over atoms in mask
     // Accumulate sums
-    hist.accumulate( Sum );
+    if (!Sum.empty())
+      hist.accumulate( Sum );
   }
 
   // Accumulate area

src/Action_Density.h

@@ -10,7 +10,8 @@
 #include "OnlineVarT.h"
 #include "ImagedAction.h"
 
-#define ROUTINE_VERSION_STRING "1.0.2"
+#define ROUTINE_VERSION_STRING "1.0.3"
+// 1.0.3 - DRR Add 'restrict'
 
 /** Calculate density along a coordinate.
   * \author Hannes H. Loeffler.
@@ -58,6 +59,7 @@ class Action_Density : public Action {
   enum DirectionType {DX = 0, DY, DZ};
   enum PropertyType {NUMBER = 0, MASS, CHARGE, ELECTRON};
   enum BinCoordType {CENTER = 0, EDGE};
+  enum RestrictType { NONE=0, CYLINDER, SQUARE };
 
   DirectionType axis_;          ///< Which axis to bin along.
   DirectionType area_coord_[2]; ///< Hold which two axes used to calc. area
@@ -77,6 +79,8 @@ class Action_Density : public Action {
 
   DataSet* density_;            ///< Hold total system density (if not binning)
   ImagedAction image_;          ///< Used to calculate system volume for total density.
+  RestrictType restrictType_;   ///< Used to restrict calculation to a certain shape.
+  double cutVal_;               ///< Cutoff to use if shape restriction in use.
 # ifdef MPI
   Parallel::Comm trajComm_;
 # endif

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 "V4.32.1"
+#define CPPTRAJ_INTERNAL_VERSION "V4.32.2"
 /// PYTRAJ relies on this
 #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION
 #endif

test/Test_Density/RunTest.sh

@@ -8,7 +8,7 @@ out2=mass_density.dat
 out3=charge_density.dat
 out4=electron_density.dat
 
-CleanFiles $in $out1 $out2 $out3 $out4 total.dat tz2.wato.agr
+CleanFiles $in $out1 $out2 $out3 $out4 total.dat tz2.wato.agr restrict.agr
 
 INPUT="-i $in"
 TESTNAME='Density tests'
@@ -49,9 +49,11 @@ parm ../tz2.ortho.parm7
 trajin ../tz2.ortho.nc
 autoimage origin
 density :WAT@O out tz2.wato.agr xydy
+density :WAT@O out restrict.agr xydy delta 0.5 restrict cylinder cutoff 5.0
 EOF
   RunCpptraj "$UNITNAME"
   DoTest tz2.wato.agr.save tz2.wato.agr
+  DoTest restrict.agr.save restrict.agr
 fi
 
 # Total system density test