From aebcb15e19a24c75ea305cb648ca5eb3d789a83c Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Dec 2020 11:49:57 -0500 Subject: [PATCH 1/6] Add 'restrict' keyword to 'density' action to only bin atoms that fall within a square or cylinder from specified axis. --- src/Action_Density.cpp | 93 ++++++++++++++++++++++++++++++++++-------- src/Action_Density.h | 3 ++ 2 files changed, 80 insertions(+), 16 deletions(-) diff --git a/src/Action_Density.cpp b/src/Action_Density.cpp index 2c06d8be91..6b6cf73ba1 100644 --- a/src/Action_Density.cpp +++ b/src/Action_Density.cpp @@ -14,7 +14,9 @@ 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; } @@ -22,7 +24,8 @@ Action_Density::Action_Density() : void Action_Density::Help() const { mprintf("\t[out ] [name ]\n" "\t[ ... [delta ] [{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 ] ]\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::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 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 diff --git a/src/Action_Density.h b/src/Action_Density.h index cc5cf579e2..b769199e72 100644 --- a/src/Action_Density.h +++ b/src/Action_Density.h @@ -58,6 +58,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 +78,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 From 81c3b2dd0cc487f2d91d36a3abb36905cec6875f Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Dec 2020 11:52:03 -0500 Subject: [PATCH 2/6] Update version of density --- src/Action_Density.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Action_Density.h b/src/Action_Density.h index b769199e72..8a0eaa5352 100644 --- a/src/Action_Density.h +++ b/src/Action_Density.h @@ -10,7 +10,7 @@ #include "OnlineVarT.h" #include "ImagedAction.h" -#define ROUTINE_VERSION_STRING "1.0.2" +#define ROUTINE_VERSION_STRING "1.0.3" /** Calculate density along a coordinate. * \author Hannes H. Loeffler. From a8e4eb2521880bd12fbfa74e736c99f511fa9c15 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Dec 2020 11:52:55 -0500 Subject: [PATCH 3/6] Revision bump for adding restrict to density --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index da3e2d0f2a..85aaf9f98e 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 "V4.32.1" +#define CPPTRAJ_INTERNAL_VERSION "V4.32.2" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif From 78e825d5389e409988b3142992cff845e9c2b29b Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Dec 2020 11:57:35 -0500 Subject: [PATCH 4/6] Add test for restrict keyword --- test/Test_Density/RunTest.sh | 4 +- test/Test_Density/restrict.agr.save | 84 +++++++++++++++++++++++++++++ 2 files changed, 87 insertions(+), 1 deletion(-) create mode 100644 test/Test_Density/restrict.agr.save diff --git a/test/Test_Density/RunTest.sh b/test/Test_Density/RunTest.sh index 3f46c1a039..97ef5fd0d5 100755 --- a/test/Test_Density/RunTest.sh +++ b/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 diff --git a/test/Test_Density/restrict.agr.save b/test/Test_Density/restrict.agr.save new file mode 100644 index 0000000000..b97a602d77 --- /dev/null +++ b/test/Test_Density/restrict.agr.save @@ -0,0 +1,84 @@ +@with g0 +@ xaxis label "Z" +@ yaxis label "" +@ legend 0.2, 0.995 +@ legend char size 0.60 +@ s0 legend ":WAT@O" +@target G0.S0 +@type xydy + -18.750 0.0000 0.0000 + -18.250 1.4000 0.4830 + -17.750 2.2000 1.1005 + -17.250 2.8000 1.0750 + -16.750 3.0000 1.3540 + -16.250 2.4000 0.7888 + -15.750 2.8000 0.9661 + -15.250 1.2000 0.6992 + -14.750 2.4000 0.7888 + -14.250 2.6000 0.8233 + -13.750 1.8000 0.9944 + -13.250 4.6000 1.4944 + -12.750 2.2000 0.9944 + -12.250 2.8000 1.0750 + -11.750 2.4000 1.3984 + -11.250 1.6000 0.7888 + -10.750 2.8000 0.8433 + -10.250 2.6000 0.9487 + -9.750 1.2000 0.9661 + -9.250 4.0000 1.6997 + -8.750 2.8000 1.0750 + -8.250 2.4000 0.6325 + -7.750 1.8000 0.8756 + -7.250 1.8000 0.8756 + -6.750 1.4000 0.8233 + -6.250 1.6000 0.4216 + -5.750 0.6000 0.6749 + -5.250 0.8000 0.5164 + -4.750 2.0000 0.8165 + -4.250 1.2000 0.6992 + -3.750 0.4000 0.4216 + -3.250 0.2000 0.3162 + -2.750 0.2000 0.3162 + -2.250 0.2000 0.3162 + -1.750 0.0000 0.0000 + -1.250 0.0000 0.0000 + -0.750 0.0000 0.0000 + -0.250 0.0000 0.0000 + 0.250 0.0000 0.0000 + 0.750 0.0000 0.0000 + 1.250 0.0000 0.0000 + 1.750 0.0000 0.0000 + 2.250 0.0000 0.0000 + 2.750 0.4000 0.4216 + 3.250 0.0000 0.0000 + 3.750 1.0000 0.5270 + 4.250 0.2000 0.3162 + 4.750 0.2000 0.3162 + 5.250 1.0000 0.8498 + 5.750 1.4000 0.6749 + 6.250 1.4000 0.6749 + 6.750 1.0000 0.7071 + 7.250 1.8000 0.7379 + 7.750 1.4000 0.6749 + 8.250 1.6000 1.2293 + 8.750 2.6000 1.3375 + 9.250 2.8000 0.9661 + 9.750 2.2000 1.1005 + 10.250 3.2000 1.2649 + 10.750 2.2000 0.8756 + 11.250 3.0000 1.0801 + 11.750 3.2000 0.8433 + 12.250 2.6000 1.1595 + 12.750 3.2000 1.1738 + 13.250 2.0000 1.0541 + 13.750 2.2000 0.9944 + 14.250 1.8000 0.7379 + 14.750 1.6000 0.6325 + 15.250 2.4000 0.7888 + 15.750 2.8000 0.9661 + 16.250 3.0000 1.0801 + 16.750 2.4000 0.9189 + 17.250 2.2000 0.9944 + 17.750 2.4000 0.7888 + 18.250 0.2000 0.3162 + 18.750 0.0000 0.0000 From 9c8a1da13a4606c7e70a024837327523141eefc4 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Dec 2020 18:54:02 -0500 Subject: [PATCH 5/6] Update manual with 'restrict' keyword. --- doc/cpptraj.lyx | 79 +++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 70 insertions(+), 9 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 87b2407735..acd243591f 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -21500,29 +21500,33 @@ density \end_layout \begin_layout LyX-Code -density [out ] +density [out ] [name ] \end_layout \begin_layout LyX-Code [ ... - [name ] [delta ] [{x|y|z}] + [delta ] [{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 ] ] +\end_layout + \begin_deeper \begin_layout Description [out \begin_inset space ~ \end_inset -] Output file for histogram (relative distances vs. +] 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 + + 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. +[avg]: Average density over coordinate for mask number . +\end_layout + +\begin_layout Description +[sd]: Standard deviation of density over coordinate for mask number + . +\end_layout + +\begin_layout Standard +DataSets Created if no masks specified: \end_layout \begin_layout Description -[sd] Standard deviation of density over coordinate. + 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 From 60ea7189bd8e343a51df271d71f01cdd115abdbc Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Tue, 15 Dec 2020 19:01:50 -0500 Subject: [PATCH 6/6] Ensure both cutoffs satisfied for SQUARE in order to bin. Add note about version --- src/Action_Density.cpp | 2 +- src/Action_Density.h | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Action_Density.cpp b/src/Action_Density.cpp index 6b6cf73ba1..4e13e49edb 100644 --- a/src/Action_Density.cpp +++ b/src/Action_Density.cpp @@ -273,7 +273,7 @@ Action::RetType Action_Density::HistAction(int frameNum, ActionFrame& frm) { binValues = true; //mprintf("DEBUG: binValues= %i\n", (int)binValues); } else if (restrictType_ == SQUARE) { - if (fabs(XYZ[area_coord_[0]]) < cutVal_ || + if (fabs(XYZ[area_coord_[0]]) < cutVal_ && fabs(XYZ[area_coord_[1]]) < cutVal_) binValues = true; } else diff --git a/src/Action_Density.h b/src/Action_Density.h index 8a0eaa5352..d5be443bec 100644 --- a/src/Action_Density.h +++ b/src/Action_Density.h @@ -11,6 +11,7 @@ #include "ImagedAction.h" #define ROUTINE_VERSION_STRING "1.0.3" +// 1.0.3 - DRR Add 'restrict' /** Calculate density along a coordinate. * \author Hannes H. Loeffler.