Amber-MD · drroe · May 22, 2019 · May 21, 2019 · May 21, 2019 · May 21, 2019
diff --git a/doc/CpptrajManual.pdf b/doc/CpptrajManual.pdf
diff --git a/doc/cpptraj.lyx b/doc/cpptraj.lyx
@@ -18812,19 +18812,15 @@ closest
 
 \end_inset
 
- <# to keep> <mask> [noimage] [first | oxygen] [center]
+ <# to keep> <mask> [solventmask <solvent mask>] [noimage]
 \end_layout
 
 \begin_layout LyX-Code
-        [closestout <filename>] [name <setname>] [outprefix <parmprefix>]
-\begin_inset Separator latexpar
-\end_inset
-
-
+    	[first | oxygen] [center] [closestout <filename> [name <setname>]]
 \end_layout
 
 \begin_layout LyX-Code
-        [parmout <file>]
+    	[outprefix <parmprefix>] [parmout <file>]
 \end_layout
 
 \begin_deeper
@@ -18844,6 +18840,27 @@ keep> Number of solvent molecules to keep around <mask>
 <mask> Mask of atoms to search for closest waters around.
 \end_layout
 
+\begin_layout Description
+[solventmask
+\begin_inset space ~
+\end_inset
+
+<solvent
+\begin_inset space ~
+\end_inset
+
+mask>] 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
 <mask>
 \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
@@ -39,9 +39,9 @@ Action_Closest::Action_Closest() :
 {} 
 
 void Action_Closest::Help() const {
-  mprintf("\t<# to keep> <mask> [noimage] [first | oxygen] [center]\n"
-          "\t[closestout <filename> [name <setname>]] [outprefix <parmprefix>]\n"
-          "\t[parmout <file>]\n"
+  mprintf("\t<# to keep> <mask> [solventmask <solvent mask>] [noimage]\n"
+          "\t[first | oxygen] [center] [closestout <filename> [name <setname>]]\n"
+          "\t[outprefix <parmprefix>] [parmout <file>]\n"
           "  Keep only the closest <# to keep> solvent molecules to atoms in <mask>.\n"
           "  Molecules can be marked as solvent with the 'solvent' command.\n"
           "  If 'center' specified use geometric center of atoms in <mask>.\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<bool> 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<MolDist>::iterator mdist = SolventMols_.begin();
+  std::vector<bool>::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
@@ -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
@@ -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
@@ -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<char> 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
@@ -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.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
@@ -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 <<EOF
+parm ../tz2.ortho.parm7
+trajin ../tz2.ortho.nc
+closest 10 :2,4 solventmask :WAT center closestout solventmask.dat name C10
+EOF
+RunCpptraj "Closest command test, using mask center and solvent mask"
+DoTest closest10.mols.dat.save solventmask.dat
+
+# Test 5 - Solventmask, oxygen atoms only
+  cat > closest.in <<EOF
+parm ../tz2.truncoct.parm7
+trajin ../tz2.truncoct.nc
+closest 10 :1-13 solventmask :WAT@O closestout wato.dat name CL
+EOF
+RunCpptraj "Closest command test, solvent mask, water oxygen only"
+DoTest wato.dat.save wato.dat
+
 EndTest
 
 exit 0