From ba66676316d0c17b99736f69850d92f0b43afebd Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Fri, 22 Jan 2016 21:23:33 +0000
Subject: [PATCH] Improvements to multicolvar density to ensure it is updated
 in wcsurface

---
 regtest/adjmat/rt-dfg-wcsurf/surface.xyz.reference       | 2 +-
 .../crystallization/rt-wcsurface/surface.xyz.reference   | 2 +-
 src/gridtools/FindContour.cpp                            | 2 ++
 src/gridtools/GridVessel.cpp                             | 8 +++++---
 src/gridtools/GridVessel.h                               | 9 +++++++++
 src/multicolvar/MultiColvarDensity.cpp                   | 7 ++-----
 6 files changed, 20 insertions(+), 10 deletions(-)

diff --git a/regtest/adjmat/rt-dfg-wcsurf/surface.xyz.reference b/regtest/adjmat/rt-dfg-wcsurf/surface.xyz.reference
index b2fbdebe6..7e094a654 100644
--- a/regtest/adjmat/rt-dfg-wcsurf/surface.xyz.reference
+++ b/regtest/adjmat/rt-dfg-wcsurf/surface.xyz.reference
@@ -1,5 +1,5 @@
 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
diff --git a/regtest/crystallization/rt-wcsurface/surface.xyz.reference b/regtest/crystallization/rt-wcsurface/surface.xyz.reference
index 46d7591ed..2cf006387 100644
--- a/regtest/crystallization/rt-wcsurface/surface.xyz.reference
+++ b/regtest/crystallization/rt-wcsurface/surface.xyz.reference
@@ -1,5 +1,5 @@
 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
diff --git a/src/gridtools/FindContour.cpp b/src/gridtools/FindContour.cpp
index 001ec6ce9..ba1238444 100644
--- a/src/gridtools/FindContour.cpp
+++ b/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");
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index 9039b94d0..2b3fb4354 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/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 );
 }
 
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index cf303df71..f18680b18 100644
--- a/src/gridtools/GridVessel.h
+++ b/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
diff --git a/src/multicolvar/MultiColvarDensity.cpp b/src/multicolvar/MultiColvarDensity.cpp
index 697bc822b..9d731faf6 100644
--- a/src/multicolvar/MultiColvarDensity.cpp
+++ b/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 );
-- 
GitLab