From 5e789009ff971420355653fbd9c0dc287e1288b0 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 25 Sep 2020 09:38:33 -0400 Subject: [PATCH 1/8] For DataSet OP DataSet, use X values from first DataSet --- src/RPNcalc.cpp | 28 +++++++++++++++++++++------- src/cpptrajdepend | 2 +- 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/src/RPNcalc.cpp b/src/RPNcalc.cpp index d591a6022a..fd72667a87 100644 --- a/src/RPNcalc.cpp +++ b/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( *tempDS ); - D0.Allocate( DataSet::SizeArray(1, ds1->Size()) ); - DataSet_1D const& D1 = static_cast( *ds1 ); - DataSet_1D const& D2 = static_cast( *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( *tempDS ); + D0.Allocate( DataSet::SizeArray(1, ds1->Size()) ); + DataSet_1D const& D1 = static_cast( *ds1 ); + DataSet_1D const& D2 = static_cast( *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( *tempDS ); + D0.Allocate( DataSet::SizeArray(1, ds1->Size()) ); + DataSet_1D const& D1 = static_cast( *ds1 ); + DataSet_1D const& D2 = static_cast( *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) { diff --git a/src/cpptrajdepend b/src/cpptrajdepend index 8447092034..58f05e58d1 100644 --- a/src/cpptrajdepend +++ b/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 From 469db2d02589d3a50b468631ac436ed4f88a3d64 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 25 Sep 2020 09:47:26 -0400 Subject: [PATCH 2/8] Preserve 1D X values for dataset OP value and value OP dataset --- src/RPNcalc.cpp | 49 +++++++++++++++++++++++++++++++++++++------------ 1 file changed, 37 insertions(+), 12 deletions(-) diff --git a/src/RPNcalc.cpp b/src/RPNcalc.cpp index fd72667a87..f5ef18931a 100644 --- a/src/RPNcalc.cpp +++ b/src/RPNcalc.cpp @@ -820,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( *tempDS ); - D0.Allocate( DataSet::SizeArray(1,ds2->Size()) ); - DataSet_1D const& D2 = static_cast( *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( *tempDS ); + D0.Allocate( DataSet::SizeArray(1, ds2->Size()) ); + DataSet_1D const& D2 = static_cast( *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( *tempDS ); + D0.Allocate( DataSet::SizeArray(1,ds2->Size()) ); + DataSet_1D const& D2 = static_cast( *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()); @@ -840,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( *tempDS ); - D0.Allocate( DataSet::SizeArray(1, ds1->Size()) ); - DataSet_1D const& D1 = static_cast( *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( *tempDS ); + D0.Allocate( DataSet::SizeArray(1, ds1->Size()) ); + DataSet_1D const& D1 = static_cast( *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( *tempDS ); + D0.Allocate( DataSet::SizeArray(1, ds1->Size()) ); + DataSet_1D const& D1 = static_cast( *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 ) { From 09ffa0e6117d9543edcfeba640b38accd44a216b Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 25 Sep 2020 10:59:57 -0400 Subject: [PATCH 3/8] Add dataset to mask command to hold number of atoms selected per frame --- src/Action_Mask.cpp | 10 ++++++++-- src/Action_Mask.h | 1 + 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/Action_Mask.cpp b/src/Action_Mask.cpp index 896c39dcc9..128df40ff1 100644 --- a/src/Action_Mask.cpp +++ b/src/Action_Mask.cpp @@ -4,6 +4,7 @@ // CONSTRUCTOR Action_Mask::Action_Mask() : outfile_(0), + nselected_(0), fnum_(0), anum_(0), aname_(0), @@ -82,9 +83,12 @@ 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; + // 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 +145,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)) { diff --git a/src/Action_Mask.h b/src/Action_Mask.h index ce737b8cba..ec81c7b149 100644 --- a/src/Action_Mask.h +++ b/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 From c8d134fd600f3b6c1194609a83261cfa804a0d8d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 25 Sep 2020 11:10:15 -0400 Subject: [PATCH 4/8] Add nselectedout for mask command --- src/Action_Mask.cpp | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/src/Action_Mask.cpp b/src/Action_Mask.cpp index 128df40ff1..c7a6c49f51 100644 --- a/src/Action_Mask.cpp +++ b/src/Action_Mask.cpp @@ -18,13 +18,17 @@ Action_Mask::Action_Mask() : {} void Action_Mask::Help() const { - mprintf("\t [maskout ]\n" + mprintf("\t [maskout ] [out ] [nselectedout ]\n" "\t[ {maskpdb | maskmol2 }\n" "\t [trajargs ] ]\n" " Print atoms selected by 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() @@ -47,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()) { @@ -87,6 +87,11 @@ Action::RetType Action_Mask::Init(ArgList& actionArgs, ActionInit& init, int deb 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()) { MetaData::tsType ts = MetaData::NOT_TS; // None are a straight time series From 26c5344fcbd350c85a13ab530e9e8a603c741269 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 25 Sep 2020 11:14:52 -0400 Subject: [PATCH 5/8] Add name to mask help --- src/Action_Mask.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Action_Mask.cpp b/src/Action_Mask.cpp index c7a6c49f51..61854ff32a 100644 --- a/src/Action_Mask.cpp +++ b/src/Action_Mask.cpp @@ -19,8 +19,8 @@ Action_Mask::Action_Mask() : void Action_Mask::Help() const { mprintf("\t [maskout ] [out ] [nselectedout ]\n" - "\t[ {maskpdb | maskmol2 }\n" - "\t [trajargs ] ]\n" + "\t[name ] [ {maskpdb | maskmol2 }\n" + "\t [trajargs ] ]\n" " Print atoms selected by 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" From 79b3674b2e083d133dfac9867eece9ee52bf8d6a Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 25 Sep 2020 11:17:45 -0400 Subject: [PATCH 6/8] Fix up mask manual entry. --- doc/cpptraj.lyx | 63 ++++++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 60 insertions(+), 3 deletions(-) diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index ab71a0d6a5..34c535d1ae 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -27806,15 +27806,15 @@ mask \end_inset - [maskout ] + [maskout ] [out ] [nselectedout ] \end_layout \begin_layout LyX-Code - [ {maskpdb | maskmol2 } + [name ] [ {maskpdb | maskmol2 } \end_layout \begin_layout LyX-Code - [trajargs ] ] + [trajargs ] ] \end_layout \begin_deeper @@ -27830,6 +27830,31 @@ maskout Write information on atoms in to . \end_layout +\begin_layout Description +out +\begin_inset space ~ +\end_inset + + 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 + + Write the total number of selected atoms to file. +\end_layout + +\begin_layout Description +name +\begin_inset space ~ +\end_inset + + 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 + Number of atoms selected each frame. +\end_layout + +\begin_layout Description +[Frm] Frame number for each selected atom. +\end_layout + +\begin_layout Description +[AtNum] Atom number for each selected atom. +\end_layout + +\begin_layout Description +[Aname] Atom name for each selected atom. +\end_layout + +\begin_layout Description +[Rnum] Residue number for each selected atom. +\end_layout + +\begin_layout Description +[Rname] Residue name for each selected atom. +\end_layout + +\begin_layout Description +[Mnum] Molecule number for each selected atom. +\end_layout + \end_deeper \begin_layout Standard For each frame determine all atoms that correspond to From f7cbe15f99cb9f8909fafe1e4f63e6fd36efcc56 Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 25 Sep 2020 11:18:21 -0400 Subject: [PATCH 7/8] Revision bump for fixed handling of dataset X values during calc and mask command nselectedout option. --- src/Version.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Version.h b/src/Version.h index 7b95b0741d..db6e9b0408 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.30.1" +#define CPPTRAJ_INTERNAL_VERSION "V4.30.2" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif From 9005fc0000ea1b10dc5d32751e65da952b90235d Mon Sep 17 00:00:00 2001 From: "Daniel R. Roe" Date: Fri, 25 Sep 2020 11:21:10 -0400 Subject: [PATCH 8/8] Add test for the nselectedout option --- test/Test_Mask/Nselected.dat.save | 11 +++++++++++ test/Test_Mask/RunTest.sh | 5 +++-- 2 files changed, 14 insertions(+), 2 deletions(-) create mode 100644 test/Test_Mask/Nselected.dat.save diff --git a/test/Test_Mask/Nselected.dat.save b/test/Test_Mask/Nselected.dat.save new file mode 100644 index 0000000000..b214a60735 --- /dev/null +++ b/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 diff --git a/test/Test_Mask/RunTest.sh b/test/Test_Mask/RunTest.sh index 65830b6cf5..04f331405f 100755 --- a/test/Test_Mask/RunTest.sh +++ b/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 <