From e23214d803dfa67aeaf4b8c7e7c8e9f85b20a2e4 Mon Sep 17 00:00:00 2001
From: Gareth Aneurin Tribello <garethtribello@Gareths-MacBook-Pro.local>
Date: Tue, 30 Jul 2024 10:19:38 +0100
Subject: [PATCH] Changed the code so that values are transferred directly to
 the Value rather than to the buffer

---
 .../p_entropies.xyz.reference                 |  4 +-
 src/colvar/MultiColvarTemplate.h              |  6 +-
 src/colvar/RMSDVector.cpp                     |  4 +-
 src/core/ActionWithMatrix.cpp                 | 93 ++++++-------------
 src/core/ActionWithMatrix.h                   | 16 +---
 src/core/ActionWithVector.cpp                 | 76 +++++++++------
 src/core/ActionWithVector.h                   | 14 +--
 src/core/Value.cpp                            |  4 -
 src/core/Value.h                              | 16 +---
 src/function/FunctionOfMatrix.h               |  6 +-
 src/function/FunctionOfVector.h               |  6 +-
 src/matrixtools/MatrixTimesVector.cpp         | 13 +--
 src/matrixtools/Voronoi.cpp                   |  3 +-
 src/refdist/MatrixProductDiagonal.cpp         |  2 +-
 src/tools/MultiValue.cpp                      | 10 +-
 src/tools/MultiValue.h                        | 27 +-----
 16 files changed, 118 insertions(+), 182 deletions(-)

diff --git a/regtest/gridtools/rt-pairentropies/p_entropies.xyz.reference b/regtest/gridtools/rt-pairentropies/p_entropies.xyz.reference
index 21d072bcc7..13527ae324 100644
--- a/regtest/gridtools/rt-pairentropies/p_entropies.xyz.reference
+++ b/regtest/gridtools/rt-pairentropies/p_entropies.xyz.reference
@@ -148,7 +148,7 @@ X    0.1056    1.6641   -0.0477   -1.8340
 X    0.7283    1.7274    0.7849   -1.1006
 X    0.9103    2.4770   -0.1044   -1.6791
 X    0.0688    2.5236    0.7617   -1.7107
-X   -0.1428    1.6741    1.8280   -2.2931
+X   -0.1428    1.6741    1.8280   -2.2930
 X    0.7031    1.5181    2.4387   -1.2983
 X    0.8450    2.4859    1.8490   -1.4526
 X    0.0682    2.5773    2.6448   -1.7259
@@ -238,7 +238,7 @@ X    1.5343   -0.0986   -0.1431   -1.1706
 X    2.4293   -0.2020    1.0139   -2.1738
 X    2.6885    0.7142    0.0622   -2.1087
 X    1.9617    0.7380    0.8137   -1.3226
-X    1.6673    0.0998    1.7622   -1.7197
+X    1.6673    0.0998    1.7622   -1.7196
 X    2.4352    0.0073    2.5561   -2.1955
 X    2.5243    0.9413    1.6602   -2.0844
 X    1.8709    0.8862    2.6985   -1.8413
diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h
index 3c15e05bec..27ce4eb512 100644
--- a/src/colvar/MultiColvarTemplate.h
+++ b/src/colvar/MultiColvarTemplate.h
@@ -44,7 +44,7 @@ class MultiColvarTemplate : public ActionWithVector {
   unsigned getNumberOfDerivatives() override ;
   void addValueWithDerivatives( const std::vector<unsigned>& shape=std::vector<unsigned>() ) override ;
   void addComponentWithDerivatives( const std::string& name, const std::vector<unsigned>& shape=std::vector<unsigned>() ) override ;
-  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) override ;
+  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) override ;
   void performTask( const unsigned&, MultiValue& ) const override ;
   void calculate() override;
 };
@@ -127,9 +127,9 @@ void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& nam
 }
 
 template <class T>
-void MultiColvarTemplate<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
+void MultiColvarTemplate<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) {
   if( wholemolecules ) makeWhole();
-  ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
+  ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol );
 }
 
 template <class T>
diff --git a/src/colvar/RMSDVector.cpp b/src/colvar/RMSDVector.cpp
index b183dd5e1e..2c214bdd99 100644
--- a/src/colvar/RMSDVector.cpp
+++ b/src/colvar/RMSDVector.cpp
@@ -269,9 +269,9 @@ void RMSDVector::gatherStoredValue( const unsigned& valindex, const unsigned& co
                                     const unsigned& bufstart, std::vector<double>& buffer ) const {
   if( getConstPntrToComponent(valindex)->getRank()==1 ) { ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer ); return; }
   const std::vector<Vector>& direction( myvals.getConstFirstAtomDerivativeVector()[1] );
-  unsigned natoms = direction.size(); unsigned vindex = bufstart + 3*code*natoms;
+  unsigned natoms = direction.size(); unsigned vindex = 3*code*natoms; Value* myval = const_cast<Value*>( getConstPntrToComponent(valindex) ); 
   for(unsigned i=0; i<natoms; ++i) {
-    for(unsigned j=0; j<3; ++j ) buffer[vindex + j*natoms + i] += direction[i][j];
+    for(unsigned j=0; j<3; ++j ) myval->set(vindex + j*natoms + i, direction[i][j] );
   }
 }
 
diff --git a/src/core/ActionWithMatrix.cpp b/src/core/ActionWithMatrix.cpp
index 7bdaeee5b8..a16ab86461 100644
--- a/src/core/ActionWithMatrix.cpp
+++ b/src/core/ActionWithMatrix.cpp
@@ -51,8 +51,8 @@ void ActionWithMatrix::getAllActionLabelsInMatrixChain( std::vector<std::string>
   if( matrix_to_do_after ) matrix_to_do_after->getAllActionLabelsInMatrixChain( mylabels );
 }
 
-void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
-  ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
+void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) {
+  ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol );
 
   for(int i=0; i<getNumberOfComponents(); ++i) {
     Value* myval=getPntrToComponent(i);
@@ -60,8 +60,6 @@ void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsi
     myval->setPositionInMatrixStash(nmat); nmat++;
     if( !myval->valueIsStored() ) continue;
     if( myval->getShape()[1]>maxcol ) maxcol=myval->getShape()[1];
-    myval->setMatrixBookeepingStart(nbookeeping);
-    nbookeeping += myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
   }
   // Turn off clearning of derivatives after each matrix run if there are no matrices in the output of this action
   clearOnEachCycle = false;
@@ -91,14 +89,13 @@ const ActionWithMatrix* ActionWithMatrix::getFirstMatrixInChain() const {
   return matrix_to_do_before->getFirstMatrixInChain();
 }
 
-void ActionWithMatrix::getTotalMatrixBookeeping( unsigned& nbookeeping ) {
+void ActionWithMatrix::setupMatrixStore() {
   for(int i=0; i<getNumberOfComponents(); ++i) {
     Value* myval=getPntrToComponent(i);
     if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) continue;
     myval->reshapeMatrixStore( getNumberOfColumns() );
-    nbookeeping += myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
   }
-  if( next_action_in_chain ) next_action_in_chain->getTotalMatrixBookeeping( nbookeeping );
+  if( next_action_in_chain ) next_action_in_chain->setupMatrixStore();
 }
 
 void ActionWithMatrix::calculate() {
@@ -106,9 +103,7 @@ void ActionWithMatrix::calculate() {
   // Update all the neighbour lists
   updateAllNeighbourLists();
   // Setup the matrix indices
-  unsigned nbookeeping=0; getTotalMatrixBookeeping( nbookeeping );
-  if( matrix_bookeeping.size()!=nbookeeping ) matrix_bookeeping.resize( nbookeeping );
-  std::fill( matrix_bookeeping.begin(), matrix_bookeeping.end(), 0 );
+  setupMatrixStore();
   // And run all the tasks
   runAllTasks();
 }
@@ -118,6 +113,16 @@ void ActionWithMatrix::updateAllNeighbourLists() {
   if( next_action_in_chain ) next_action_in_chain->updateAllNeighbourLists();
 }
 
+void ActionWithMatrix::clearBookeepingBeforeTask( const unsigned& task_index ) const {
+  // Reset the bookeeping elements for storage
+  for(unsigned i=0; i<getNumberOfComponents(); ++i) {
+      Value* myval = const_cast<Value*>( getConstPntrToComponent(i) ); unsigned ncols = myval->getNumberOfColumns();
+      if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() || ncols>=myval->getShape()[1] ) continue;
+      myval->matrix_bookeeping[task_index*(1+ncols)]=0;
+  }
+  if( matrix_to_do_after ) matrix_to_do_after->clearBookeepingBeforeTask( task_index );
+}
+
 void ActionWithMatrix::performTask( const unsigned& task_index, MultiValue& myvals ) const {
   std::vector<unsigned> & indices( myvals.getIndices() );
   if( matrix_to_do_before ) {
@@ -127,6 +132,9 @@ void ActionWithMatrix::performTask( const unsigned& task_index, MultiValue& myva
   }
   setupForTask( task_index, indices, myvals );
 
+  // Reset the bookeeping elements for storage
+  clearBookeepingBeforeTask( task_index );
+
   // Now loop over the row of the matrix
   unsigned ntwo_atoms = myvals.getSplitIndex();
   for(unsigned i=1; i<ntwo_atoms; ++i) {
@@ -149,72 +157,31 @@ void ActionWithMatrix::runTask( const std::string& controller, const unsigned& c
 
   if( hasval ) {
     for(int i=0; i<getNumberOfComponents(); ++i) {
-      const Value* myval=getConstPntrToComponent(i);
+      const Value* myval=getConstPntrToComponent(i); unsigned ncols = myval->getNumberOfColumns();
       if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() ) continue;
-      unsigned matindex = myval->getPositionInMatrixStash(), matbook_start = myval->getMatrixBookeepingStart(), col_stash_index = colno;
+      unsigned matindex = myval->getPositionInMatrixStash(), col_stash_index = colno;
       if( colno>=myval->getShape()[0] ) col_stash_index = colno - myval->getShape()[0];
-      unsigned rowstart = matbook_start+current*(1+myval->getNumberOfColumns());
       if( myval->forcesWereAdded() ) {
-        unsigned sind = myval->getPositionInStream(), find = myvals.getMatrixBookeeping()[rowstart];
-        double fforce = myval->getForce( myvals.getTaskIndex()*myval->getNumberOfColumns() + find );
-        if( getNumberOfColumns()>=myval->getShape()[1] ) fforce = myval->getForce( myvals.getTaskIndex()*myval->getShape()[1] + col_stash_index );
+        unsigned sind = myval->getPositionInStream();
+        double fforce = myval->getForce( myvals.getTaskIndex()*ncols + myval->matrix_bookeeping[current*(1+ncols)] );
+        if( ncols>=myval->getShape()[1] ) fforce = myval->getForce( myvals.getTaskIndex()*myval->getShape()[1] + col_stash_index );
         for(unsigned j=0; j<myvals.getNumberActive(sind); ++j) {
           unsigned kindex = myvals.getActiveIndex(sind,j); myvals.addMatrixForce( matindex, kindex, fforce*myvals.getDerivative(sind,kindex ) );
         }
       }
       double finalval = myvals.get( myval->getPositionInStream() );
-      if( fabs(finalval)>0 ) myvals.stashMatrixElement( matindex, rowstart, col_stash_index, finalval );
+      if( fabs(finalval)>0 ) {
+          Value* myv = const_cast<Value*>( myval );
+          if( ncols<myval->getShape()[1] ) { 
+              myv->set( current*ncols + myval->matrix_bookeeping[current*(1+ncols)], finalval );
+              myv->matrix_bookeeping[current*(1+ncols)]++; myv->matrix_bookeeping[current*(1+ncols)+myval->matrix_bookeeping[current*(1+ncols)]] = col_stash_index;
+          } else myv->set( current*myval->getShape()[1] + col_stash_index, finalval );
+      }
     }
   }
   if( matrix_to_do_after ) matrix_to_do_after->runTask( controller, current, colno, myvals );
 }
 
-void ActionWithMatrix::gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) {
-  ActionWithVector::gatherThreads( nt, bufsize, omp_buffer, buffer, myvals );
-  for(unsigned i=0; i<matrix_bookeeping.size(); ++i) matrix_bookeeping[i] += myvals.getMatrixBookeeping()[i];
-}
-
-void ActionWithMatrix::gatherProcesses( std::vector<double>& buffer ) {
-  ActionWithVector::gatherProcesses( buffer );
-  if( matrix_bookeeping.size()>0 && !runInSerial() ) comm.Sum( matrix_bookeeping );
-  unsigned nval=0; transferNonZeroMatrixElementsToValues( nval, matrix_bookeeping );
-}
-
-void ActionWithMatrix::transferNonZeroMatrixElementsToValues( unsigned& nval, const std::vector<unsigned>& matbook ) {
-  for(int i=0; i<getNumberOfComponents(); ++i) {
-    Value* myval=getPntrToComponent(i);
-    if( myval->getRank()!=2 || myval->hasDerivatives() || !myval->valueIsStored() || getNumberOfColumns()>=myval->getShape()[1] ) continue;
-    unsigned nelements = myval->getShape()[0]*( 1 + myval->getNumberOfColumns() );
-    for(unsigned j=0; j<nelements; ++j) myval->setMatrixBookeepingElement( j, matbook[nval+j] );
-    nval += nelements;
-  }
-  if( next_action_in_chain ) next_action_in_chain->transferNonZeroMatrixElementsToValues( nval, matbook );
-}
-
-void ActionWithMatrix::gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals,
-    const unsigned& bufstart, std::vector<double>& buffer ) const {
-  if( getConstPntrToComponent(valindex)->getRank()==1 ) { ActionWithVector::gatherStoredValue( valindex, code, myvals, bufstart, buffer ); return; }
-  const Value* myval=getConstPntrToComponent(valindex);
-  unsigned ncols = myval->getNumberOfColumns(), matind = myval->getPositionInMatrixStash();
-  unsigned matbook_start = myval->getMatrixBookeepingStart(), vindex = bufstart + code*myval->getNumberOfColumns();
-  const std::vector<unsigned> & matbook( myvals.getMatrixBookeeping() ); unsigned nelements = matbook[matbook_start+code*(1+ncols)];
-  if( ncols>=myval->getShape()[1] ) {
-    // In this case we store the full matrix
-    for(unsigned j=0; j<nelements; ++j) {
-      unsigned jind = matbook[matbook_start+code*(1+ncols)+1+j];
-      plumed_dbg_massert( vindex+j<buffer.size(), "failing in " + getLabel() + " on value " + myval->getName() );
-      buffer[vindex + jind] += myvals.getStashedMatrixElement( matind, jind );
-    }
-  } else {
-    // This is for storing sparse matrices when we can
-    for(unsigned j=0; j<nelements; ++j) {
-      unsigned jind = matbook[matbook_start+code*(1+ncols)+1+j];
-      plumed_dbg_massert( vindex+j<buffer.size(), "failing in " + getLabel() + " on value " + myval->getName() );
-      buffer[vindex + j] += myvals.getStashedMatrixElement( matind, jind );
-    }
-  }
-}
-
 bool ActionWithMatrix::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
   if( myval->getRank()<2 ) return ActionWithVector::checkForTaskForce( itask, myval );
   unsigned nelements = myval->getRowLength(itask), startr = itask*myval->getNumberOfColumns();
diff --git a/src/core/ActionWithMatrix.h b/src/core/ActionWithMatrix.h
index 00a96ac10c..b5a870ba2c 100644
--- a/src/core/ActionWithMatrix.h
+++ b/src/core/ActionWithMatrix.h
@@ -31,16 +31,14 @@ class ActionWithMatrix : public ActionWithVector {
   ActionWithMatrix* next_action_in_chain;
   ActionWithMatrix* matrix_to_do_before;
   ActionWithMatrix* matrix_to_do_after;
-/// This holds the bookeeping arrays for sparse matrices
-  std::vector<unsigned> matrix_bookeeping;
 /// Update all the neighbour lists in the chain
   void updateAllNeighbourLists();
+/// This clears all bookeeping arrays before the ith task
+  void clearBookeepingBeforeTask( const unsigned& task_index ) const ;
 /// This is used to clear up the matrix elements
   void clearMatrixElements( MultiValue& myvals ) const ;
 /// This is used to find the total amount of space we need for storing matrix elements
-  void getTotalMatrixBookeeping( unsigned& stashsize );
-/// This transfers the non-zero elements to the Value
-  void transferNonZeroMatrixElementsToValues( unsigned& nval, const std::vector<unsigned>& matbook );
+  void setupMatrixStore();
 /// This does the calculation of a particular matrix element
   void runTask( const std::string& controller, const unsigned& current, const unsigned colno, MultiValue& myvals ) const ;
 protected:
@@ -71,17 +69,11 @@ class ActionWithMatrix : public ActionWithVector {
 /// This should return the number of columns to help with sparse storage of matrices
   virtual unsigned getNumberOfColumns() const = 0;
 /// This requires some thought
-  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) override;
+  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) override;
 //// This does some setup before we run over the row of the matrix
   virtual void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const = 0;
 /// Run over one row of the matrix
   virtual void performTask( const unsigned& task_index, MultiValue& myvals ) const ;
-/// Gather a row of the matrix
-  void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, const unsigned& bufstart, std::vector<double>& buffer ) const override;
-/// Gather all the data from the threads
-  void gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) override ;
-/// Gather all the data from the MPI processes
-  void gatherProcesses( std::vector<double>& buffer ) override;
 /// This is the virtual that will do the calculation of the task for a particular matrix element
   virtual void performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const = 0;
 /// This is the jobs that need to be done when we have run all the jobs in a row of the matrix
diff --git a/src/core/ActionWithVector.cpp b/src/core/ActionWithVector.cpp
index fa19b58ca2..47942aff98 100644
--- a/src/core/ActionWithVector.cpp
+++ b/src/core/ActionWithVector.cpp
@@ -457,6 +457,16 @@ bool ActionWithVector::doNotCalculateDerivatives() const {
   return ActionWithValue::doNotCalculateDerivatives();
 }
 
+void ActionWithVector::clearMatrixBookeeping() {
+  for(unsigned i=0; i<getNumberOfComponents(); ++i) {
+      Value* myval = getPntrToComponent(i); if( !myval->storedata ) continue;
+      if( myval->getRank()==2 && myval->getNumberOfColumns()<myval->getShape()[1] ) {
+          std::fill(myval->matrix_bookeeping.begin(), myval->matrix_bookeeping.end(), 0);
+      }
+      myval->set(0);
+  }
+}
+
 void ActionWithVector::runAllTasks() {
 // Skip this if this is done elsewhere
   if( action_to_do_before ) return;
@@ -465,6 +475,9 @@ void ActionWithVector::runAllTasks() {
   unsigned rank=comm.Get_rank();
   if(serial) { stride=1; rank=0; }
 
+  // Clear matrix bookeeping arrays
+  if( stride>1 ) clearMatrixBookeeping();
+
   // Get the list of active tasks
   std::vector<unsigned> & partialTaskList( getListOfActiveTasks( this ) );
   unsigned nactive_tasks=partialTaskList.size();
@@ -484,8 +497,8 @@ void ActionWithVector::runAllTasks() {
     }
   }
   // Get the total number of streamed quantities that we need
-  unsigned nquants=0, nmatrices=0, maxcol=0, nbooks=0;
-  getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol, nbooks );
+  unsigned nquants=0, nmatrices=0, maxcol=0;
+  getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol );
   // Get size for buffer
   unsigned bufsize=0; getSizeOfBuffer( nactive_tasks, bufsize );
   if( buffer.size()!=bufsize ) buffer.resize( bufsize );
@@ -500,7 +513,7 @@ void ActionWithVector::runAllTasks() {
   {
     std::vector<double> omp_buffer;
     if( nt>1 ) omp_buffer.resize( bufsize, 0.0 );
-    MultiValue myvals( nquants, nderivatives, nmatrices, maxcol, nbooks );
+    MultiValue myvals( nquants, nderivatives, nmatrices, maxcol );
     myvals.clearAll();
 
     #pragma omp for nowait
@@ -516,20 +529,25 @@ void ActionWithVector::runAllTasks() {
       myvals.clearAll();
     }
     #pragma omp critical
-    gatherThreads( nt, bufsize, omp_buffer, buffer, myvals );
+    if( nt>1 ) for(unsigned i=0; i<bufsize; ++i) buffer[i]+=omp_buffer[i];
   }
 
   // MPI Gather everything
-  if( !serial && buffer.size()>0 ) gatherProcesses( buffer );
+  if( !serial ) { 
+    if( buffer.size()>0 ) comm.Sum( buffer ); 
+    gatherProcesses(); 
+  }
   finishComputations( buffer ); forwardPass=false;
 }
 
-void ActionWithVector::gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) {
-  if( nt>1 ) for(unsigned i=0; i<bufsize; ++i) buffer[i]+=omp_buffer[i];
-}
-
-void ActionWithVector::gatherProcesses( std::vector<double>& buffer ) {
-  comm.Sum( buffer );
+void ActionWithVector::gatherProcesses() {
+  for(unsigned i=0; i<getNumberOfComponents(); ++i) {
+      Value* myval = getPntrToComponent(i);
+      if( myval->storedata && !myval->hasDeriv ) {
+          comm.Sum( myval->data ); if( myval->getRank()==2 && myval->getNumberOfColumns()<myval->getShape()[1] ) comm.Sum( myval->matrix_bookeeping );
+      }
+  }
+  if( action_to_do_after ) action_to_do_after->gatherProcesses();
 }
 
 bool ActionWithVector::checkForGrids( unsigned& nder ) const {
@@ -564,20 +582,23 @@ void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) {
   if( action_to_do_after ) action_to_do_after->getNumberOfTasks( ntasks );
 }
 
-void ActionWithVector::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
+void ActionWithVector::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) {
   for(unsigned i=0; i<getNumberOfArguments(); ++i) {
     if( !getPntrToArgument(i)->ignoreStoredValue(headstr) ) { getPntrToArgument(i)->streampos=nquants; nquants++; }
   }
   for(int i=0; i<getNumberOfComponents(); ++i) { getPntrToComponent(i)->streampos=nquants; nquants++; }
 }
 
-void ActionWithVector::getNumberOfStreamedQuantities( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
-  setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
-  if( action_to_do_after ) action_to_do_after->getNumberOfStreamedQuantities( headstr, nquants, nmat, maxcol, nbookeeping );
+void ActionWithVector::getNumberOfStreamedQuantities( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) {
+  setupStreamedComponents( headstr, nquants, nmat, maxcol );
+  if( action_to_do_after ) action_to_do_after->getNumberOfStreamedQuantities( headstr, nquants, nmat, maxcol );
 }
 
 void ActionWithVector::getSizeOfBuffer( const unsigned& nactive_tasks, unsigned& bufsize ) {
-  for(int i=0; i<getNumberOfComponents(); ++i) { getPntrToComponent(i)->bufstart=bufsize; bufsize += getPntrToComponent(i)->data.size(); }
+  for(int i=0; i<getNumberOfComponents(); ++i) { 
+      getPntrToComponent(i)->bufstart=bufsize; 
+      if( getPntrToComponent(i)->hasDerivatives() || getPntrToComponent(i)->getRank()==0 ) bufsize += getPntrToComponent(i)->data.size();
+  }
   if( action_to_do_after ) action_to_do_after->getSizeOfBuffer( nactive_tasks, bufsize );
 }
 
@@ -618,6 +639,11 @@ bool ActionWithVector::getNumberOfStoredValues( Value* startat, unsigned& nvals,
 
 void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) const {
   myvals.setTaskIndex(current); myvals.vector_call=true; performTask( current, myvals );
+  for(unsigned i=0; i<getNumberOfComponents(); ++i) {
+      const Value* myval = getConstPntrToComponent(i); 
+      if( myval->getRank()!=1 || myval->hasDerivatives() || !myval->valueIsStored() ) continue;
+      Value* myv = const_cast<Value*>( myval ); myv->set( current, myvals.get( myval->getPositionInStream() ) );
+  }
   if( action_to_do_after ) action_to_do_after->runTask( current, myvals );
 }
 
@@ -641,20 +667,12 @@ void ActionWithVector::gatherAccumulators( const unsigned& taskCode, const Multi
   if( action_to_do_after ) action_to_do_after->gatherAccumulators( taskCode, myvals, buffer );
 }
 
-void ActionWithVector::gatherStoredValue( const unsigned& valindex, const unsigned& taskCode, const MultiValue& myvals,
-    const unsigned& bufstart, std::vector<double>& buffer ) const {
-  plumed_dbg_assert( getConstPntrToComponent(valindex)->getRank()==1 && !getConstPntrToComponent(valindex)->hasDeriv );
-  unsigned vindex = getConstPntrToComponent(valindex)->bufstart + taskCode; plumed_dbg_massert( vindex<buffer.size(), "failing in " + getLabel() );
-  buffer[vindex] += myvals.get(getConstPntrToComponent(valindex)->streampos);
-}
-
 void ActionWithVector::finishComputations( const std::vector<double>& buf ) {
   for(int i=0; i<getNumberOfComponents(); ++i) {
     // This gathers vectors and grids at the end of the calculation
     unsigned bufstart = getPntrToComponent(i)->bufstart;
-    getPntrToComponent(i)->data.assign( getPntrToComponent(i)->data.size(), 0 );
-    if( (getPntrToComponent(i)->getRank()>0 && getPntrToComponent(i)->hasDerivatives()) || getPntrToComponent(i)->storedata ) {
-      unsigned sz_v = getPntrToComponent(i)->data.size();
+    if( (getPntrToComponent(i)->getRank()>0 && getPntrToComponent(i)->hasDerivatives()) ) {
+      unsigned sz_v = getPntrToComponent(i)->data.size(); getPntrToComponent(i)->data.assign( getPntrToComponent(i)->data.size(), 0 ); 
       for(unsigned j=0; j<sz_v; ++j) {
         plumed_dbg_assert( bufstart+j<buf.size() );
         getPntrToComponent(i)->add( j, buf[bufstart+j] );
@@ -699,8 +717,8 @@ bool ActionWithVector::checkForForces() {
   if( nt==0 ) nt=1;
 
   // Now determine how big the multivalue needs to be
-  unsigned nquants=0, nmatrices=0, maxcol=0, nbooks=0;
-  getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol, nbooks );
+  unsigned nquants=0, nmatrices=0, maxcol=0;
+  getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol );
   // Recover the number of derivatives we require (this should be equal to the number of forces)
   unsigned nderiv=0; getNumberOfStreamedDerivatives( nderiv, NULL );
   if( !action_to_do_after && arg_deriv_starts.size()>0 ) {
@@ -717,7 +735,7 @@ bool ActionWithVector::checkForForces() {
   {
     std::vector<double> omp_forces;
     if( nt>1 ) omp_forces.resize( forcesForApply.size(), 0.0 );
-    MultiValue myvals( nquants, nderiv, nmatrices, maxcol, nbooks );
+    MultiValue myvals( nquants, nderiv, nmatrices, maxcol );
     myvals.clearAll();
 
     #pragma omp for nowait
diff --git a/src/core/ActionWithVector.h b/src/core/ActionWithVector.h
index 4b8676ebd0..6c4f0036de 100644
--- a/src/core/ActionWithVector.h
+++ b/src/core/ActionWithVector.h
@@ -53,6 +53,10 @@ class ActionWithVector:
   std::vector<ActionWithVector*> task_control_list;
 /// Work backwards through the chain to find an action that has either stored arguments or derivatives
   ActionWithVector* getActionWithDerivatives( ActionWithVector* depaction );
+/// Gather all the data from the MPI processes
+  void gatherProcesses();
+/// Clear all the bookeeping arrays
+  void clearMatrixBookeeping();
 /// Check if there are any grids in the stream
   bool checkForGrids(unsigned& nder) const ;
 ///  Run the task
@@ -64,7 +68,7 @@ class ActionWithVector:
 /// Get the size of the buffer array that holds the data we are gathering over the MPI loop
   void getSizeOfBuffer( const unsigned& nactive_tasks, unsigned& bufsize );
 /// Get the number of quantities in the stream
-  void getNumberOfStreamedQuantities( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping );
+  void getNumberOfStreamedQuantities( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol );
 /// Get the number of stored values in the stream
   bool getNumberOfStoredValues( Value* startat, unsigned& nvals, const unsigned& astart, const std::vector<Value*>& stopat );
 /// Add this action to the recursive chain
@@ -155,19 +159,15 @@ class ActionWithVector:
 /// Get the additional tasks that are required here
   virtual void getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks );
 /// setup the streamed quantities
-  virtual void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping );
+  virtual void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol );
 /// This we override to perform each individual task
   virtual void performTask( const unsigned& current, MultiValue& myvals ) const = 0;
 /// This is used to ensure that all indices are updated when you do local average
   virtual void updateAdditionalIndices( const unsigned& ostrn, MultiValue& myvals ) const {}
-/// Gather the data from all the OpenMP threads
-  virtual void gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals );
 /// Can be used to reduce the number of tasks that are performed when you use an ation from elsewhere
   virtual void switchTaskReduction( const bool& task_reduction, ActionWithVector* aselect ) {}
-/// Gather all the data from the MPI processes
-  virtual void gatherProcesses( std::vector<double>& buffer );
 /// Gather the values that we intend to store in the buffer
-  virtual void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, const unsigned& bufstart, std::vector<double>& buffer ) const ;
+  virtual void gatherStoredValue( const unsigned& valindex, const unsigned& code, const MultiValue& myvals, const unsigned& bufstart, std::vector<double>& buffer ) const {}
 /// Get the force tasks that are active for this action
   virtual void updateForceTasksFromValue( const Value* myval, std::vector<unsigned>& force_tasks ) const ;
 /// Check if there is a force that needs to be accumulated on the ith task
diff --git a/src/core/Value.cpp b/src/core/Value.cpp
index 6f4466bc3f..5cd438cf89 100644
--- a/src/core/Value.cpp
+++ b/src/core/Value.cpp
@@ -44,7 +44,6 @@ Value::Value():
   matpos(0),
   ngrid_der(0),
   ncols(0),
-  book_start(0),
   symmetric(false),
   periodicity(unset),
   min(0.0),
@@ -69,7 +68,6 @@ Value::Value(const std::string& name):
   ngrid_der(0),
   matpos(0),
   ncols(0),
-  book_start(0),
   symmetric(false),
   periodicity(unset),
   min(0.0),
@@ -94,7 +92,6 @@ Value::Value(ActionWithValue* av, const std::string& name, const bool withderiv,
   ngrid_der(0),
   matpos(0),
   ncols(0),
-  book_start(0),
   symmetric(false),
   periodicity(unset),
   min(0.0),
@@ -307,7 +304,6 @@ void Value::reshapeMatrixStore( const unsigned& n ) {
       }
     }
   }
-  if( ncols<shape[1] ) std::fill(matrix_bookeeping.begin(), matrix_bookeeping.end(), 0);
 }
 
 void Value::setPositionInMatrixStash( const unsigned& p ) {
diff --git a/src/core/Value.h b/src/core/Value.h
index 20cbfeeee6..edc550d866 100644
--- a/src/core/Value.h
+++ b/src/core/Value.h
@@ -48,6 +48,7 @@ class ActionAtomistic;
 class Value {
   friend class ActionWithValue;
   friend class ActionWithVector;
+  friend class ActionWithMatrix;
   friend class ActionAtomistic;
   friend class ActionWithArguments;
   friend class ActionWithVirtualAtom;
@@ -86,7 +87,7 @@ class Value {
 /// Does this quanity have derivatives
   bool hasDeriv;
 /// Variables for storing data
-  unsigned bufstart, streampos, matpos, ngrid_der, ncols, book_start;
+  unsigned bufstart, streampos, matpos, ngrid_der, ncols;
 /// If we are storing a matrix is it symmetric?
   bool symmetric;
 /// This is a bookeeping array that holds the non-zero elements of the "sparse" matrix
@@ -213,9 +214,6 @@ class Value {
 /// This stuff handles where to look for the start of the row that contains the row of the matrix
   void setPositionInMatrixStash( const unsigned& p );
   unsigned getPositionInMatrixStash() const ;
-/// This stuff handles where to keep the bookeeping stuff for storing the sparse matrix
-  void setMatrixBookeepingStart( const unsigned& b );
-  unsigned getMatrixBookeepingStart() const ;
 /// Convert the input index to its corresponding indices
   void convertIndexToindices(const std::size_t& index, std::vector<unsigned>& indices ) const ;
 /// Print out all the values in this Value
@@ -436,16 +434,6 @@ unsigned Value::getPositionInMatrixStash() const {
   return matpos;
 }
 
-inline
-void Value::setMatrixBookeepingStart( const unsigned& b ) {
-  book_start = b;
-}
-
-inline
-unsigned Value::getMatrixBookeepingStart() const {
-  return book_start;
-}
-
 inline
 void Value::setMatrixBookeepingElement( const unsigned& i, const unsigned& n ) {
   plumed_dbg_assert( i<matrix_bookeeping.size() );
diff --git a/src/function/FunctionOfMatrix.h b/src/function/FunctionOfMatrix.h
index 067e83c4b8..f1a68e86b9 100644
--- a/src/function/FunctionOfMatrix.h
+++ b/src/function/FunctionOfMatrix.h
@@ -63,7 +63,7 @@ class FunctionOfMatrix : public ActionWithMatrix {
 /// This checks for tasks in the parent class
 //  void buildTaskListFromArgumentRequests( const unsigned& ntasks, bool& reduce, std::set<AtomNumber>& otasks ) override ;
 /// This ensures that we create some bookeeping stuff during the first step
-  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) override ;
+  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) override ;
 /// This sets up for the task
   void setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const ;
 /// Calculate the full matrix
@@ -248,7 +248,7 @@ void FunctionOfMatrix<T>::setupForTask( const unsigned& task_index, std::vector<
 // }
 
 template <class T>
-void FunctionOfMatrix<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
+void FunctionOfMatrix<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) {
   if( firststep ) {
     stored_arguments.resize( getNumberOfArguments() );
     update_arguments.resize( getNumberOfArguments(), true );
@@ -260,7 +260,7 @@ void FunctionOfMatrix<T>::setupStreamedComponents( const std::string& headstr, u
     }
     firststep=false;
   }
-  ActionWithMatrix::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
+  ActionWithMatrix::setupStreamedComponents( headstr, nquants, nmat, maxcol );
 }
 
 template <class T>
diff --git a/src/function/FunctionOfVector.h b/src/function/FunctionOfVector.h
index bb962cf924..98780b8f3b 100644
--- a/src/function/FunctionOfVector.h
+++ b/src/function/FunctionOfVector.h
@@ -67,7 +67,7 @@ class FunctionOfVector : public ActionWithVector {
 /// This builds the task list for the action
   void calculate() override;
 /// This ensures that we create some bookeeping stuff during the first step
-  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) override ;
+  void setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) override ;
 /// Calculate the function
   void performTask( const unsigned& current, MultiValue& myvals ) const override ;
 };
@@ -206,7 +206,7 @@ void FunctionOfVector<T>::prepare() {
 }
 
 template <class T>
-void FunctionOfVector<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
+void FunctionOfVector<T>::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol ) {
   if( firststep ) {
     stored_arguments.resize( getNumberOfArguments() );
     std::string control = getFirstActionInChain()->getLabel();
@@ -216,7 +216,7 @@ void FunctionOfVector<T>::setupStreamedComponents( const std::string& headstr, u
     }
     firststep=false;
   }
-  ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
+  ActionWithVector::setupStreamedComponents( headstr, nquants, nmat, maxcol );
 }
 
 template <class T>
diff --git a/src/matrixtools/MatrixTimesVector.cpp b/src/matrixtools/MatrixTimesVector.cpp
index 6e54bf65f3..ec99a4262a 100644
--- a/src/matrixtools/MatrixTimesVector.cpp
+++ b/src/matrixtools/MatrixTimesVector.cpp
@@ -249,23 +249,24 @@ void MatrixTimesVector::setupForTask( const unsigned& task_index, std::vector<un
 }
 
 void MatrixTimesVector::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
+  std::string headstr=getFirstActionInChain()->getLabel();
   unsigned ind2 = index2; if( index2>=getPntrToArgument(0)->getShape()[0] ) ind2 = index2 - getPntrToArgument(0)->getShape()[0];
   if( sumrows ) {
     unsigned n=getNumberOfArguments()-1; double matval = 0;
     for(unsigned i=0; i<getNumberOfArguments()-1; ++i) {
       unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
       Value* myarg = getPntrToArgument(i);
-      if( !myarg->valueHasBeenSet() ) myvals.addValue( ostrn, myvals.get( myarg->getPositionInStream() ) );
-      else myvals.addValue( ostrn, myarg->get( index1*myarg->getNumberOfColumns() + ind2, false ) );
+      if( !myarg->valueHasBeenSet() || myarg->ignoreStoredValue(headstr) ) myvals.addValue( ostrn, myvals.get( myarg->getPositionInStream() ) );
+      else myvals.addValue( ostrn, myarg->get( index1*myarg->getNumberOfColumns() + ind2, true ) );
       // Now lets work out the derivatives
       if( doNotCalculateDerivatives() ) continue;
       addDerivativeOnMatrixArgument( stored_arg[i], i, i, index1, ind2, 1.0, myvals );
     }
   } else if( getPntrToArgument(1)->getRank()==1 ) {
     double matval = 0; Value* myarg = getPntrToArgument(0); unsigned vcol = ind2;
-    if( !myarg->valueHasBeenSet() ) matval = myvals.get( myarg->getPositionInStream() );
+    if( !myarg->valueHasBeenSet() || myarg->ignoreStoredValue(headstr) ) matval = myvals.get( myarg->getPositionInStream() );
     else {
-      matval = myarg->get( index1*myarg->getNumberOfColumns() + ind2, false );
+      matval = myarg->get( index1*myarg->getNumberOfColumns() + ind2, true );
       vcol = getPntrToArgument(0)->getRowIndex( index1, ind2 );
     }
     for(unsigned i=0; i<getNumberOfArguments()-1; ++i) {
@@ -282,9 +283,9 @@ void MatrixTimesVector::performTask( const std::string& controller, const unsign
     for(unsigned i=0; i<getNumberOfArguments()-1; ++i) {
       unsigned ostrn = getConstPntrToComponent(i)->getPositionInStream();
       Value* myarg = getPntrToArgument(i);
-      if( !myarg->valueHasBeenSet() ) matval = myvals.get( myarg->getPositionInStream() );
+      if( !myarg->valueHasBeenSet() || myarg->ignoreStoredValue(headstr) ) matval = myvals.get( myarg->getPositionInStream() );
       else {
-        matval = myarg->get( index1*myarg->getNumberOfColumns() + ind2, false );
+        matval = myarg->get( index1*myarg->getNumberOfColumns() + ind2, true );
         vcol = getPntrToArgument(i)->getRowIndex( index1, ind2 );
       }
       double vecval=getArgumentElement( n, vcol, myvals );
diff --git a/src/matrixtools/Voronoi.cpp b/src/matrixtools/Voronoi.cpp
index a5ab13090b..64083c2c96 100644
--- a/src/matrixtools/Voronoi.cpp
+++ b/src/matrixtools/Voronoi.cpp
@@ -79,7 +79,8 @@ void Voronoi::gatherStoredValue( const unsigned& valindex, const unsigned& code,
     double value = arg0->get( code*arg0->getShape()[1] + i );
     if( value<minmax ) { minmax = value; nv = i; }
   }
-  buffer[bufstart + code*arg0->getShape()[1] + nv] = 1;
+  Value* myval = const_cast<Value*>( getConstPntrToComponent(0) );
+  myval->set( code*arg0->getShape()[1] + nv, 1 ); 
 }
 
 }
diff --git a/src/refdist/MatrixProductDiagonal.cpp b/src/refdist/MatrixProductDiagonal.cpp
index 74e32fc7fd..e52e535d9c 100644
--- a/src/refdist/MatrixProductDiagonal.cpp
+++ b/src/refdist/MatrixProductDiagonal.cpp
@@ -129,7 +129,7 @@ void MatrixProductDiagonal::performTask( const unsigned& task_index, MultiValue&
 void MatrixProductDiagonal::calculate() {
   if( getPntrToArgument(1)->getRank()==1 ) {
     unsigned nder = getNumberOfDerivatives();
-    MultiValue myvals( 1, nder, 0, 0, 0 ); performTask( 0, myvals );
+    MultiValue myvals( 1, nder, 0, 0 ); performTask( 0, myvals );
 
     Value* myval=getPntrToComponent(0); myval->set( myvals.get(0) );
     for(unsigned i=0; i<nder; ++i) myval->setDerivative( i, myvals.getDerivative(0,i) );
diff --git a/src/tools/MultiValue.cpp b/src/tools/MultiValue.cpp
index c374a5d286..2e07d58385 100644
--- a/src/tools/MultiValue.cpp
+++ b/src/tools/MultiValue.cpp
@@ -24,7 +24,7 @@
 
 namespace PLMD {
 
-MultiValue::MultiValue( const size_t& nvals, const size_t& nder, const size_t& nmat, const size_t& maxcol, const size_t& nbook ):
+MultiValue::MultiValue( const size_t& nvals, const size_t& nder, const size_t& nmat, const size_t& maxcol ):
   task_index(0),
   task2_index(0),
   values(nvals),
@@ -40,9 +40,7 @@ MultiValue::MultiValue( const size_t& nvals, const size_t& nder, const size_t& n
   nindices(0),
   nsplit(0),
   nmatrix_cols(maxcol),
-  matrix_row_stash(nmat*maxcol,0),
   matrix_force_stash(nder*nmat),
-  matrix_bookeeping(nbook,0),
   matrix_row_nderivatives(nmat,0),
   matrix_row_derivative_indices(nmat)
 {
@@ -52,10 +50,10 @@ MultiValue::MultiValue( const size_t& nvals, const size_t& nder, const size_t& n
   for(unsigned i=0; i<nder; ++i) myind[i]=i;
 }
 
-void MultiValue::resize( const size_t& nvals, const size_t& nder, const size_t& nmat, const size_t& maxcol, const size_t& nbook ) {
+void MultiValue::resize( const size_t& nvals, const size_t& nder, const size_t& nmat, const size_t& maxcol ) {
   values.resize(nvals); nderivatives=nder; derivatives.resize( nvals*nder );
   hasderiv.resize(nvals*nder,false); nactive.resize(nvals); active_list.resize(nvals*nder);
-  nmatrix_cols=maxcol; matrix_row_stash.resize(nmat*maxcol,0); matrix_force_stash.resize(nmat*nder,0); matrix_bookeeping.resize(nbook, 0);
+  nmatrix_cols=maxcol; matrix_force_stash.resize(nmat*nder,0); 
   matrix_row_nderivatives.resize(nmat,0); matrix_row_derivative_indices.resize(nmat); atLeastOneSet=false;
   for(unsigned i=0; i<nmat; ++i) matrix_row_derivative_indices[i].resize( nder );
   // All crap from here onwards
@@ -65,8 +63,6 @@ void MultiValue::resize( const size_t& nvals, const size_t& nder, const size_t&
 
 void MultiValue::clearAll() {
   for(unsigned i=0; i<values.size(); ++i) values[i]=0;
-  // Clear matrix row
-  std::fill( matrix_row_stash.begin(), matrix_row_stash.end(), 0 );
   // Clear matrix derivative indices
   std::fill( matrix_row_nderivatives.begin(), matrix_row_nderivatives.end(), 0 );
   // Clear matrix forces
diff --git a/src/tools/MultiValue.h b/src/tools/MultiValue.h
index 4e786d2fb0..0fc32c7451 100644
--- a/src/tools/MultiValue.h
+++ b/src/tools/MultiValue.h
@@ -56,9 +56,7 @@ class MultiValue {
   unsigned nindices, nsplit;
 /// This allows us to store matrix elements
   unsigned nmatrix_cols;
-  std::vector<double> matrix_row_stash;
   std::vector<double> matrix_force_stash;
-  std::vector<unsigned> matrix_bookeeping;
 /// These are used to store the indices that have derivatives wrt to at least one
 /// of the elements in a matrix
   std::vector<unsigned> matrix_row_nderivatives;
@@ -70,8 +68,8 @@ class MultiValue {
   std::vector<Tensor> tmp_atom_virial;
   std::vector<std::vector<double> > tmp_vectors;
 public:
-  MultiValue( const std::size_t& nvals, const std::size_t& nder, const std::size_t& nmat=0, const std::size_t& maxcol=0, const std::size_t& nbook=0 );
-  void resize( const std::size_t& nvals, const std::size_t& nder, const std::size_t& nmat=0, const std::size_t& maxcol=0, const std::size_t& nbook=0 );
+  MultiValue( const std::size_t& nvals, const std::size_t& nder, const std::size_t& nmat=0, const std::size_t& maxcol=0 );
+  void resize( const std::size_t& nvals, const std::size_t& nder, const std::size_t& nmat=0, const std::size_t& maxcol=0 );
 /// Set the task index prior to the loop
   void setTaskIndex( const std::size_t& tindex );
 ///
@@ -136,10 +134,6 @@ class MultiValue {
   unsigned getNumberActive( const std::size_t& ival ) const ;
 ///
   unsigned getActiveIndex( const unsigned& ) const ;
-/// Get the matrix bookeeping array
-  const std::vector<unsigned> & getMatrixBookeeping() const ;
-  void stashMatrixElement( const unsigned& nmat, const unsigned& rowstart, const unsigned& jcol, const double& val );
-  double getStashedMatrixElement( const unsigned& nmat, const unsigned& jcol ) const ;
 /// Get the bookeeping stuff for the derivatives wrt to rows of matrix
   void setNumberOfMatrixRowDerivatives( const unsigned& nmat, const unsigned& nind );
   unsigned getNumberOfMatrixRowDerivatives( const unsigned& nmat ) const ;
@@ -303,23 +297,6 @@ std::vector<Tensor>& MultiValue::getFirstAtomVirialVector() {
   return tmp_atom_virial;
 }
 
-inline
-void MultiValue::stashMatrixElement( const unsigned& nmat, const unsigned& rowstart, const unsigned& jcol, const double& val ) {
-  plumed_dbg_assert( jcol<nmatrix_cols && rowstart + matrix_bookeeping[rowstart]<matrix_bookeeping.size() && nmatrix_cols*nmat + matrix_bookeeping[rowstart]<matrix_row_stash.size() );
-  matrix_bookeeping[rowstart]++; matrix_bookeeping[rowstart + matrix_bookeeping[rowstart]]=jcol; matrix_row_stash[ nmatrix_cols*nmat + jcol] = val;
-}
-
-inline
-double MultiValue::getStashedMatrixElement( const unsigned& nmat, const unsigned& jcol ) const {
-  plumed_dbg_assert( nmatrix_cols*nmat + jcol<matrix_row_stash.size() );
-  return matrix_row_stash[ nmatrix_cols*nmat + jcol ];
-}
-
-inline
-const std::vector<unsigned> & MultiValue::getMatrixBookeeping() const {
-  return matrix_bookeeping;
-}
-
 inline
 void MultiValue::setNumberOfMatrixRowDerivatives( const unsigned& nmat, const unsigned& nind ) {
   plumed_dbg_assert( nmat<matrix_row_nderivatives.size() && nind<=matrix_row_derivative_indices[nmat].size() );