diff --git a/regtest/adjmat/rt-dfg-wcsurf/surface.xyz.reference b/regtest/adjmat/rt-dfg-wcsurf/surface.xyz.reference
index b2fbdebe61386efbe22626c926c3d7004162a289..7e094a65483ff9fb018c35cdc145df1bf2729e23 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 46d7591ed91903ddcf6bbedb511d6753c65ee632..2cf006387e91c0746a16dc57dc6625bf70164cd7 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 001ec6ce9700a74330d934167d8aeb176ab1b8f7..ba1238444135c435ddf8d8ba5ad49631e7b6e502 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 9039b94d057caa6bbe7bd25645608a2a378ce3d9..2b3fb435414f849c5b94bad3a48f279929f2bf08 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 cf303df71bce0931ab25bae76760dcf479d4e6f3..f18680b180d403438e0c879baf7657b63e6a00c5 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 697bc822b0db55a31c0f2ba1ebd175a6d050075f..9d731faf63a1d686a3aaa582b7e70ae720fce701 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 );