Improvements to multicolvar density to ensure it is updated in wcsurface

regtest/adjmat/rt-dfg-wcsurf/surface.xyz.reference

 676
-   28.6600   28.6600   28.6600
+Points found on isocontour
 X   -7.6427   -0.0443  -14.3300
 X   -6.6873   -0.2252  -14.3300
 X   -5.7320   -0.2733  -14.3300

regtest/crystallization/rt-wcsurface/surface.xyz.reference

 582
-   14.5594   14.5594   28.1255
+Points found on isocontour
 X    5.1998   -3.1199  -10.1655
 X    6.2398   -3.1199  -10.0459
 X    4.1598   -2.0799  -10.1258

src/gridtools/FindContour.cpp

@@ -153,6 +153,8 @@ void FindContour::performOperationsWithGrid( const bool& from_update ){
      } 
    
   }
+  // Clear the grid ready for next time
+  if( from_update ) mygrid->clear();
 
   of.printf("%u\n",npoints);
   of.printf("Points found on isocontour\n");

src/gridtools/GridVessel.cpp

@@ -37,6 +37,7 @@ void GridVessel::registerKeywords( Keywords& keys ){
 GridVessel::GridVessel( const vesselbase::VesselOptions& da ):
 Vessel(da),
 noderiv(false),
+wascleared(true), 
 bounds_set(false),
 cube_units(1.0)
 {
@@ -76,7 +77,7 @@ void GridVessel::setNoDerivatives(){
 void GridVessel::setBounds( const std::vector<std::string>& smin, const std::vector<std::string>& smax,
                             const std::vector<unsigned>& binsin, const std::vector<double>& spacing ){
   plumed_dbg_assert( smin.size()==dimension && smax.size()==dimension );
-  plumed_assert( spacing.size()==dimension || binsin.size()==dimension );
+  plumed_assert( wascleared && spacing.size()==dimension || binsin.size()==dimension );
 
   npoints=1; bounds_set=true;
   stride.resize( dimension ); max.resize( dimension );
@@ -194,12 +195,12 @@ double GridVessel::getGridElement( const unsigned& ipoint, const unsigned& jelem
 
 void GridVessel::setGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ){
   plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
-  data[ nper*ipoint + jelement ] = value;
+  wascleared=false; data[ nper*ipoint + jelement ] = value;
 }
 
 void GridVessel::addToGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value ){
   plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
-  data[ nper*ipoint + jelement ] += value;
+  wascleared=false; data[ nper*ipoint + jelement ] += value;
 }
 
 double GridVessel::getGridElement( const std::vector<unsigned>& indices, const unsigned& jelement ) const {
@@ -267,6 +268,7 @@ void GridVessel::getNeighbors( const std::vector<double>& pp, const std::vector<
 
 void GridVessel::clear(){
   if( !nomemory ) return ;
+  wascleared=true;
   data.assign( data.size(), 0.0 );
 }
 

src/gridtools/GridVessel.h

@@ -34,6 +34,8 @@ class GridVessel : public vesselbase::Vessel {
 private:
 /// Do we have derivatives
  bool noderiv;
+/// The grid was recently cleared and bounds can be set
+ bool wascleared;
 /// Have the minimum and maximum for the grid been set
  bool bounds_set;
 /// The number of points in the grid
@@ -146,6 +148,8 @@ public:
  bool noDerivatives() const ;
 /// Get the value and derivatives at a particular location using spline interpolation
  double getValueAndDerivatives( const std::vector<double>& x, const unsigned& ind, std::vector<double>& der ) const ; 
+/// Was the grid cleared on the last step
+ bool wasreset() const ;
 };
 
 inline
@@ -205,6 +209,11 @@ bool GridVessel::noDerivatives() const {
   return noderiv;
 }
 
+inline
+bool GridVessel::wasreset() const {
+  return wascleared;
+}
+
 }
 }
 #endif

src/multicolvar/MultiColvarDensity.cpp

@@ -59,7 +59,6 @@ class MultiColvarDensity :
   public vesselbase::ActionWithInputVessel
 {
   std::string kerneltype;
-  bool firststep;
   bool fractional;
   unsigned rstride;
   MultiColvarBase* mycolv; 
@@ -107,8 +106,7 @@ MultiColvarDensity::MultiColvarDensity(const ActionOptions&ao):
   ActionPilot(ao),
   ActionAtomistic(ao),
   ActionWithVessel(ao),
-  ActionWithInputVessel(ao),
-  firststep(true)
+  ActionWithInputVessel(ao)
 {
 
   std::vector<AtomNumber> atom;
@@ -191,7 +189,7 @@ unsigned MultiColvarDensity::getNumberOfQuantities() const {
 }
 
 void MultiColvarDensity::update(){
-  if(firststep){
+  if( mygrid->wasreset() ){
      std::vector<double> min(directions.size()), max(directions.size());
      std::vector<std::string> gmin(directions.size()), gmax(directions.size());;
      for(unsigned i=0;i<directions.size();++i){ min[i]=-0.5; max[i]=0.5; }
@@ -210,7 +208,6 @@ void MultiColvarDensity::update(){
      } else {
         mygrid->setBounds( gmin, gmax, nbins, gspacing ); resizeFunctions();
      }
-     firststep=false;    // We only have the first step once
   } else {
       for(unsigned i=0;i<directions.size();++i){
           double max; Tools::convert( mygrid->getMax()[i], max );