diff --git a/src/analysis/Histogram.cpp b/src/analysis/Histogram.cpp
index 16abb4a6dcb9ffa1791edd7a7721c8ca6dda1d3b..58e9815981cd83873586ca0022e8ddf60a4fdb1b 100644
--- a/src/analysis/Histogram.cpp
+++ b/src/analysis/Histogram.cpp
@@ -311,6 +311,8 @@ void Histogram::turnOnDerivatives() {
     if( !mbase ) error("do not know how to get histogram derivatives for actions of type " + myvessels[i]->getName() );
     tmp_atoms = mbase->getAbsoluteIndexes();
     for(unsigned j=0; j<tmp_atoms.size(); ++j) all_atoms.push_back( tmp_atoms[j] );
+    // Make a tempory multi value so we can avoid vector resizing
+    stashes[i]->resizeTemporyMultiValues( 1 );
   }
   ActionAtomistic::requestAtoms( all_atoms );
   finalForces.resize( 3*all_atoms.size() + 9 );
@@ -386,7 +388,10 @@ void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
     for(unsigned i=2; i<myvessels[0]->getNumberOfQuantities(); ++i) myvals.setValue( i-1, cvals[i] );
     myvals.setValue( 0, cvals[0] ); myvals.setValue( myvessels[0]->getNumberOfQuantities() - 1, ww );
     if( in_apply ) {
-      MultiValue tmpval( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
+      MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
+      if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
+          tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
+        tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
       stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), true, tmpval );
       for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
         unsigned jder=tmpval.getActiveIndex(j); myvals.addDerivative( 0, jder, tmpval.getDerivative(0, jder) );
@@ -406,7 +411,10 @@ void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
         stashes[j]->retrieveSequentialValue( current, false, cvals ); totweight *= cvals[0];
       }
       // And this bit the derivatives
-      MultiValue tmpval( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
+      MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
+      if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
+          tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
+        tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
       stashes[0]->retrieveDerivatives( stashes[0]->getTrueIndex(current), false, tmpval );
       for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
         unsigned jder=tmpval.getActiveIndex(j);
@@ -421,7 +429,10 @@ void Histogram::compute( const unsigned& current, MultiValue& myvals ) const {
       tnorm *= cvals[0]; myvals.setValue( 1+i, cvals[1] );
       // Get the derivatives as well if we are in apply
       if( in_apply ) {
-        MultiValue tmpval( myvessels[i]->getNumberOfQuantities(), myvessels[i]->getNumberOfDerivatives() );
+        MultiValue& tmpval = stashes[0]->getTemporyMultiValue(0);
+        if( tmpval.getNumberOfValues()!=myvessels[0]->getNumberOfQuantities() ||
+            tmpval.getNumberOfDerivatives()!=myvessels[0]->getNumberOfDerivatives() )
+          tmpval.resize( myvessels[0]->getNumberOfQuantities(), myvessels[0]->getNumberOfDerivatives() );
         stashes[i]->retrieveDerivatives( stashes[i]->getTrueIndex(current), false, tmpval );
         for(unsigned j=0; j<tmpval.getNumberActive(); ++j) {
           unsigned jder=tmpval.getActiveIndex(j);
diff --git a/src/gridtools/GridVessel.cpp b/src/gridtools/GridVessel.cpp
index 5282667a4c9c9146d20e64642ad8f7e92e70a2af..3b89af0f2d996e97431edaed5398f051291de031 100644
--- a/src/gridtools/GridVessel.cpp
+++ b/src/gridtools/GridVessel.cpp
@@ -191,9 +191,13 @@ void GridVessel::getIndices( const unsigned& index, std::vector<unsigned>& indic
 }
 
 void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<double>& x ) const {
-  plumed_dbg_assert( bounds_set && x.size()==dimension && ipoint<npoints );
+  std::vector<unsigned> tindices( dimension ); getGridPointCoordinates( ipoint, tindices, x );
+}
+
+void GridVessel::getGridPointCoordinates( const unsigned& ipoint, std::vector<unsigned>& tindices, std::vector<double>& x ) const {
+  plumed_dbg_assert( bounds_set && x.size()==dimension && tindices.size()==dimension && ipoint<npoints );
   if( gtype==flat ) {
-    std::vector<unsigned> tindices( dimension ); getIndices( ipoint, tindices );
+    getIndices( ipoint, tindices );
     for(unsigned i=0; i<dimension; ++i) x[i] = min[i] + dx[i]*tindices[i];
   } else if( gtype==fibonacci ) {
     x[1] = ((ipoint*fib_offset) - 1) + (fib_offset/2);
diff --git a/src/gridtools/GridVessel.h b/src/gridtools/GridVessel.h
index 4cb3c7beeb5c3e10c711749afd1a051172478987..120b9bba4981f352da9b16c54a39ada0ad3cf667 100644
--- a/src/gridtools/GridVessel.h
+++ b/src/gridtools/GridVessel.h
@@ -117,6 +117,7 @@ public:
   unsigned getNumberOfPoints() const;
 /// Get the coordinates for a point in the grid
   void getGridPointCoordinates( const unsigned&, std::vector<double>& ) const ;
+  void getGridPointCoordinates( const unsigned&, std::vector<unsigned>&, std::vector<double>& ) const ;
 /// Get the dimensionality of the function
   unsigned getDimension() const ;
 /// Get the number of components in the vector stored on each grid point
diff --git a/src/gridtools/HistogramOnGrid.cpp b/src/gridtools/HistogramOnGrid.cpp
index 43fe891c27f5b62d34fdfce21e1d489bb0605c0e..1ddae05a4c157e1eadc397e30f4b6a80b07c5b8f 100644
--- a/src/gridtools/HistogramOnGrid.cpp
+++ b/src/gridtools/HistogramOnGrid.cpp
@@ -121,11 +121,11 @@ void HistogramOnGrid::calculate( const unsigned& current, MultiValue& myvals, st
       std::vector<double> intforce( 2*dimension, 0.0 );
       std::vector<Value*> vv( getVectorOfValues() );
 
-      double newval; std::vector<double> xx( dimension );
+      double newval; std::vector<unsigned> tindices( dimension ); std::vector<double> xx( dimension );
       for(unsigned i=0; i<num_neigh; ++i) {
         unsigned ineigh=neighbors[i];
         if( inactive( ineigh ) ) continue ;
-        getGridPointCoordinates( ineigh, xx );
+        getGridPointCoordinates( ineigh, tindices, xx );
         if( kernel ) {
           for(unsigned j=0; j<dimension; ++j) vv[j]->set(xx[j]);
           newval = kernel->evaluate( vv, der, true );