Fix X values in data set calc. Add 'nselectedout' option to mask. (#848)

* For DataSet OP DataSet, use X values from first DataSet

* Preserve 1D X values for dataset OP value and value OP dataset

* Add dataset to mask command to hold number of atoms selected per frame

* Add nselectedout for mask command

* Add name to mask help

* Fix up mask manual entry.

* Revision bump for fixed handling of dataset X values during calc and mask command nselectedout option.

* Add test for the nselectedout option

doc/cpptraj.lyx

@@ -27806,15 +27806,15 @@ mask
 
 \end_inset
 
- <mask> [maskout <filename>]
+ <mask> [maskout <filename>] [out <filename>] [nselectedout <filename>]
 \end_layout
 
 \begin_layout LyX-Code
-     [ {maskpdb <filename> | maskmol2 <filename>}
+     [name <setname>] [ {maskpdb <filename> | maskmol2 <filename>}
 \end_layout
 
 \begin_layout LyX-Code
-       [trajargs <comma-separated args>] ]
+                        [trajargs <comma-separated args>] ]
 \end_layout
 
 \begin_deeper
@@ -27830,6 +27830,31 @@ maskout
 <filename> Write information on atoms in <mask> to <filename>.
 \end_layout
 
+\begin_layout Description
+out
+\begin_inset space ~
+\end_inset
+
+<filename> Write the frame, atom number, atom name, residue number, residue
+ name, and molecule number for each selected atom to file.
+\end_layout
+
+\begin_layout Description
+nselectedout
+\begin_inset space ~
+\end_inset
+
+<filename> Write the total number of selected atoms to file.
+\end_layout
+
+\begin_layout Description
+name
+\begin_inset space ~
+\end_inset
+
+<setname> Name for output data sets.
+\end_layout
+
 \begin_layout Description
 maskpdb
 \begin_inset space ~
@@ -27859,6 +27884,38 @@ args> When writing output PDB/Mol2, additional trajectory arguments to pass
  to the output trajectory.
 \end_layout
 
+\begin_layout Standard
+DataSets Created
+\end_layout
+
+\begin_layout Description
+<name> Number of atoms selected each frame.
+\end_layout
+
+\begin_layout Description
+<name>[Frm] Frame number for each selected atom.
+\end_layout
+
+\begin_layout Description
+<name>[AtNum] Atom number for each selected atom.
+\end_layout
+
+\begin_layout Description
+<name>[Aname] Atom name for each selected atom.
+\end_layout
+
+\begin_layout Description
+<name>[Rnum] Residue number for each selected atom.
+\end_layout
+
+\begin_layout Description
+<name>[Rname] Residue name for each selected atom.
+\end_layout
+
+\begin_layout Description
+<name>[Mnum] Molecule number for each selected atom.
+\end_layout
+
 \end_deeper
 \begin_layout Standard
 For each frame determine all atoms that correspond to 

src/Action_Mask.cpp

@@ -4,6 +4,7 @@
 // CONSTRUCTOR
 Action_Mask::Action_Mask() :
   outfile_(0),
+  nselected_(0),
   fnum_(0),
   anum_(0),
   aname_(0),
@@ -17,13 +18,17 @@ Action_Mask::Action_Mask() :
 {} 
 
 void Action_Mask::Help() const {
-  mprintf("\t<mask1> [maskout <filename>]\n"
-          "\t[ {maskpdb <filename> | maskmol2 <filename>}\n"
-          "\t  [trajargs <comma-separated args>] ]\n"
+  mprintf("\t<mask1> [maskout <filename>] [out <filename>] [nselectedout <filename>]\n"
+          "\t[name <setname>] [ {maskpdb <filename> | maskmol2 <filename>}\n"
+          "\t                   [trajargs <comma-separated args>] ]\n"
           "  Print atoms selected by <mask1> to file specified by 'maskout' and/or\n"
           "  the PDB or Mol2 file specified by 'maskpdb' or 'maskmol2'. Additional\n"
           "  trajectory arguments can be specified in a comma-separated list via\n"
-          "  the 'trajargs' keyword. Good for distance-based masks.\n");
+          "  the 'trajargs' keyword. Good for distance-based masks.\n"
+          "  If 'out' is specified, the file will contain for each selected atom\n"
+          "  the frame, atom number, atom name, residue number, residue name, and\n"
+          "  molecule number. If 'nselectedout' is specified, the file will contain\n"
+          "  the total number of atoms selected each frame.\n");
 }
 
 // Action_Mask::Init()
@@ -46,12 +51,8 @@ Action::RetType Action_Mask::Init(ArgList& actionArgs, ActionInit& init, int deb
   std::string maskmol2 = actionArgs.GetStringKey("maskmol2");
   std::string dsname = actionArgs.GetStringKey("name");
   std::string dsout = actionArgs.GetStringKey("out");
+  std::string nSelectedArg = actionArgs.GetStringKey("nselectedout");
   std::string additionalTrajArgs = actionArgs.GetStringKey("trajargs");
-  // At least 1 of maskout, maskpdb, maskmol2, or name must be specified.
-  if (outfile_ == 0 && maskpdb.empty() && maskmol2.empty() && dsname.empty()) {
-    mprinterr("Error: At least one of maskout, maskpdb, maskmol2, or name must be specified.\n");
-    return Action::ERR;
-  }
   // Set up any trajectory options
   TrajectoryFile::TrajFormatType trajFmt = TrajectoryFile::PDBFILE;
   if (!maskpdb.empty() || !maskmol2.empty()) {
@@ -82,9 +83,17 @@ Action::RetType Action_Mask::Init(ArgList& actionArgs, ActionInit& init, int deb
     writeTraj_ = false;
   // Get Mask
   if (Mask1_.SetMaskString( actionArgs.GetMaskNext() )) return Action::ERR;
-  // Set up data sets
+  // Set up Nselected data set
+  if (dsname.empty()) dsname = init.DSL().GenerateDefaultName("MASK");
+  nselected_ = init.DSL().AddSet(DataSet::INTEGER, MetaData(dsname));
+  if (nselected_ == 0) return Action::ERR;
+  DataFile* nSelectedOut = 0;
+  if (!nSelectedArg.empty()) {
+    nSelectedOut = init.DFL().AddDataFile( nSelectedArg, actionArgs );
+    nSelectedOut->AddDataSet( nselected_ );
+  }
+  // Set up additional data sets
   if (!dsname.empty() || !dsout.empty()) {
-    if (dsname.empty()) dsname = init.DSL().GenerateDefaultName("MASK");
     MetaData::tsType ts = MetaData::NOT_TS; // None are a straight time series
     fnum_  = init.DSL().AddSet(DataSet::INTEGER, MetaData(dsname, "Frm",   ts));
     anum_  = init.DSL().AddSet(DataSet::INTEGER, MetaData(dsname, "AtNum", ts));
@@ -141,6 +150,8 @@ Action::RetType Action_Mask::DoAction(int frameNum, ActionFrame& frm) {
             Mask1_.MaskString());
     return Action::ERR;
   }
+  int nAt = Mask1_.Nselected();
+  nselected_->Add(frameNum, &nAt);
   // Print out information for every atom in the mask
   for (int atom=0; atom < CurrentParm_->Natom(); atom++) {
     if (Mask1_.AtomInCharMask(atom)) {

src/Action_Mask.h

@@ -25,6 +25,7 @@ class Action_Mask: public Action {
 
     CharMask Mask1_;         ///< Atoms which will be selected each frame
     CpptrajFile* outfile_;   ///< File to write selected atom info to
+    DataSet* nselected_;     ///< Hold number of atoms selected per frame
     DataSet* fnum_;          ///< Hold frame numbers for selections
     DataSet* anum_;          ///< Hold selected atom numbers
     DataSet* aname_;         ///< Hold selected atom names

src/RPNcalc.cpp

@@ -6,6 +6,7 @@
 #include "DataSetList.h"
 #include "DataSet_Vector.h"
 #include "DataSet_double.h"
+#include "DataSet_Mesh.h"
 #include "DataSet_MatrixDbl.h"
 #include "DataSet_GridFlt.h"
 #include "CpptrajStdio.h"
@@ -679,13 +680,26 @@ int RPNcalc::TokenLoop(DataSetList& DSL) const {
                           ds1->legend(), ds2->legend(), T->name());
                 return 1;
               }
-              tempDS = LocalList.AddSet(DataSet::DOUBLE, MetaData("TEMP", T-tokens_.begin()));
-              DataSet_double& D0 = static_cast<DataSet_double&>( *tempDS );
-              D0.Allocate( DataSet::SizeArray(1, ds1->Size()) );
-              DataSet_1D const& D1 = static_cast<DataSet_1D const&>( *ds1 );
-              DataSet_1D const& D2 = static_cast<DataSet_1D const&>( *ds2 );
-              for (unsigned int n = 0; n != D1.Size(); n++)
-                D0.AddElement( DoOperation(D1.Dval(n), D2.Dval(n), T->Type()) );
+              // Determine how X values will be handled. Default to using X values from D1.
+              // For mesh need to AddXY. For all others just set the dimension.
+              if (ds1->Type() == DataSet::XYMESH) {
+                tempDS = LocalList.AddSet(DataSet::XYMESH, MetaData("TEMP", T-tokens_.begin()));
+                DataSet_Mesh& D0 = static_cast<DataSet_Mesh&>( *tempDS );
+                D0.Allocate( DataSet::SizeArray(1, ds1->Size()) );
+                DataSet_1D const& D1 = static_cast<DataSet_1D const&>( *ds1 );
+                DataSet_1D const& D2 = static_cast<DataSet_1D const&>( *ds2 );
+                for (unsigned int n = 0; n != D1.Size(); n++)
+                  D0.AddXY( D1.Xcrd(n), DoOperation(D1.Dval(n), D2.Dval(n), T->Type()) );
+              } else {
+                tempDS = LocalList.AddSet(DataSet::DOUBLE, MetaData("TEMP", T-tokens_.begin()));
+                tempDS->SetDim(Dimension::X, ds1->Dim(0));
+                DataSet_double& D0 = static_cast<DataSet_double&>( *tempDS );
+                D0.Allocate( DataSet::SizeArray(1, ds1->Size()) );
+                DataSet_1D const& D1 = static_cast<DataSet_1D const&>( *ds1 );
+                DataSet_1D const& D2 = static_cast<DataSet_1D const&>( *ds2 );
+                for (unsigned int n = 0; n != D1.Size(); n++)
+                  D0.AddElement( DoOperation(D1.Dval(n), D2.Dval(n), T->Type()) );
+              }
             }
             else if (ds1->Group() == DataSet::VECTOR_1D && ds2->Group() == DataSet::VECTOR_1D)
             {
@@ -806,12 +820,25 @@ int RPNcalc::TokenLoop(DataSetList& DSL) const {
                         ds2->legend(), T-tokens_.begin());
               if (ScalarTimeSeries( ds2 ))
               {
-                tempDS = LocalList.AddSet(DataSet::DOUBLE, MetaData("TEMP", T-tokens_.begin()));
-                DataSet_double& D0 = static_cast<DataSet_double&>( *tempDS );
-                D0.Allocate( DataSet::SizeArray(1,ds2->Size()) );
-                DataSet_1D const& D2 = static_cast<DataSet_1D const&>( *ds2 );
-                for (unsigned int n = 0; n != D2.Size(); n++)
-                  D0.AddElement( DoOperation(D2.Dval(n), d1, T->Type()) );
+                // Determine how X values will be handled. Default to using X values from D2.
+                // For mesh need to AddXY. For all others just set the dimension.
+                // TODO does the order of D2 and d1 matter here?
+                if ( ds2->Type() == DataSet::XYMESH ) {
+                  tempDS = LocalList.AddSet(DataSet::XYMESH, MetaData("TEMP", T-tokens_.begin()));
+                  DataSet_Mesh& D0 = static_cast<DataSet_Mesh&>( *tempDS );
+                  D0.Allocate( DataSet::SizeArray(1, ds2->Size()) );
+                  DataSet_1D const& D2 = static_cast<DataSet_1D const&>( *ds2 );
+                  for (unsigned int n = 0; n != D2.Size(); n++)
+                    D0.AddXY( D2.Xcrd(n), DoOperation(D2.Dval(n), d1, T->Type()) );
+                } else {
+                  tempDS = LocalList.AddSet(DataSet::DOUBLE, MetaData("TEMP", T-tokens_.begin()));
+                  tempDS->SetDim(Dimension::X, ds2->Dim(0));
+                  DataSet_double& D0 = static_cast<DataSet_double&>( *tempDS );
+                  D0.Allocate( DataSet::SizeArray(1,ds2->Size()) );
+                  DataSet_1D const& D2 = static_cast<DataSet_1D const&>( *ds2 );
+                  for (unsigned int n = 0; n != D2.Size(); n++)
+                    D0.AddElement( DoOperation(D2.Dval(n), d1, T->Type()) );
+                }
               } else {
                 mprinterr("Error: Operation '%s' between value and set %s not yet permitted.\n",
                           T->Description(), ds2->legend());
@@ -826,12 +853,24 @@ int RPNcalc::TokenLoop(DataSetList& DSL) const {
                         d2, T-tokens_.begin());
               if (ScalarTimeSeries( ds1 ))
               {
-                tempDS = LocalList.AddSet(DataSet::DOUBLE, MetaData("TEMP", T-tokens_.begin()));
-                DataSet_double& D0 = static_cast<DataSet_double&>( *tempDS );
-                D0.Allocate( DataSet::SizeArray(1, ds1->Size()) );
-                DataSet_1D const& D1 = static_cast<DataSet_1D const&>( *ds1 );
-                for (unsigned int n = 0; n != D1.Size(); n++)
-                  D0.AddElement( DoOperation(d2, D1.Dval(n), T->Type()) );
+                // Determine how X values will be handled. Default to using X values from D1.
+                // For mesh need to AddXY. For all others just set the dimension.
+                if ( ds1->Type() == DataSet::XYMESH ) {
+                  tempDS = LocalList.AddSet(DataSet::XYMESH, MetaData("TEMP", T-tokens_.begin()));
+                  DataSet_Mesh& D0 = static_cast<DataSet_Mesh&>( *tempDS );
+                  D0.Allocate( DataSet::SizeArray(1, ds1->Size()) );
+                  DataSet_1D const& D1 = static_cast<DataSet_1D const&>( *ds1 );
+                  for (unsigned int n = 0; n != D1.Size(); n++)
+                    D0.AddXY( D1.Xcrd(n), DoOperation(d2, D1.Dval(n), T->Type()) );
+                } else {
+                  tempDS = LocalList.AddSet(DataSet::DOUBLE, MetaData("TEMP", T-tokens_.begin()));
+                  tempDS->SetDim(Dimension::X, ds1->Dim(0));
+                  DataSet_double& D0 = static_cast<DataSet_double&>( *tempDS );
+                  D0.Allocate( DataSet::SizeArray(1, ds1->Size()) );
+                  DataSet_1D const& D1 = static_cast<DataSet_1D const&>( *ds1 );
+                  for (unsigned int n = 0; n != D1.Size(); n++)
+                    D0.AddElement( DoOperation(d2, D1.Dval(n), T->Type()) );
+                }
               }
               else if ( ds1->Group() == DataSet::VECTOR_1D )
               {

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

src/cpptrajdepend

@@ -348,7 +348,7 @@ PotentialTerm_LJ_Coulomb.o : PotentialTerm_LJ_Coulomb.cpp Atom.h AtomExtra.h Ato
 ProgressBar.o : ProgressBar.cpp CpptrajStdio.h ProgressBar.h
 ProgressTimer.o : ProgressTimer.cpp CpptrajStdio.h ProgressTimer.h Timer.h
 PubFFT.o : PubFFT.cpp ArrayIterator.h ComplexArray.h CpptrajStdio.h PubFFT.h
-RPNcalc.o : RPNcalc.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h Box.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_MatrixDbl.h DataSet_Vector.h DataSet_double.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h RPNcalc.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
+RPNcalc.o : RPNcalc.cpp ArgList.h ArrayIterator.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h Box.h ComplexArray.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajStdio.h DataSet.h DataSetList.h DataSet_1D.h DataSet_2D.h DataSet_3D.h DataSet_Coords.h DataSet_Coords_REF.h DataSet_GridFlt.h DataSet_MatrixDbl.h DataSet_Mesh.h DataSet_Vector.h DataSet_double.h Dimension.h FileIO.h FileName.h Frame.h Grid.h GridBin.h MaskToken.h Matrix.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h RPNcalc.h Range.h ReferenceFrame.h ReplicaDimArray.h Residue.h Segment.h Spline.h StringRoutines.h SymbolExporting.h TextFormat.h Timer.h Topology.h TypeNameHolder.h Unit.h Vec3.h
 Random.o : Random.cpp CpptrajStdio.h Random.h
 Range.o : Range.cpp ArgList.h CpptrajStdio.h Range.h
 ReadLine.o : ReadLine.cpp Action.h ActionList.h ActionState.h Analysis.h AnalysisList.h AnalysisState.h ArgList.h AssociatedData.h Atom.h AtomExtra.h AtomMask.h AtomType.h BaseIOtype.h Box.h Cmd.h CmdInput.h CmdList.h Command.h Constants.h CoordinateInfo.h CpptrajFile.h CpptrajState.h CpptrajStdio.h DataFile.h DataFileList.h DataSet.h DataSetList.h DataSet_Coords.h DataSet_Coords_REF.h Dimension.h DispatchObject.h EnsembleIn.h EnsembleOutList.h FileIO.h FileName.h FileTypes.h Frame.h FramePtrArray.h InputTrajCommon.h MaskToken.h Matrix_3x3.h MetaData.h Molecule.h NameType.h Parallel.h ParameterHolders.h ParameterSet.h ParameterTypes.h Range.h ReadLine.h ReferenceFrame.h ReplicaDimArray.h ReplicaInfo.h Residue.h Segment.h SymbolExporting.h TextFormat.h Timer.h Topology.h TrajFrameCounter.h Trajin.h TrajinList.h TrajoutList.h TypeNameHolder.h Unit.h Vec3.h

test/Test_Mask/Nselected.dat.save

@@ -0,0 +1,11 @@
+#Frame              M
+       1            3
+       2            2
+       3            4
+       4            3
+       5            3
+       6            2
+       7            4
+       8            3
+       9            2
+      10            3

test/Test_Mask/RunTest.sh

@@ -3,7 +3,7 @@
 . ../MasterTest.sh
 
 # Clean
-CleanFiles mask.in mask.out mask.pdb.1 mask.mol2.1 M.dat
+CleanFiles mask.in mask.out mask.pdb.1 mask.mol2.1 M.dat Nselected.dat
 
 TESTNAME='Mask command tests'
 Requires netcdf maxthreads 10
@@ -30,10 +30,11 @@ cat > mask.in <<EOF
 noprogress
 parm ../tz2.ortho.parm7
 trajin ../tz2.ortho.nc
-mask "(:8@NZ <:3.0) & :WAT@O" name M out M.dat
+mask "(:8@NZ <:3.0) & :WAT@O" name M out M.dat nselectedout Nselected.dat
 EOF
 RunCpptraj "Mask longer command test."
 DoTest M.dat.save M.dat
+DoTest Nselected.dat.save Nselected.dat
 
 EndTest