From 264af882870f9338fd4bd9fe1bab4edcbf0669cd Mon Sep 17 00:00:00 2001
From: Gareth Tribello <gareth.tribello@gmail.com>
Date: Thu, 20 Apr 2017 22:00:08 -0300
Subject: [PATCH] Fixed small bug in discrete histogram that causes nans
---
src/analysis/Histogram.cpp | 2 +-
src/gridtools/GridVessel.cpp | 5 +++++
src/gridtools/GridVessel.h | 4 ++--
src/gridtools/HistogramOnGrid.cpp | 2 +-
4 files changed, 9 insertions(+), 4 deletions(-)
diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp
index a2c9c7e23..cee0a87e3 100644
--- a/src/analysis/Histogram.cpp
+++ b/src/analysis/Histogram.cpp
@@ -298,7 +298,7 @@ void Histogram::prepareForAveraging(){
lockContributors();
} else {
// This is used when we are not doing kernel density evaluation
- mygrid->setGridElement( neighbors[0], 0, mygrid->getGridElement( neighbors[0], 0 ) + cweight );
+ mygrid->addToGridElement( neighbors[0], 0, cweight );
}
}
}
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index 55def546f..82a6a6ce6 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/src/gridtools/GridVessel.cpp
@@ -201,6 +201,11 @@ void GridVessel::setGridElement( const unsigned& ipoint, const unsigned& jelemen
setDataElement( 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 );
+ addDataElement( nper*ipoint + jelement, value );
+}
+
void GridVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
plumed_dbg_assert( myvals.getNumberOfValues()==(nper+1) );
for(unsigned i=0;i<nper;++i) buffer[bufstart + nper*current + i] += myvals.get(i+1);
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index a0aacf47c..6fdf49b73 100644
--- a/src/gridtools/GridVessel.h
+++ b/src/gridtools/GridVessel.h
@@ -89,10 +89,10 @@ public:
void getIndices( const unsigned& index, std::vector<unsigned>& indices ) const ;
/// Get the indices of a particular point
void getIndices( const std::vector<double>& point, std::vector<unsigned>& indices ) const ;
-
/// Operations on one of the elements of grid point i
void setGridElement( const unsigned&, const unsigned&, const double& );
-
+/// Add data to an element of the grid
+ void addToGridElement( const unsigned& ipoint, const unsigned& jelement, const double& value );
/// Operations on one of the elements of grid point specified by vector
double getGridElement( const std::vector<unsigned>&, const unsigned& ) const ;
void setGridElement( const std::vector<unsigned>&, const unsigned&, const double& );
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index b227ef693..50962b916 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/src/gridtools/HistogramOnGrid.cpp
@@ -102,7 +102,7 @@ void HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, st
KernelFunctions* kernel=getKernelAndNeighbors( point, num_neigh, neighbors );
if( !kernel ){
- plumed_dbg_assert( num_neigh==1 );
+ plumed_dbg_assert( num_neigh==1 ); der.resize(0);
accumulate( neighbors[0], weight, 1.0, der, buffer );
} else {
std::vector<Value*> vv( getVectorOfValues() );
--
GitLab