From e9492f9d16b6afb3c4c86931d5c96db23f12ffc5 Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Sat, 27 Feb 2016 16:51:54 +0000
Subject: [PATCH] Ensured that PLUMED crashes with an error if user tries to
 use too many atoms in multicolvar

---
 src/multicolvar/MultiColvar.cpp         | 16 ++-----------
 src/multicolvar/MultiColvarBase.cpp     | 30 ++++++++++++++-----------
 src/multicolvar/MultiColvarBase.h       |  2 ++
 src/multicolvar/MultiColvarFunction.cpp | 13 ++---------
 4 files changed, 23 insertions(+), 38 deletions(-)

diff --git a/src/multicolvar/MultiColvar.cpp b/src/multicolvar/MultiColvar.cpp
index 3a21f4520..b906e1a7a 100644
--- a/src/multicolvar/MultiColvar.cpp
+++ b/src/multicolvar/MultiColvar.cpp
@@ -102,20 +102,8 @@ void MultiColvar::readAtomsLikeKeyword( const std::string & key, int& natoms, st
      t.resize(0); 
   }
   if( all_atoms.size()>0 ){
-     nblock=ablocks[0].size(); 
-     if( natoms<4 ) resizeBookeepingArray( nblock, nblock ); 
-
-     for(unsigned i=0;i<nblock;++i){
-         if( natoms<4 ){
-            unsigned cvcode=0, tmpc=1; 
-            for(unsigned j=0;j<natoms;++j){ cvcode += i*tmpc; tmpc *= nblock; }
-            bookeeping(i,i).first=getFullNumberOfTasks();
-            addTaskToList( cvcode ); 
-            bookeeping(i,i).second=getFullNumberOfTasks();
-         } else {
-            addTaskToList( i );
-         }
-     }
+     nblock=0; 
+     for(unsigned i=0;i<ablocks[0].size();++i) addTaskToList( i );
   }
 }
 
diff --git a/src/multicolvar/MultiColvarBase.cpp b/src/multicolvar/MultiColvarBase.cpp
index ff8a2478e..bac4526d1 100644
--- a/src/multicolvar/MultiColvarBase.cpp
+++ b/src/multicolvar/MultiColvarBase.cpp
@@ -67,7 +67,8 @@ ActionWithValue(ao),
 ActionWithVessel(ao),
 usepbc(false),
 linkcells(comm),
-usespecies(false)
+usespecies(false),
+nblock(0)
 {
   if( keywords.exists("NOPBC") ){ 
     bool nopbc=!usepbc; parseFlag("NOPBC",nopbc);
@@ -90,13 +91,16 @@ void MultiColvarBase::resizeBookeepingArray( const unsigned& num1, const unsigne
 
 void MultiColvarBase::setupMultiColvarBase(){
   // Setup decoder array
-  if( !usespecies && ablocks.size()<4 ){
+  if( !usespecies && nblock>0 ){
+     // Check if number of atoms is too large
+     if( pow( nblock, ablocks.size()+1)>std::numeric_limits<unsigned>::max() ) error("number of atoms in groups is too big for PLUMED to handle"); 
+
      decoder.resize( ablocks.size() ); unsigned code=1;
      for(unsigned i=0;i<ablocks.size();++i){ decoder[ablocks.size()-1-i]=code; code *= nblock; } 
-     use_for_central_atom.resize( ablocks.size(), true );
+     ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
      numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
   } else if( !usespecies ){
-     use_for_central_atom.resize( ablocks.size(), true );
+     ncentral=ablocks.size(); use_for_central_atom.resize( ablocks.size(), true );
      numberForCentralAtom = 1.0 / static_cast<double>( ablocks.size() );
   }
 
@@ -110,7 +114,8 @@ void MultiColvarBase::setAtomsForCentralAtom( const std::vector<bool>& catom_ind
       use_for_central_atom[i]=catom_ind[i]; 
       if( use_for_central_atom[i] ) nat++;
   }
-  plumed_dbg_assert( nat>0 ); numberForCentralAtom = 1.0 / static_cast<double>( nat );
+  plumed_dbg_assert( nat>0 ); ncentral=nat;
+  numberForCentralAtom = 1.0 / static_cast<double>( nat );
 }
 
 void MultiColvarBase::turnOnDerivatives(){
@@ -125,7 +130,7 @@ void MultiColvarBase::setLinkCellCutoff( const double& lcut ){
 }
 
 void MultiColvarBase::setupLinkCells(){
-  if( !linkcells.enabled() ) return ;
+  if( (!usespecies && nblock==0) || !linkcells.enabled() ) return ;
 
   unsigned iblock;
   if( usespecies ){
@@ -234,7 +239,7 @@ void MultiColvarBase::setupLinkCells(){
 }
 
 void MultiColvarBase::decodeIndexToAtoms( const unsigned& taskCode, std::vector<unsigned>& atoms ) const {
-  plumed_dbg_assert( atoms.size()==ablocks.size() && !usespecies && ablocks.size()<4 );
+  plumed_dbg_assert( atoms.size()==ablocks.size() && !usespecies && nblock>0 );
   unsigned scode = taskCode;
   for(unsigned i=0;i<ablocks.size();++i){
       unsigned ind=( scode / decoder[i] );
@@ -248,7 +253,7 @@ bool MultiColvarBase::setupCurrentAtomList( const unsigned& taskCode, AtomValueP
      if( isDensity() ) return true;
      unsigned natomsper=myatoms.setupAtomsFromLinkCells( taskCode, getPositionOfAtomForLinkCells(taskCode), linkcells );
      return natomsper>1;
-  } else if( ablocks.size()<4 ){
+  } else if( nblock>0 ){
      std::vector<unsigned> atoms( ablocks.size() );
      decodeIndexToAtoms( taskCode, atoms ); myatoms.setNumberOfAtoms( ablocks.size() );
      for(unsigned i=0;i<ablocks.size();++i) myatoms.setAtom( i, atoms[i] ); 
@@ -292,7 +297,7 @@ Vector MultiColvarBase::getCentralAtomPos( const unsigned& taskIndex ){
 
   if( usespecies || isDensity() ){
      return getPositionOfAtomForLinkCells(curr);
-  } else if( ablocks.size()<4 ){
+  } else if( nblock>0 ){
      // double factor=1.0/static_cast<double>( ablocks.size() );
      Vector mypos; mypos.zero(); 
      std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
@@ -317,9 +322,8 @@ CatomPack MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsig
      mypack.resize(1);
      mypack.setIndex( 0, basn + curr );
      mypack.setDerivative( 0, Tensor::identity() );
-  } else if( ablocks.size()<4 ){
-     mypack.resize(ablocks.size());
-     unsigned k=0;
+  } else if( nblock>0 ){
+     mypack.resize(ncentral); unsigned k=0;
      std::vector<unsigned> atoms( ablocks.size() ); decodeIndexToAtoms( curr, atoms );
      for(unsigned i=0;i<ablocks.size();++i){
          if( use_for_central_atom[i] ){
@@ -329,7 +333,7 @@ CatomPack MultiColvarBase::getCentralAtomPack( const unsigned& basn, const unsig
          }
      }
   } else {
-     unsigned k=0;
+     mypack.resize(ncentral); unsigned k=0;
      for(unsigned i=0;i<ablocks.size();++i){
          if( use_for_central_atom[i] ){
              mypack.setIndex( k, basn + ablocks[i][curr] );
diff --git a/src/multicolvar/MultiColvarBase.h b/src/multicolvar/MultiColvarBase.h
index 227e9ef6e..c47c89b38 100644
--- a/src/multicolvar/MultiColvarBase.h
+++ b/src/multicolvar/MultiColvarBase.h
@@ -56,6 +56,8 @@ private:
   LinkCells linkcells;
 /// This remembers where the boundaries are for the tasks. It makes link cells work fast
   Matrix<std::pair<unsigned,unsigned> > bookeeping;
+/// Number of atoms that are being used for central atom position
+  unsigned ncentral;
 /// Bool vector telling us which atoms are required to calculate central atom position
   std::vector<bool> use_for_central_atom;
 /// 1/number of atoms involved in central atoms
diff --git a/src/multicolvar/MultiColvarFunction.cpp b/src/multicolvar/MultiColvarFunction.cpp
index cbcd2aa8c..da48b45d1 100644
--- a/src/multicolvar/MultiColvarFunction.cpp
+++ b/src/multicolvar/MultiColvarFunction.cpp
@@ -129,23 +129,14 @@ void MultiColvarFunction::buildSets(){
      if( mybasemulticolvars[i]->getFullNumberOfTasks()!=nblock ){
           error("mismatch between numbers of tasks in various base multicolvars");
      }
-//     mybasemulticolvars[i]->buildDataStashes( false, 0.0 );
   }
-  ablocks.resize( mybasemulticolvars.size() );
+  nblock=0; ablocks.resize( mybasemulticolvars.size() );
   usespecies=false; // current_atoms.resize( mybasemulticolvars.size() );
   for(unsigned i=0;i<mybasemulticolvars.size();++i){
       ablocks[i].resize( nblock ); 
       for(unsigned j=0;j<nblock;++j) ablocks[i][j]=i*nblock+j;  
   }
-  for(unsigned i=0;i<nblock;++i){
-      if( mybasemulticolvars.size()<4 ){
-          unsigned cvcode=0, tmpc=1;
-          for(unsigned j=0;j<ablocks.size();++j){ cvcode +=i*tmpc; tmpc *= nblock; }
-          addTaskToList( cvcode );
-      } else {
-          addTaskToList( i );
-      }
-  }
+  for(unsigned i=0;i<nblock;++i) addTaskToList( i );
   setupAtomLists();
 }
 
-- 
GitLab