Handle non-contiguous molecules

src/Action_AreaPerMol.cpp

@@ -77,7 +77,7 @@ Action::RetType Action_AreaPerMol::Setup(ActionSetup& setup) {
     for (Topology::mol_iterator mol = setup.Top().MolStart();
                                 mol != setup.Top().MolEnd(); ++mol)
     {
-      if (Mask1_.AtomsInCharMask(mol->BeginAtom(), mol->EndAtom()))
+      if (Mask1_.AtomsInCharMask(mol->MolUnit()))
         Nmols_ += 1.0;
     }
     mprintf("\tMask '%s' selects %.0f molecules.\n", Mask1_.MaskString(), Nmols_);

src/Action_AutoImage.cpp

@@ -4,6 +4,7 @@
 #include "DistRoutines.h"
 #include "ImageRoutines.h"
 #include "CharMask.h"
+#include "Image_List_Unit.h"
 
 // CONSTRUCTOR
 Action_AutoImage::Action_AutoImage() :
@@ -14,9 +15,17 @@ Action_AutoImage::Action_AutoImage() :
   truncoct_(false),
   useMass_(false),
   movingAnchor_(false),
-  triclinic_(OFF)
+  triclinic_(OFF),
+  fixedList_(0),
+  mobileList_(0)
 {}
 
+/** DESTRUCTOR */
+Action_AutoImage::~Action_AutoImage() {
+  if (fixedList_ != 0) delete fixedList_;
+  if (mobileList_ != 0) delete mobileList_;
+}
+
 void Action_AutoImage::Help() const {
   mprintf("\t[<mask> | anchor <mask> [fixed <fmask>] [mobile <mmask>]]\n"
           "\t[origin] [firstatom] [familiar | triclinic] [moveanchor]\n"
@@ -74,42 +83,6 @@ Action::RetType Action_AutoImage::Init(ArgList& actionArgs, ActionInit& init, in
   return Action::OK;
 }
 
-// Action_AutoImage::SetupAtomRanges()
-/** Based on the given atom mask expression determine what molecules are
-  * selected by the mask. If a mask selects any part of a molecule the
-  * entire molecule will be selected.
-  * \return A list of atom pairs that mark the beginning and end of each
-  *         selected molecule.
-  */
-Action_AutoImage::pairList
-  Action_AutoImage::SetupAtomRanges(Topology const& currentParm, std::string const& maskexpr)
-{
-  pairList imageList;
-  CharMask Mask1( maskexpr.c_str() );
-
-  if (currentParm.SetupCharMask( Mask1 )) return imageList;
-  if (Mask1.None()) return imageList;
-  for (Topology::mol_iterator mol = currentParm.MolStart(); mol != currentParm.MolEnd(); mol++)
-  {
-    int firstAtom = mol->BeginAtom();
-    int lastAtom = mol->EndAtom();
-    bool rangeIsValid = false;
-    // Check that any atom in the range is in Mask1
-    for (int atom = firstAtom; atom < lastAtom; ++atom) {
-      if (Mask1.AtomInCharMask(atom)) {
-        rangeIsValid = true;
-        break;
-      }
-    }
-    if (rangeIsValid) {
-      imageList.push_back( firstAtom );
-      imageList.push_back( lastAtom );
-    }
-  }
-  mprintf("\tMask [%s] corresponds to %zu molecules\n", Mask1.MaskString(), imageList.size()/2);
-  return imageList;
-}
-
 // Action_AutoImage::Setup()
 Action::RetType Action_AutoImage::Setup(ActionSetup& setup) {
   bool fixedauto = false;
@@ -159,23 +132,37 @@ Action::RetType Action_AutoImage::Setup(ActionSetup& setup) {
     // No anchor specified. Use first molecule as anchor.
     anchormolnum = 0;
     mprintf("\tUsing first molecule as anchor.\n");
-    anchorMask_.AddAtomRange( setup.Top().Mol(0).BeginAtom(),
-                              setup.Top().Mol(0).EndAtom()    );
+    anchorMask_.AddUnit( setup.Top().Mol(0).MolUnit() );
   }
 
+  if (fixedList_ != 0) delete fixedList_;
+  if (mobileList_ != 0) delete mobileList_;
   // Set up fixed region
-  if (!fixed_.empty()) 
-    fixedList_ = SetupAtomRanges( setup.Top(), fixed_ );
+  // NOTE: Always use molecule center when imaging fixed list
+  if (!fixed_.empty())
+    fixedList_ = (Image::List_Unit*)
+                 Image::CreateImageList(setup.Top(), Image::BYMOL, fixed_, useMass_, true); 
   else { 
     fixedauto = true;
-    fixedList_.clear();
+    fixedList_ = (Image::List_Unit*)
+                 Image::CreateImageList(Image::BYMOL, useMass_, true);
+  }
+  if (fixedList_ == 0) {
+    mprinterr("Internal Error: Could not allocate fixed list.\n");
+    return Action::ERR;
   }
   // Set up mobile region
   if (!mobile_.empty())
-    mobileList_ = SetupAtomRanges( setup.Top(), mobile_ );
+    mobileList_ = (Image::List_Unit*)
+                  Image::CreateImageList(setup.Top(), Image::BYMOL, mobile_, useMass_, usecom_);
   else {
     mobileauto = true;
-    mobileList_.clear();
+    mobileList_ = (Image::List_Unit*)
+                  Image::CreateImageList(Image::BYMOL, useMass_, usecom_);
+  }
+  if (mobileList_ == 0) {
+    mprinterr("Internal Error: Could not allocate mobile list.\n");
+    return Action::ERR;
   }
   // Automatic search through molecules for fixed/mobile
   if (fixedauto || mobileauto) {
@@ -189,33 +176,31 @@ Action::RetType Action_AutoImage::Setup(ActionSetup& setup) {
         // everything else into fixed list.
         if ( mol->IsSolvent() || mol->NumAtoms() == 1 ) {
           if (mobileauto) {
-            mobileList_.push_back( mol->BeginAtom() );
-            mobileList_.push_back( mol->EndAtom()   );
+            mobileList_->AddUnit( mol->MolUnit() );
           }
         } else {
           if (fixedauto) {
-            fixedList_.push_back( mol->BeginAtom() );
-            fixedList_.push_back( mol->EndAtom()   );
+            fixedList_->AddUnit( mol->MolUnit() );
           }
         }
       }
       ++molnum;
     }
   }
   // Print fixed and mobile lists
-  if (!fixedList_.empty()) {
-    mprintf("\t%zu molecules are fixed to anchor:", fixedList_.size() / 2);
-    for (pairList::const_iterator atom = fixedList_.begin();
-                                  atom != fixedList_.end(); atom += 2)
-      mprintf(" %i", setup.Top()[ *atom ].MolNum()+1 );
+  if (!fixedList_->empty()) {
+    mprintf("\t%u molecules are fixed to anchor:", fixedList_->nEntities());
+    for (Image::List_Unit::const_iterator it = fixedList_->begin();
+                                it != fixedList_->end(); ++it)
+      mprintf(" %i", setup.Top()[ it->Front() ].MolNum()+1 );
     mprintf("\n");
   }
-  mprintf("\t%zu molecules are mobile.\n", mobileList_.size() / 2 );
+  mprintf("\t%u molecules are mobile.\n", mobileList_->nEntities() );
   if (debug_ > 1) {
     mprintf("\tThe following molecules are mobile:\n");
-    for (pairList::const_iterator atom = mobileList_.begin();
-                                  atom != mobileList_.end(); atom += 2)
-      mprintf(" %i\n", setup.Top()[ *atom ].MolNum()+1 );
+    for (Image::List_Unit::const_iterator it = mobileList_->begin();
+                                it != mobileList_->end(); ++it)
+      mprintf(" %i\n", setup.Top()[ it->Front() ].MolNum()+1 );
     mprintf("\n");
   }
 
@@ -264,28 +249,22 @@ Action::RetType Action_AutoImage::DoAction(int frameNum, ActionFrame& frm) {
       // TODO: Return OK for now so next frame is tried; eventually indicate SKIP?
       return Action::OK; // FIXME return MODIFY_COORDS instead?
     }
-    Image::Ortho(frm.ModifyFrm(), bp, bm, offset, usecom_, useMass_, mobileList_);
+    Image::Ortho(frm.ModifyFrm(), bp, bm, offset, *mobileList_);
   } else {
     if (truncoct_)
       fcom = Image::SetupTruncoct( frm.Frm(), 0, useMass_, origin_ );
     Image::Nonortho(frm.ModifyFrm(), origin_, fcom, offset, ucell, recip, truncoct_,
-                    usecom_, useMass_, mobileList_);
+                    *mobileList_);
   }
 
   if (movingAnchor_) {
     // TODO I think the way the translation is calculated here is robust and
     //      more efficient than the !movingAnchor_ case but more testing is
     //      needed.
     // Loop over fixed molecules
-    for (pairList::const_iterator atom1 = fixedList_.begin();
-                                  atom1 != fixedList_.end(); atom1 += 2)
+    for (unsigned int idx = 0; idx != fixedList_->nEntities(); ++idx)
     {
-      int firstAtom = *atom1;
-      int lastAtom = *(atom1+1);
-      if (useMass_)
-        framecenter = frm.Frm().VCenterOfMass(firstAtom, lastAtom);
-      else
-        framecenter = frm.Frm().VGeometricCenter(firstAtom, lastAtom);
+      framecenter = fixedList_->GetCoord(idx, frm.Frm());
 
       // Determine distance in terms of box lengths
       if (ortho_) {
@@ -300,7 +279,7 @@ Action::RetType Action_AutoImage::DoAction(int frameNum, ActionFrame& frm) {
         //        frameNum, (atom1-fixedList_.begin())/2, firstAtom+1, lastAtom,
         //        minTrans[0], minTrans[1], minTrans[2]);
         // Move atoms closer to anchor. Update coords in currentFrame.
-        frm.ModifyFrm().Translate(minTrans, firstAtom, lastAtom);
+        fixedList_->DoTranslation(frm.ModifyFrm(), idx, minTrans);
         // New anchor is previous fixed mol
         anchorcenter = minImage;
       } else {
@@ -317,7 +296,7 @@ Action::RetType Action_AutoImage::DoAction(int frameNum, ActionFrame& frm) {
           //        Trans[0], Trans[1], Trans[2], sqrt(framedist2), sqrt(imageddist2));
           if (imageddist2 < framedist2) {
             // Imaging these atoms moved them closer to anchor. Update coords in currentFrame.
-            frm.ModifyFrm().Translate(Trans, firstAtom, lastAtom);
+            fixedList_->DoTranslation(frm.ModifyFrm(), idx, Trans);
             newAnchor = imagedcenter;
           }
         }
@@ -328,15 +307,10 @@ Action::RetType Action_AutoImage::DoAction(int frameNum, ActionFrame& frm) {
     // For each molecule defined by atom pairs in fixedList, determine if the
     // imaged position is closer to anchor center than the current position.
     // Always use molecule center when imaging fixedList.
-    for (pairList::const_iterator atom1 = fixedList_.begin();
-                                  atom1 != fixedList_.end(); atom1 += 2)
+    for (unsigned int idx = 0; idx != fixedList_->nEntities(); ++idx)
     {
-      int firstAtom = *atom1;
-      int lastAtom = *(atom1+1);
-      if (useMass_)
-        framecenter = frm.Frm().VCenterOfMass(firstAtom, lastAtom);
-      else
-        framecenter = frm.Frm().VGeometricCenter(firstAtom, lastAtom);
+      framecenter = fixedList_->GetCoord(idx, frm.Frm());
+
       // Determine if molecule would be imaged.
       if (ortho_)
         Trans = Image::Ortho(framecenter, bp, bm, box);
@@ -352,7 +326,7 @@ Action::RetType Action_AutoImage::DoAction(int frameNum, ActionFrame& frm) {
 //                Trans[0], Trans[1], Trans[2], sqrt(framedist2), sqrt(imageddist2));
         if (imageddist2 < framedist2) {
           // Imaging these atoms moved them closer to anchor. Update coords in currentFrame.
-          frm.ModifyFrm().Translate(Trans, firstAtom, lastAtom);
+          fixedList_->DoTranslation(frm.ModifyFrm(), idx, Trans);
         }
       }
     }

src/Action_AutoImage.h

@@ -1,9 +1,14 @@
 #ifndef INC_ACTION_AUTOIMAGE_H
 #define INC_ACTION_AUTOIMAGE_H
 #include "Action.h"
+namespace Image {
+  class List_Unit;
+}
+/// Perform imaging, attempting to keep solute in 1 configuration
 class Action_AutoImage : public Action {
   public:
     Action_AutoImage();
+    ~Action_AutoImage();
     DispatchObject* Alloc() const { return (DispatchObject*)new Action_AutoImage(); }
     void Help() const;
   private:
@@ -12,10 +17,6 @@ class Action_AutoImage : public Action {
     Action::RetType DoAction(int, ActionFrame&);
     void Print() {}
 
-    typedef std::vector<int> pairList;
-
-    static pairList SetupAtomRanges(Topology const&, std::string const&);
-
     AtomMask anchorMask_; ///< Used to center anchor region.
     std::string anchor_;  ///< Mask expression for anchor region.
     std::string fixed_;   ///< Mask expression for fixed region.
@@ -30,7 +31,7 @@ class Action_AutoImage : public Action {
     enum TriclinicArg {OFF, FORCE, FAMILIAR};
     TriclinicArg triclinic_; ///< Determine whether triclinic code should be used.
 
-    pairList fixedList_;  ///< Contain first and last atom indices for fixed elements
-    pairList mobileList_; ///< Contain first and last atom indices for mobile elements.
+    Image::List_Unit* fixedList_;  ///< Contain atom indices for fixed elements
+    Image::List_Unit* mobileList_; ///< Contain atom indices for mobile elements.
 };
 #endif

src/Action_Closest.cpp

@@ -160,7 +160,7 @@ Action::RetType Action_Closest::Setup(ActionSetup& setup) {
     for (Topology::mol_iterator Mol = setup.Top().MolStart();
                                 Mol != setup.Top().MolEnd(); ++Mol)
     {
-      if ( solventMask_.AtomsInCharMask(Mol->BeginAtom(), Mol->EndAtom()) ) {
+      if ( solventMask_.AtomsInCharMask(Mol->MolUnit()) ) {
         IsSolventMol.push_back( true );
         nSolvent++;
       } else
@@ -214,44 +214,46 @@ Action::RetType Action_Closest::Setup(ActionSetup& setup) {
   int molnum = 1;
   int newnatom = 0;
   int nclosest = 0;
-  int NsolventAtoms = -1;
+  unsigned int NsolventAtoms = 0;
   keptWaterAtomNum_.resize(closestWaters_);
   for (Topology::mol_iterator Mol = setup.Top().MolStart();
                               Mol != setup.Top().MolEnd(); ++Mol, ++isSolvent)
   {
     if ( !(*isSolvent) ) {
       // Not solvent, add to solute mask.
-      stripMask_.AddAtomRange( Mol->BeginAtom(), Mol->EndAtom() );
+      stripMask_.AddUnit( Mol->MolUnit() );
       newnatom += Mol->NumAtoms();
     } else {
       // Solvent, check for same # of atoms.
-      if (NsolventAtoms == -1)
+      if (NsolventAtoms == 0)
         NsolventAtoms = Mol->NumAtoms();
       else if ( NsolventAtoms != Mol->NumAtoms() ) {
         mprinterr("Error: Solvent molecules in '%s' are not of uniform size.\n"
-                  "Error:   First solvent mol = %i atoms, solvent mol %i = %i atoms.\n",
-                  setup.Top().c_str(), NsolventAtoms, molnum, (*Mol).NumAtoms());
+                  "Error:   First solvent mol = %u atoms, solvent mol %i = %u atoms.\n",
+                  setup.Top().c_str(), NsolventAtoms, molnum, Mol->NumAtoms());
         return Action::ERR;
       }
       // mol here is the output molecule number which is why it starts from 1.
       mdist->mol = molnum;
       // Entire solvent molecule mask
-      mdist->mask.AddAtomRange( Mol->BeginAtom(), Mol->EndAtom() );
+      mdist->mask.AddUnit( Mol->MolUnit() );
       // Atoms in the solvent molecule to actually calculate distances to.
       if (firstAtom_) {
-        mdist->solventAtoms.assign(1, Mol->BeginAtom() );
+        mdist->solventAtoms.assign(1, Mol->MolUnit().Front() );
       } else {
         mdist->solventAtoms.clear();
         mdist->solventAtoms.reserve( Mol->NumAtoms() );
-        for (int svatom = Mol->BeginAtom(); svatom < Mol->EndAtom(); svatom++)
-          if (solventMask_.AtomInCharMask(svatom))
-            mdist->solventAtoms.push_back( svatom );
+        for (Unit::const_iterator seg = Mol->MolUnit().segBegin();
+                                  seg != Mol->MolUnit().segEnd(); ++seg)
+          for (int svatom = seg->Begin(); svatom < seg->End(); 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.
       if (nclosest < closestWaters_) {
         keptWaterAtomNum_[nclosest] = newnatom;
-        stripMask_.AddAtomRange( Mol->BeginAtom(), Mol->EndAtom() );
+        stripMask_.AddUnit( Mol->MolUnit() );
         newnatom += Mol->NumAtoms();
         ++nclosest;
       }