Added first go at faster wcsurface stuff

Now there is an option to update grids only in the region around where
the contour was found in FindContour

regtest/multicolvar/rt-wcsurface/plumed.dat

-UNITS NATURAL
-
-fcc: FCCUBIC SPECIES=1-5184 SWITCH={CUBIC D_0=1.2 D_MAX=1.5} ALPHA=27  
-
-DUMPMULTICOLVAR DATA=fcc ORIGIN=1 FILE=atoms.xyz PRECISION=4
-dens: MULTICOLVARDENS DATA=fcc ORIGIN=1 DIR=xyz NBINS=14,14,28 BANDWIDTH=1.0,1.0,1.0 
-
-PRINT_CUBE GRID=dens STRIDE=1 FILE=dens.cube
-FIND_CONTOUR GRID=dens STRIDE=1 CONTOUR=0.42 FILE=surface.xyz PRECISION=4

src/gridtools/ActionWithInputGrid.cpp

@@ -78,6 +78,12 @@ void ActionWithInputGrid::setAnalysisStride( const bool& use_all, const unsigned
 void ActionWithInputGrid::update(){
   if( unormalised ) norm = 1.0;
   else norm=1.0/mygrid->getNorm();
+
+  if( checkAllActive() ){
+     for(unsigned i=0;i<mygrid->getNumberOfPoints();++i){
+         if( mygrid->inactive(i) ) error("if FIND_CONTOUR is used with BUFFER option then other actions cannot be performed with grid");
+     }
+  }
   performOperationsWithGrid( true );
 }
 

src/gridtools/ActionWithInputGrid.h

@@ -47,6 +47,7 @@ public:
   void apply(){}
   void update();
   void runFinalJobs();
+  virtual bool checkAllActive() const { return true; }
   virtual void performOperationsWithGrid( const bool& from_update )=0;
   void setAnalysisStride( const bool& use_all, const unsigned& astride );
 };

src/gridtools/FindContour.cpp

@@ -33,6 +33,7 @@ class FindContour : public ActionWithInputGrid {
 private:
   OFile of;
   double lenunit;
+  unsigned gbuffer;
   std::string fmt_xyz;
   unsigned mycomp;
   double contour;
@@ -40,6 +41,7 @@ private:
 public:
   static void registerKeywords( Keywords& keys );
   explicit FindContour(const ActionOptions&ao);
+  bool checkAllActive() const { return gbuffer==0; }
   void performOperationsWithGrid( const bool& from_update );
   double getDifferenceFromContour( const std::vector<double>& x, std::vector<double>& der );
   unsigned getNumberOfDerivatives(){ return 0; }
@@ -61,6 +63,7 @@ void FindContour::registerKeywords( Keywords& keys ){
                                           "is thre dimensional specifying SEARCHDIR=1,2 means that only points found along grid lines parallel to the "
                                           "first and second axis of the grid will be looked for.  The code will not search for the dividing surface "
                                           "along the third axis of the grid");
+  keys.add("compulsory","BUFFER","0","number of buffer grid points around location where grid was found on last step.  If this is zero the full grid is calculated on each step");
   keys.add("optional","COMPONENT","if your input is a vector field use this to specifiy the component of the input vector field for which you wish to find the contour");
   keys.add("optional", "PRECISION","The number of digits in trajectory file");  
 }
@@ -79,8 +82,9 @@ nosearch_dirs( mygrid->getDimension() )
   }
   if( mygrid->noDerivatives() ) error("cannot find contours if input grid has no derivatives");
 
-  parse("CONTOUR",contour);
+  parse("CONTOUR",contour); parse("BUFFER",gbuffer);
   log.printf("  calculating dividing surface along which function equals %f \n", contour);
+  if( gbuffer>0 ) log.printf("  after first step a subset of only %d grid points around where the countour was found will be checked\n",gbuffer);
 
   std::string searchdir_str; parse("SEARCHDIR",searchdir_str);
   if( searchdir_str=="all" ){
@@ -149,10 +153,28 @@ void FindContour::performOperationsWithGrid( const bool& from_update ){
   }
 
   // Run over whole grid
+  std::vector<bool> active( mygrid->getNumberOfPoints(), false );
+  std::vector<unsigned> neighbours; unsigned num_neighbours;
+  std::vector<unsigned> gbuffer_vec( mygrid->getDimension(), gbuffer );
+  std::vector<unsigned> ugrid_indices( mygrid->getDimension() );
   unsigned npoints=0; RootFindingBase<FindContour> mymin( this );
   for(unsigned i=0;i<mygrid->getNumberOfPoints();++i){
      // Get the index of the current grid point
      mygrid->getIndices( i, ind ); 
+     // Ensure inactive grid points are ignored
+     bool cycle=false;
+     if( mygrid->inactive(i) ) continue ;
+     for(unsigned j=0;j<mygrid->getDimension();++j) ugrid_indices[j]=ind[j];
+     for(unsigned j=0;j<mygrid->getDimension();++j){
+        if( mygrid->isPeriodic(j) ) ugrid_indices[j]=(ugrid_indices[j]+1)%nbin[j];
+        else if( (ugrid_indices[j]+1)==nbin[j] ) continue ;
+        else ugrid_indices[j]+=1;
+        // Now check the grid is active
+        if( mygrid->inactive( mygrid->getIndex( ugrid_indices ) ) ){ cycle=true; break; }
+        // And reset ugrid_indices
+        ugrid_indices[j]=ind[j];
+     }
+     if( cycle ) continue ;
 
      // Get the value of a point on the grid
      double val1=getGridElement( i, mycomp*(mygrid->getDimension()+1) ) - contour;
@@ -173,13 +195,20 @@ void FindContour::performOperationsWithGrid( const bool& from_update ){
              direction[j]=0.999999999*dx[j];
              // And do proper search for contour point
              mymin.linesearch( direction, contour_points[npoints], &FindContour::getDifferenceFromContour );
-             direction[j]=0.0; npoints++;
+             direction[j]=0.0; 
+             // For next run through activate buffer region around this grid point
+             if( gbuffer>0 ){
+                 mygrid->getNeighbors( contour_points[npoints], gbuffer_vec, num_neighbours, neighbours );
+                 for(unsigned n=0;n<num_neighbours;++n) active[ neighbours[n] ]=true; 
+             }
+             npoints++;
          }
          if( mygrid->isPeriodic(j) && edge ){ edge=false; ind[j]=nbin[j]-1; }
          else ind[j]-=1;
      } 
    
   }
+  if( gbuffer>0 ) mygrid->activateThesePoints( active );
   // Clear the grid ready for next time
   if( from_update ) mygrid->reset();
 

src/gridtools/GridVessel.cpp

@@ -120,7 +120,7 @@ std::string GridVessel::description(){
 
 void GridVessel::resize(){
   plumed_massert( nper>0, "Number of datapoints at each grid point has not been set");
-  resizeBuffer( npoints*nper ); data.resize( npoints*nper, 0 );
+  resizeBuffer( npoints*nper ); data.resize( npoints*nper, 0 ); active.resize( npoints, true );
 }
 
 unsigned GridVessel::getIndex( const std::vector<unsigned>& indices ) const {
@@ -189,7 +189,7 @@ void GridVessel::getSplineNeighbors( const unsigned& mybox, std::vector<unsigned
 }
 
 double GridVessel::getGridElement( const unsigned& ipoint, const unsigned& jelement ) const {
-  plumed_dbg_assert( bounds_set && ipoint<npoints && jelement<nper );
+  plumed_assert( bounds_set && ipoint<npoints && jelement<nper && active[ipoint] );
   return data[ nper*ipoint + jelement ];
 }
 
@@ -344,6 +344,11 @@ double GridVessel::getValueAndDerivatives( const std::vector<double>& x, const u
   return value;
 }
 
+void GridVessel::activateThesePoints( const std::vector<bool>& to_activate ){
+  plumed_dbg_assert( to_activate.size()==npoints );
+  for(unsigned i=0;i<npoints;++i) active[i]=to_activate[i];
+}
+
 }
 }
 

src/gridtools/GridVessel.h

@@ -50,8 +50,6 @@ private:
  std::vector<unsigned> stride;
 /// The number of bins in each grid direction
  std::vector<unsigned> nbin;
-///  Flatten the grid and get the grid index for a point
- unsigned getIndex( const std::vector<unsigned>& indices ) const ;
 /// The grid point that was requested last by getGridPointCoordinates
  unsigned currentGridPoint;
 protected:
@@ -68,11 +66,10 @@ protected:
  std::vector<double> dx;
 /// The dimensionality of the grid
  unsigned dimension;
+/// Which grid points are we actively accumulating
+ std::vector<bool> active;
 /// The grid with all the data values on it
  std::vector<double> data;
-/// Get the set of points neighouring a particular location in space
- void getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh, 
-                    unsigned& num_neighbours, std::vector<unsigned>& neighbors ) const ;
 /// Get the indices of a particular point
  void getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const ;
 /// Convert a point in space the the correspoinding grid point
@@ -90,6 +87,8 @@ public:
  std::string description();
 /// Convert an index into indices
  void convertIndexToIndices( const unsigned& index, const std::vector<unsigned>& nnbin, std::vector<unsigned>& indices ) const ;
+///  Flatten the grid and get the grid index for a point
+ unsigned getIndex( const std::vector<unsigned>& indices ) const ;
 /// Get the indices fof a point
  void getIndices( const unsigned& index, std::vector<unsigned>& indices ) const ;
 
@@ -127,6 +126,9 @@ public:
  double getCellVolume() const ;
 /// Get the value of the ith grid element 
  double getGridElement( const unsigned&, const unsigned& ) const ;
+/// Get the set of points neighouring a particular location in space
+ void getNeighbors( const std::vector<double>& pp, const std::vector<unsigned>& nneigh,
+                    unsigned& num_neighbours, std::vector<unsigned>& neighbors ) const ;
 /// Get the points neighboring a particular spline point
  void getSplineNeighbors( const unsigned& mybox, std::vector<unsigned>& mysneigh ) const ;
 /// Get the spacing between grid points
@@ -152,6 +154,10 @@ public:
  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 ;
+/// Deactivate all the grid points
+ void activateThesePoints( const std::vector<bool>& to_activate );
+/// Is this point active
+ bool inactive( const unsigned& ip ) const ;
 };
 
 inline
@@ -216,6 +222,12 @@ bool GridVessel::wasreset() const {
   return wascleared;
 }
 
+inline
+bool GridVessel::inactive( const unsigned& ip ) const {
+  plumed_dbg_assert( ip<npoints );
+  return !active[ip];
+}
+
 }
 }
 #endif

src/gridtools/HistogramOnGrid.cpp

@@ -85,6 +85,7 @@ bool HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, st
       double newval; std::vector<double> der( dimension );
       for(unsigned i=0;i<num_neigh;++i){
           unsigned ineigh=neighbors[i];
+          if( inactive( ineigh ) ) continue ;
           getGridPointCoordinates( ineigh, xx );
           for(unsigned j=0;j<dimension;++j) vv[j]->set(xx[j]);
           newval = kernel.evaluate( vv, der, true );