Added facility to specify a subset of central atoms for multicolvar functions

src/multicolvar/MultiColvarBase.cpp

@@ -197,6 +197,18 @@ Vector MultiColvarBase::getSeparation( const Vector& vec1, const Vector& vec2 )
   else{ return delta( vec1, vec2 ); }
 }
 
+unsigned MultiColvarBase::getInternalIndex( const AtomNumber& iatom ) const {
+  unsigned katom; bool found=false;
+  for(unsigned i=0;i<colvar_atoms.size();++i){
+      unsigned jatom = colvar_atoms[i][0];
+      if( all_atoms[ all_atoms.linkIndex(jatom) ]==iatom ){
+         katom=i; found=true;
+      }
+  }
+  if(!true) error("could not find required atom in any of the quantities calculated by the base multicolvar");
+  return katom;
+}
+
 void MultiColvarBase::getIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& indices ){
   indices[jstore]=3*atoms_with_derivatives.getNumberActive() + 9;
   if( indices[jstore]>maxder ) error("too many derivatives to store. Run with LOWMEM");

src/multicolvar/MultiColvarBase.h

@@ -136,6 +136,8 @@ public:
   StoreCentralAtomsVessel* getCentralAtoms();
 /// Copy the list of atoms involved to a second MultiColvarBase (used by functions)
   void copyAtomListToFunction( MultiColvarBase* myfunction );
+/// Return the number of the colvar in which iatom is the first atom
+  unsigned getInternalIndex( const AtomNumber& iatom ) const ;
 /// Make sure the same list of atoms is active in a function
   void copyActiveAtomsToFunction( MultiColvarBase* myfunction );
 /// Activate the atoms that have derivatives from a storeDataVessel

src/multicolvar/MultiColvarFunction.cpp

@@ -30,6 +30,7 @@ namespace multicolvar {
 void MultiColvarFunction::registerKeywords( Keywords& keys ){
   MultiColvarBase::registerKeywords( keys );
   keys.add("compulsory","ARG","the label of the action that calculates the vectors we are interested in");
+  keys.reserve("atoms","ONLY_ATOMS","only calculate this quantity for these central atoms");
 }
 
 MultiColvarFunction::MultiColvarFunction(const ActionOptions& ao):
@@ -51,6 +52,20 @@ MultiColvarBase(ao)
 
   // Retrieve the central atoms
   catoms = mycolv->getCentralAtoms();
+
+  if( keywords.exists("ONLY_ATOMS") ){
+      std::vector<AtomNumber> atoms; parseAtomList("ONLY_ATOMS",atoms);
+      if( atoms.size()>0 ){ 
+         log.printf("  calculating function for atoms ");
+         for(unsigned i=0;i<atoms.size();++i){
+            log.printf("%d ",atoms[i].serial() );
+            subatomlist.push_back( mycolv->getInternalIndex(atoms[i]) );
+         }
+         log.printf("\n");
+      } else {
+         for(unsigned i=0;i<mycolv->colvar_atoms.size();++i) subatomlist.push_back( i );
+      }
+  } 
 }
 
 void MultiColvarFunction::completeSetup(){

src/multicolvar/MultiColvarFunction.h

@@ -30,11 +30,17 @@ namespace multicolvar {
 
 class MultiColvarFunction : public MultiColvarBase {
 private:
+/// The list of atoms that are involved in this colvar
+  std::vector<unsigned> subatomlist;
 /// The multicolvar from which we construct these quantities
   multicolvar::MultiColvarBase* mycolv;
 /// The central atom positions
   multicolvar::StoreCentralAtomsVessel* catoms;
 protected:
+/// Get the number of colvar functions we are calculating
+  unsigned getNumberOfColvarFunctions() const;
+/// Get the index of the ith colvar we are using
+  unsigned getColvarIndex( const unsigned& ) const;
 /// Get the number of functions in the multicolvar we are operating on
   unsigned getNumberOfBaseFunctions() const;
 /// Return a pointer to the multicolvar we are using as a base function
@@ -80,6 +86,18 @@ unsigned MultiColvarFunction::getNumberOfBaseFunctions() const {
   return mycolv->colvar_atoms.size();
 }
 
+inline
+unsigned MultiColvarFunction::getNumberOfColvarFunctions() const {
+  plumed_dbg_assert( subatomlist.size()>0 );
+  return subatomlist.size();
+}
+
+inline
+unsigned MultiColvarFunction::getColvarIndex( const unsigned& jindex ) const {
+  plumed_dbg_assert( subatomlist.size()>0 && jindex<subatomlist.size() );
+  return subatomlist[jindex];
+}
+
 inline
 unsigned MultiColvarFunction::getAtomIndex( const unsigned& iatom ) const {
   plumed_dbg_assert( iatom<colvar_atoms[current].getNumberActive() );