diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf index 907e1b7826..3b7084f674 100644 Binary files a/doc/CpptrajManual.pdf and b/doc/CpptrajManual.pdf differ diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx index 52f18bb580..b5c516d5d1 100644 --- a/doc/cpptraj.lyx +++ b/doc/cpptraj.lyx @@ -18812,19 +18812,15 @@ closest \end_inset - <# to keep> [noimage] [first | oxygen] [center] + <# to keep> [solventmask ] [noimage] \end_layout \begin_layout LyX-Code - [closestout ] [name ] [outprefix ] -\begin_inset Separator latexpar -\end_inset - - + [first | oxygen] [center] [closestout [name ]] \end_layout \begin_layout LyX-Code - [parmout ] + [outprefix ] [parmout ] \end_layout \begin_deeper @@ -18844,6 +18840,27 @@ keep> Number of solvent molecules to keep around Mask of atoms to search for closest waters around. \end_layout +\begin_layout Description +[solventmask +\begin_inset space ~ +\end_inset + +] Optional mask for selecting solvent atoms. + If not specified, atoms in all molecules marked as +\begin_inset Quotes eld +\end_inset + +solvent +\begin_inset Quotes erd +\end_inset + + will be used. +\end_layout + \begin_layout Description [noimage] Do not perform imaging; only recommended if trajectory has previously been imaged. @@ -18867,7 +18884,7 @@ oxygen] Calculate distances between all atoms in \end_layout \begin_layout Description -[center] Search for waters closest to center of +[center] Search for waters closest to geometric center of \series bold \series default @@ -18945,8 +18962,8 @@ strip \shape italic cpptraj \shape default - (by default residues named WAT, HOH, or TIP3) or can be specified prior - via the + (by default residues named WAT, HOH, or TIP3), can be specified prior via + the \series bold \shape italic solvent @@ -18959,8 +18976,15 @@ reference "subsec:cpptraj_solvent" \end_inset -). - The format of the +), or can be selected by +\series bold +solventmask +\series default +. +\end_layout + +\begin_layout Standard +The format of the \series bold closestout \series default diff --git a/src/Action_Closest.cpp b/src/Action_Closest.cpp index 4bd37890f9..91689b5a1e 100644 --- a/src/Action_Closest.cpp +++ b/src/Action_Closest.cpp @@ -39,9 +39,9 @@ Action_Closest::Action_Closest() : {} void Action_Closest::Help() const { - mprintf("\t<# to keep> [noimage] [first | oxygen] [center]\n" - "\t[closestout [name ]] [outprefix ]\n" - "\t[parmout ]\n" + mprintf("\t<# to keep> [solventmask ] [noimage]\n" + "\t[first | oxygen] [center] [closestout [name ]]\n" + "\t[outprefix ] [parmout ]\n" " Keep only the closest <# to keep> solvent molecules to atoms in .\n" " Molecules can be marked as solvent with the 'solvent' command.\n" " If 'center' specified use geometric center of atoms in .\n"); @@ -107,7 +107,10 @@ Action::RetType Action_Closest::Init(ArgList& actionArgs, ActionInit& init, int } // Get Masks - std::string mask1 = actionArgs.GetMaskNext(); + std::string mask1 = actionArgs.GetStringKey("solventmask"); + if (!mask1.empty()) + solventMask_.SetMaskString( mask1 ); + mask1 = actionArgs.GetMaskNext(); if (mask1.empty()) { mprinterr("Error: No mask specified.\n"); return Action::ERR; @@ -120,6 +123,8 @@ Action::RetType Action_Closest::Init(ArgList& actionArgs, ActionInit& init, int mprintf("\tGeometric center of atoms in mask will be used.\n"); if (!image_.UseImage()) mprintf("\tImaging will be turned off.\n"); + if (solventMask_.MaskStringSet()) + mprintf("\tSolvent will be selected by mask '%s'\n", solventMask_.MaskString()); if (firstAtom_) mprintf("\tOnly first atom of solvent molecule used for distance calc.\n"); if (outFile_!=0) @@ -140,20 +145,55 @@ Action::RetType Action_Closest::Init(ArgList& actionArgs, ActionInit& init, int * parm file. Atom masks for each solvent molecule will be set up. */ Action::RetType Action_Closest::Setup(ActionSetup& setup) { + if (setup.Top().Nmol() < 1) { + mprintf("Warning: 'closest' requires molecule information. Skipping.\n"); + return Action::SKIP; + } closestWaters_ = targetNclosest_; + // Determine solvent molecules + std::vector IsSolventMol; + IsSolventMol.reserve( setup.Top().Nmol() ); + int nSolvent = 0; + if (solventMask_.MaskStringSet()) { + // Use only atoms selected by solventMask + if (setup.Top().SetupCharMask( solventMask_ )) return Action::ERR; + solventMask_.MaskInfo(); + if (solventMask_.None()) { + mprintf("Warning: No solvent selected by mask '%s'. Skipping.\n", solventMask_.MaskString()); + return Action::SKIP; + } + for (Topology::mol_iterator Mol = setup.Top().MolStart(); + Mol != setup.Top().MolEnd(); ++Mol) + { + if ( solventMask_.AtomsInCharMask(Mol->BeginAtom(), Mol->EndAtom()) ) { + IsSolventMol.push_back( true ); + nSolvent++; + } else + IsSolventMol.push_back( false ); + } + } else { + // Select everything; all solvent atoms are fair game. + solventMask_.ResetMask(); + solventMask_.InitCharMask(setup.Top().Natom(), true); + for (Topology::mol_iterator Mol = setup.Top().MolStart(); + Mol != setup.Top().MolEnd(); ++Mol) + IsSolventMol.push_back( Mol->IsSolvent() ); + nSolvent = setup.Top().Nsolvent(); + } + mprintf("\t%i molecules out of %i selected as solvent.\n", nSolvent, setup.Top().Nmol()); // If there are no solvent molecules this action is not valid. - if (setup.Top().Nsolvent()==0) { - mprintf("Warning: Parm %s does not contain solvent.\n",setup.Top().c_str()); + if (nSolvent == 0) { + mprintf("Warning: Topology %s does not contain solvent.\n", setup.Top().c_str()); return Action::SKIP; } // If # solvent to keep >= solvent in this parm the action is not valid. - if (closestWaters_ == setup.Top().Nsolvent()) { + if (closestWaters_ == nSolvent) { mprintf("Warning: # solvent to keep (%i) == # solvent molecules in '%s' (%i)\n", - closestWaters_, setup.Top().c_str(), setup.Top().Nsolvent()); - } else if (closestWaters_ > setup.Top().Nsolvent()) { + closestWaters_, setup.Top().c_str(), nSolvent); + } else if (closestWaters_ > nSolvent) { mprintf("Warning: # solvent to keep (%i) > # solvent molecules in '%s' (%i)\n", - closestWaters_, setup.Top().c_str(), setup.Top().Nsolvent()); - closestWaters_ = setup.Top().Nsolvent(); + closestWaters_, setup.Top().c_str(), nSolvent); + closestWaters_ = nSolvent; mprintf("Warning: Keeping %i solvent molecules.\n", closestWaters_); } image_.SetupImaging( setup.CoordInfo().TrajBox().Type() ); @@ -171,8 +211,9 @@ Action::RetType Action_Closest::Setup(ActionSetup& setup) { MolDist solvent; solvent.D = 0.0; solvent.mol = 0; - SolventMols_.resize(setup.Top().Nsolvent(), solvent); + SolventMols_.resize(nSolvent, solvent); std::vector::iterator mdist = SolventMols_.begin(); + std::vector::iterator isSolvent = IsSolventMol.begin(); // 3: Set up the soluteMask for all non-solvent molecules. stripMask_.ResetMask(); int molnum = 1; @@ -181,12 +222,14 @@ Action::RetType Action_Closest::Setup(ActionSetup& setup) { int NsolventAtoms = -1; keptWaterAtomNum_.resize(closestWaters_); for (Topology::mol_iterator Mol = setup.Top().MolStart(); - Mol != setup.Top().MolEnd(); ++Mol) + Mol != setup.Top().MolEnd(); ++Mol, ++isSolvent) { - if ( !Mol->IsSolvent() ) { // Not solvent, add to solute mask. + if ( !(*isSolvent) ) { + // Not solvent, add to solute mask. stripMask_.AddAtomRange( Mol->BeginAtom(), Mol->EndAtom() ); newnatom += Mol->NumAtoms(); - } else { // Solvent, check for same # of atoms. + } else { + // Solvent, check for same # of atoms. if (NsolventAtoms == -1) NsolventAtoms = Mol->NumAtoms(); else if ( NsolventAtoms != Mol->NumAtoms() ) { @@ -197,7 +240,7 @@ Action::RetType Action_Closest::Setup(ActionSetup& setup) { } // mol here is the output molecule number which is why it starts from 1. mdist->mol = molnum; - // Solvent molecule mask + // Entire solvent molecule mask mdist->mask.AddAtomRange( Mol->BeginAtom(), Mol->EndAtom() ); // Atoms in the solvent molecule to actually calculate distances to. if (firstAtom_) { @@ -206,7 +249,8 @@ Action::RetType Action_Closest::Setup(ActionSetup& setup) { mdist->solventAtoms.clear(); mdist->solventAtoms.reserve( Mol->NumAtoms() ); for (int svatom = Mol->BeginAtom(); svatom < Mol->EndAtom(); svatom++) - mdist->solventAtoms.push_back( svatom ); + if (solventMask_.AtomInCharMask(svatom)) + mdist->solventAtoms.push_back( svatom ); } // For solvent molecules that will be kept, record what the atom number // will be in the new stripped parm. @@ -242,7 +286,7 @@ Action::RetType Action_Closest::Setup(ActionSetup& setup) { // Store original # of molecules. // NOTE: This is stored so that it can be used in the OpenMP section // of action. I dont think iterators are thread-safe. - NsolventMolecules_ = setup.Top().Nsolvent(); + NsolventMolecules_ = nSolvent; // Create stripped Parm if (newParm_!=0) delete newParm_; newParm_ = setup.Top().modifyStateByMask(stripMask_); diff --git a/src/Action_Closest.h b/src/Action_Closest.h index e34398317e..84d7dbf224 100644 --- a/src/Action_Closest.h +++ b/src/Action_Closest.h @@ -43,6 +43,7 @@ class Action_Closest: public Action { bool useMaskCenter_; ///< If true use geometric center of mask. AtomMask stripMask_; ///< Mask including all solute and closest molecules. AtomMask distanceMask_; ///< Mask of atoms to calculate distance from solvent to. + CharMask solventMask_; ///< Optional mask selecting solvent. Topology *newParm_; ///< New topology with solute and closest molecules. int NsolventMolecules_; ///< # of solvent molecules in SolventMols. int debug_; @@ -53,10 +54,10 @@ class Action_Closest: public Action { /** The moldist structure is used in order to preserve the original * solvent molecule numbers after sorting. */ struct MolDist { - int mol; ///< Original solvent molecule number (starts from 1). - double D; ///< Closest distance of solvent molecule to atoms in distanceMask. - AtomMask mask; ///< Original topology solvent molecule atom mask. - Iarray solventAtoms; ///< Actual solvent atom #s to loop over. + int mol; ///< Original solvent molecule number (starts from 1). + double D; ///< Closest distance of solvent molecule to atoms in distanceMask. + AtomMask mask; ///< Entire solvent molecule atom mask. + Iarray solventAtoms; ///< Actual solvent atom #s to calc distance to. }; /// Return true if the first molecule is closer than the second struct moldist_cmp { diff --git a/src/CharMask.cpp b/src/CharMask.cpp index b3696602f0..dae907376c 100644 --- a/src/CharMask.cpp +++ b/src/CharMask.cpp @@ -1,6 +1,14 @@ #include "CharMask.h" #include "CpptrajStdio.h" // PrintMaskAtoms +/** Use to initialize the mask without a mask expression. */ +void CharMask::InitCharMask(int natoms, bool initSelected) { + if (initSelected) + CharMask_.assign(natoms, SelectedChar_); + else + CharMask_.assign(natoms, UnselectedChar_); +} + /** Given atom and residue info and coordinates, setup character mask * based on current mask tokens. */ diff --git a/src/CharMask.h b/src/CharMask.h index 6454b35094..f9cf5f0b44 100644 --- a/src/CharMask.h +++ b/src/CharMask.h @@ -15,6 +15,8 @@ class CharMask : public MaskTokenArray { CharMask(std::string const& e) : nselected_(0) { SetMaskString(e); } /// Construct from given char array and # selected atoms. CharMask(std::vector const& c, int n) : CharMask_(c), nselected_(n) {} + /// Initialize mask with SelectedChar_ if true or UnselectedChar_ if false + void InitCharMask(int, bool); /// \return true if given atom is selected. bool AtomInCharMask(int) const; /// \return true if any atoms within given range are selected. diff --git a/src/Version.h b/src/Version.h index cd2e723a24..92da390c43 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.14.4" +#define CPPTRAJ_INTERNAL_VERSION "V4.14.5" /// PYTRAJ relies on this #define CPPTRAJ_VERSION_STRING CPPTRAJ_INTERNAL_VERSION #endif diff --git a/test/Test_Closest/RunTest.sh b/test/Test_Closest/RunTest.sh index 4281d6d4a9..e32ed23668 100755 --- a/test/Test_Closest/RunTest.sh +++ b/test/Test_Closest/RunTest.sh @@ -5,7 +5,8 @@ # Clean CleanFiles closest.in first.Closest.pdb.1 closestmols.dat \ closest.tz2.truncoct.parm7 all.Closest.pdb.1 \ - closest10.center2_4.crd closest10.mols.dat + closest10.center2_4.crd closest10.mols.dat \ + solventmask.dat wato.dat INPUT="-i closest.in" TESTNAME='Closest tests' @@ -50,6 +51,24 @@ RunCpptraj "Closest command test, using mask center" DoTest closest10.center2_4.crd.save closest10.center2_4.crd DoTest closest10.mols.dat.save closest10.mols.dat +# Test 4 - Select solvent via solventmask +cat > closest.in < closest.in <