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 diff --git a/src/Action_Mask.cpp b/src/Action_Mask.cpp index 896c39dcc9..61854ff32a 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), @@ -17,13 +18,17 @@ Action_Mask::Action_Mask() : {} void Action_Mask::Help() const { - mprintf("\t [maskout ]\n" - "\t[ {maskpdb | maskmol2 }\n" - "\t [trajargs ] ]\n" + mprintf("\t [maskout ] [out ] [nselectedout ]\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" - " 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)) { 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 diff --git a/src/RPNcalc.cpp b/src/RPNcalc.cpp index d591a6022a..f5ef18931a 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) { @@ -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( *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()); @@ -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( *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 ) { 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 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 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 <