diff --git a/src/crystallization/DFSClustering.cpp b/src/crystallization/DFSClustering.cpp
index 26f5e7aa105d7f91c6f6dc818e0e79efd5acf5c0..174ebb4433fd01575ff1dc4b60ab987d4afe81b7 100644
--- a/src/crystallization/DFSClustering.cpp
+++ b/src/crystallization/DFSClustering.cpp
@@ -150,7 +150,7 @@ void DFSClustering::completeCalculation(){
    ActionWithVessel::doJobsRequiredBeforeTaskList();  // Note we loose adjacency data by doing this
    // Get size for buffer
    unsigned bsize=0; std::vector<double> buffer( getSizeOfBuffer( bsize ), 0.0 );
-   std::vector<double> vals( getNumberOfQuantities() );
+   std::vector<double> vals( getNumberOfQuantities() ); std::vector<unsigned> der_index;
    vesselbase::MultiValue myvals( getNumberOfQuantities(), getNumberOfDerivatives() );
    vesselbase::MultiValue bvals( getNumberOfQuantities(), getNumberOfDerivatives() );
    // Get rid of bogus derivatives
@@ -163,7 +163,7 @@ void DFSClustering::completeCalculation(){
        for(unsigned i=0;i<vals.size();++i) myvals.setValue( i, vals[i] );
        if( !doNotCalculateDerivatives() ) getVectorDerivatives( i, false, myvals );
        // Run calculate all vessels
-       calculateAllVessels( i, myvals, bvals, buffer );
+       calculateAllVessels( i, myvals, bvals, buffer, der_index );
        myvals.clearAll();
    }
    // MPI Gather everything
diff --git a/src/crystallization/GradientVessel.cpp b/src/crystallization/GradientVessel.cpp
index 0f2d6372dbfebd5719bac50621cdfc11e52120f7..fbbb2adbced0d282228130a6fb31d712d52d6c2d 100644
--- a/src/crystallization/GradientVessel.cpp
+++ b/src/crystallization/GradientVessel.cpp
@@ -34,15 +34,13 @@ private:
   bool isdens;
   unsigned nweights, ncomponents;
   std::vector<unsigned> starts;
-//  std::vector<double> val_interm;
-//  Matrix<double> der_interm;
 public:
   static void registerKeywords( Keywords& keys );
   static void reserveKeyword( Keywords& keys );
   GradientVessel( const vesselbase::VesselOptions& da );
   std::string function_description();
   void resize();
-  bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer );
+  bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
   void finish( const std::vector<double>& buffer );
 };
 
@@ -68,7 +66,6 @@ FunctionVessel(da)
    } else {
        ncomponents = 1;
    }
-   // ncomponents = vg->vend;
 
    starts.push_back(0);
    if( vg->nbins[0]>0 ){
@@ -96,34 +93,26 @@ void GradientVessel::resize(){
      unsigned nder=getAction()->getNumberOfDerivatives();
      resizeBuffer( (1+nder)*(ncomponents+1)*nweights );
      setNumberOfDerivatives( nder );
-//     val_interm.resize( ncomponents*nweights );
-//     der_interm.resize( ncomponents*nweights, nder );
   } else {
      setNumberOfDerivatives(0); 
      resizeBuffer( (ncomponents+1)*nweights );
-//     val_interm.resize( ncomponents*nweights );
   }
 }
 
-bool GradientVessel::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer ){
+bool GradientVessel::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
   unsigned nder=getAction()->getNumberOfDerivatives();
   unsigned wstart, cstart; if( ncomponents==1 ){ cstart=1; wstart=2; } else { cstart=2; wstart=2+ncomponents; }
 
   for(unsigned iw=0;iw<nweights;++iw){
       unsigned xx = (ncomponents+1)*iw;
-      double weight=myvals.get(wstart+iw); // getAction()->getElementValue(ncomponents + iw);
+      double weight=myvals.get(wstart+iw); 
       buffer[bufstart+xx*(nder+1)] += weight;
       myvals.chainRule( wstart + iw, xx, 1, 0, 1.0, bufstart, buffer );
-      // addValueIgnoringTolerance( xx, weight ); 
-      // getAction()->chainRuleForElementDerivatives( xx , ncomponents + iw, 1.0, bufstart, buffer );
       for(unsigned jc=0;jc<ncomponents;++jc){
-          double colvar=myvals.get( cstart + jc );   // getAction()->getElementValue( jc );
-          // addValueIgnoringTolerance( xx + 1 + jc, weight*colvar );
+          double colvar=myvals.get( cstart + jc );   
           buffer[bufstart+(xx+1+jc)*(nder+1) ] += weight*colvar;
           myvals.chainRule( cstart + jc, xx + 1 + jc, 1, 0, weight, bufstart, buffer );
           myvals.chainRule( wstart + iw, xx + 1 + jc, 1, 0, colvar, bufstart, buffer );
-          // getAction()->chainRuleForElementDerivatives( xx + 1 + jc, jc, weight, bufstart, buffer );
-          // getAction()->chainRuleForElementDerivatives( xx + 1 + jc, ncomponents + iw, colvar, bufstart, buffer );   
       }
   }
 
@@ -137,25 +126,24 @@ void GradientVessel::finish( const std::vector<double>& buffer ){
 
   if( isdens ){
       for(unsigned iw=0;iw<nweights;++iw){
-          val_interm[iw] = buffer[bufstart + 2*iw*(1+nder)]; // getFinalValue( 2*iw );
+          val_interm[iw] = buffer[bufstart + 2*iw*(1+nder)]; 
           if( getAction()->derivativesAreRequired() ){
               unsigned wstart = bufstart + 2*iw*(nder+1) + 1;
-              for(unsigned jder=0;jder<nder;++jder) der_interm( iw, jder ) += buffer[ wstart + jder ]; // getBufferElement( wstart + jder );
+              for(unsigned jder=0;jder<nder;++jder) der_interm( iw, jder ) += buffer[ wstart + jder ]; 
           }
       }
   } else {
       for(unsigned iw=0;iw<nweights;++iw){
           unsigned xx = (ncomponents+1)*iw;
-          double ww=buffer[bufstart + xx*(1+nder)];  // getFinalValue( xx );
-          for(unsigned jc=0;jc<ncomponents;++jc) val_interm[ iw*ncomponents + jc ] = buffer[bufstart + (xx+1+jc)*(1+nder)] / ww; //getFinalValue( xx + 1 + jc ) / ww;
+          double ww=buffer[bufstart + xx*(1+nder)];  
+          for(unsigned jc=0;jc<ncomponents;++jc) val_interm[ iw*ncomponents + jc ] = buffer[bufstart + (xx+1+jc)*(1+nder)] / ww; 
           if( getAction()->derivativesAreRequired() ){
               unsigned wstart = bufstart + xx*(nder+1) + 1;
               for(unsigned jc=0;jc<ncomponents;++jc){
                   unsigned bstart = bufstart + ( xx + 1 + jc )*(nder+1) + 1;
-                  double val = buffer[bufstart + (nder+1)*(xx+1+jc)]; // getFinalValue( xx + 1 + jc );
+                  double val = buffer[bufstart + (nder+1)*(xx+1+jc)]; 
                   for(unsigned jder=0;jder<nder;++jder) 
                      der_interm( iw*ncomponents + jc, jder ) = (1.0/ww)*buffer[bstart + jder] - (val/(ww*ww))*buffer[wstart + jder]; 
-                     // (1.0/ww)*getBufferElement( bstart + jder ) - (val/(ww*ww))*getBufferElement( wstart + jder );
               }
           }
       }
diff --git a/src/crystallization/VectorMean.cpp b/src/crystallization/VectorMean.cpp
index 64c1927014970844564478f04cd044b22123d5e5..965ab55eb9195b73608060f2974225c016946e44 100644
--- a/src/crystallization/VectorMean.cpp
+++ b/src/crystallization/VectorMean.cpp
@@ -35,7 +35,7 @@ public:
   VectorMean( const vesselbase::VesselOptions& da );
   std::string function_description();
   void resize();
-  bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer );
+  bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
   void finish( const std::vector<double>& buffer );
 };
 
@@ -73,7 +73,7 @@ void VectorMean::resize(){
   }
 }
 
-bool VectorMean::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer ){
+bool VectorMean::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
   unsigned ncomp=getAction()->getNumberOfQuantities()-2, nder=getAction()->getNumberOfDerivatives();
 
   double weight=myvals.get(0); plumed_dbg_assert( weight>=getTolerance() ); 
diff --git a/src/crystallization/VectorSum.cpp b/src/crystallization/VectorSum.cpp
index f12e098f988bd470e7d47e4ee6ae10c85bd01ad2..1c6b3bf7d8c5727213994653dfb2329baa9c451b 100644
--- a/src/crystallization/VectorSum.cpp
+++ b/src/crystallization/VectorSum.cpp
@@ -35,7 +35,7 @@ public:
   VectorSum( const vesselbase::VesselOptions& da );
   std::string function_description();
   void resize();
-  bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer );
+  bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
   void finish( const std::vector<double>& buffer );
 };
 
@@ -73,7 +73,7 @@ void VectorSum::resize(){
   }
 }
 
-bool VectorSum::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer ){
+bool VectorSum::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
   unsigned ncomp=getAction()->getNumberOfQuantities()-2, nder=getAction()->getNumberOfDerivatives();
 
   double weight=myvals.get(0); 
diff --git a/src/mapping/SpathVessel.cpp b/src/mapping/SpathVessel.cpp
index 525244eb2fc36b68638ce0225554f98fe9d2a5d9..2ca3ad40c908c680c5e8dcf95a44f568d695f24a 100644
--- a/src/mapping/SpathVessel.cpp
+++ b/src/mapping/SpathVessel.cpp
@@ -38,7 +38,7 @@ public:
   SpathVessel( const vesselbase::VesselOptions& da );
   std::string function_description();
   void prepare();
-  bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer );
+  bool calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const ;
 };
 
 PLUMED_REGISTER_VESSEL(SpathVessel,"SPATH")
@@ -74,7 +74,7 @@ void SpathVessel::prepare(){
   foundoneclose=false;
 }
 
-bool SpathVessel::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer ){
+bool SpathVessel::calculate( const unsigned& current, vesselbase::MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const {
   double pp=mymap->getPropertyValue( current, mycoordnumber ), weight=myvals.get(0);
   if( weight<getTolerance() ) return false;
   buffer[bufstart] += weight*pp; buffer[bufstart+1+nderiv] += weight; 
diff --git a/src/mapping/ZpathVessel.cpp b/src/mapping/ZpathVessel.cpp
index c5ee11d54c8121a674a14ae5501cec3ec6c2de6b..b22e89f1a1c6dd99725f7fb26221a4c898a52574 100644
--- a/src/mapping/ZpathVessel.cpp
+++ b/src/mapping/ZpathVessel.cpp
@@ -34,7 +34,7 @@ public:
   static void reserveKeyword( Keywords& keys );
   ZpathVessel( const vesselbase::VesselOptions& da );
   std::string function_description();
-  double calcTransform( const double& val, double& dv );
+  double calcTransform( const double& val, double& dv ) const ;
   double finalTransform( const double& val, double& dv );
 };
 
@@ -61,7 +61,7 @@ std::string ZpathVessel::function_description(){
   return "the distance from the low-dimensional manifold";
 }
 
-double ZpathVessel::calcTransform( const double& val, double& dv ){
+double ZpathVessel::calcTransform( const double& val, double& dv ) const {
   dv=0.0; return 1.0;
 }
 
diff --git a/src/multicolvar/BridgedMultiColvarFunction.cpp b/src/multicolvar/BridgedMultiColvarFunction.cpp
index 1cdb192ef2183b9d40f3b9b953d417f77fff415e..6b82ce32298ac15bf99387eef3e65c8f9095b699 100644
--- a/src/multicolvar/BridgedMultiColvarFunction.cpp
+++ b/src/multicolvar/BridgedMultiColvarFunction.cpp
@@ -51,10 +51,6 @@ MultiColvarBase(ao)
   for(unsigned i=0;i<mycolv->getFullNumberOfTasks();++i) addTaskToList( mycolv->getTaskCode(i) );
 }
 
-//void BridgedMultiColvarFunction::getIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& indices ){
-//  mycolv->getIndexList( ntotal, jstore, maxder, indices );
-//}
-
 void BridgedMultiColvarFunction::transformBridgedDerivatives( const unsigned& current, vesselbase::MultiValue& invals, vesselbase::MultiValue& outvals ){
   completeTask( current, invals, outvals );
   
@@ -70,44 +66,12 @@ void BridgedMultiColvarFunction::transformBridgedDerivatives( const unsigned& cu
   outvals.sortActiveList(); 
 }
 
-/// These two functions (commented out one underneath) are needed and are a bit tricky
 void BridgedMultiColvarFunction::performTask( const unsigned& taskIndex, const unsigned& current, vesselbase::MultiValue& myvals ){
-  //atoms_with_derivatives.deactivateAll();
-
-  // GAT this is for recomputing - this needs to work -- how?
-  if( !myBridgeVessel->prerequisitsCalculated() ){
-      mycolv->performTask( taskIndex, current, myvals );
-  } 
-
-  // completeTask( myvals, outvals );  -- also need this
-  //atoms_with_derivatives.updateActiveMembers();
+  vesselbase::MultiValue invals( mycolv->getNumberOfQuantities(), mycolv->getNumberOfDerivatives() );
+  mycolv->performTask( taskIndex, current, invals );
+  transformBridgedDerivatives( taskIndex, invals, myvals ); 
 }
 
-// Vector BridgedMultiColvarFunction::retrieveCentralAtomPos(){
-//  if( atomsWithCatomDer.getNumberActive()==0 ){
-//      Vector cvec = mycolv->retrieveCentralAtomPos();
-//
-//      // Copy the value and derivatives from the MultiColvar
-//      atomsWithCatomDer.emptyActiveMembers();
-//      for(unsigned i=0;i<3;++i){
-//         setElementValue( getCentralAtomElementIndex() + i, mycolv->getElementValue( mycolv->getCentralAtomElementIndex() + i ) );
-//         unsigned nbase = ( getCentralAtomElementIndex() + i)*getNumberOfDerivatives();
-//         unsigned nbas2 = ( mycolv->getCentralAtomElementIndex() + i )*mycolv->getNumberOfDerivatives();
-//         for(unsigned j=0;j<mycolv->atomsWithCatomDer.getNumberActive();++j){
-//             unsigned n=mycolv->atomsWithCatomDer[j], nx=3*n; atomsWithCatomDer.activate(n);
-//             addElementDerivative(nbase + nx + 0, mycolv->getElementDerivative(nbas2 + nx + 0) );
-//             addElementDerivative(nbase + nx + 1, mycolv->getElementDerivative(nbas2 + nx + 1) );
-//             addElementDerivative(nbase + nx + 2, mycolv->getElementDerivative(nbas2 + nx + 2) ); 
-//         } 
-//      }
-//      atomsWithCatomDer.updateActiveMembers();  // This can perhaps be faster
-//      return cvec;
-//  }
-//  Vector cvec;
-//  for(unsigned i=0;i<3;++i) cvec[i]=getElementValue(1+i);
-//  return cvec;
-// }
-
 void BridgedMultiColvarFunction::calculateNumericalDerivatives( ActionWithValue* a ){
   if(!a){
     a=dynamic_cast<ActionWithValue*>(this);
diff --git a/src/vesselbase/ActionWithVessel.cpp b/src/vesselbase/ActionWithVessel.cpp
index 47529851960ac982eb8a39d9b39632ff2d9e6760..9f1639406cdddd7946c4e5d1c5c556b0a47d1951 100644
--- a/src/vesselbase/ActionWithVessel.cpp
+++ b/src/vesselbase/ActionWithVessel.cpp
@@ -56,6 +56,7 @@ ActionWithVessel::ActionWithVessel(const ActionOptions&ao):
   lowmem(false),
   noderiv(true),
   actionIsBridged(false),
+  mydata(NULL),
   contributorsAreUnlocked(false),
   weightHasDerivatives(false)
 {
@@ -104,6 +105,10 @@ void ActionWithVessel::addVessel( Vessel* vv ){
   ShortcutVessel* sv=dynamic_cast<ShortcutVessel*>(vv);
   if(!sv){ vv->checkRead(); functions.push_back(vv); }
   else { delete sv; }
+
+  StoreDataVessel* mm=dynamic_cast<StoreDataVessel*>( vv );
+  if( mydata && mm ) error("cannot have more than one StoreDataVessel in one action");
+  else mydata=mm;
 }
 
 BridgeVessel* ActionWithVessel::addBridgingVessel( ActionWithVessel* tome ){
@@ -117,15 +122,12 @@ BridgeVessel* ActionWithVessel::addBridgingVessel( ActionWithVessel* tome ){
 }
 
 StoreDataVessel* ActionWithVessel::buildDataStashes( const bool& allow_wcutoff, const double& wtol ){
-  for(unsigned i=0;i<functions.size();++i){
-      StoreDataVessel* vsv=dynamic_cast<StoreDataVessel*>( functions[i] );
-      if( vsv ) return vsv;
-  }
+  if(mydata) return mydata;
   
   VesselOptions da("","",0,"",this);
-  StoreDataVessel* mydata=new StoreDataVessel(da);
-  if( allow_wcutoff ) mydata->setHardCutoffOnWeight( wtol );
-  addVessel(mydata);
+  StoreDataVessel* mm=new StoreDataVessel(da);
+  if( allow_wcutoff ) mm->setHardCutoffOnWeight( wtol );
+  addVessel(mm);
 
   // Make sure resizing of vessels is done
   resizeFunctions();
@@ -278,6 +280,9 @@ void ActionWithVessel::runAllTasks(){
   // Get size for buffer
   unsigned bsize=0, bufsize=getSizeOfBuffer( bsize ); 
 
+  std::vector<unsigned> der_list;
+  if( mydata ) der_list.resize( mydata->getSizeOfDerivativeList(), 0 ); 
+
   // Build storage stuff for loop
   std::vector<double> buffer( bufsize, 0.0 );
   MultiValue myvals( getNumberOfQuantities(), getNumberOfDerivatives() );
@@ -304,11 +309,17 @@ void ActionWithVessel::runAllTasks(){
       // Now calculate all the functions
       // If the contribution of this quantity is very small at neighbour list time ignore it
       // untill next neighbour list time
-      if( !calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, buffer ) && contributorsAreUnlocked ) deactivate_task( indexOfTaskInFullList[i] );
+      if( !calculateAllVessels( indexOfTaskInFullList[i], myvals, bvals, buffer, der_list ) && contributorsAreUnlocked ) deactivate_task( indexOfTaskInFullList[i] );
+
+      // Clear the value
       myvals.clearAll();
   }
   // MPI Gather everything
   if( !serial && buffer.size()>0 ) comm.Sum( buffer );
+  // MPI Gather index stores
+  if( mydata && !lowmem && !noderiv ){ 
+     comm.Sum( der_list ); mydata->setActiveValsAndDerivatives( der_list ); 
+  }
   // Update the elements that are makign contributions to the sum here
   // this causes problems if we do it in prepare
   if( !serial && contributorsAreUnlocked ) comm.Sum( taskFlags );
@@ -320,20 +331,12 @@ void ActionWithVessel::transformBridgedDerivatives( const unsigned& current, Mul
   plumed_error();
 }
 
-// void ActionWithVessel::getIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& indices ){
-//   indices[jstore]=getNumberOfDerivatives();
-//   if( indices[jstore]>maxder ) error("too many derivatives to store. Run with LOWMEM");
-// 
-//   unsigned kder = ntotal + jstore*getNumberOfDerivatives();
-//   for(unsigned jder=0;jder<getNumberOfDerivatives();++jder){ indices[ kder ] = jder; kder++; }
-// }
-
-bool ActionWithVessel::calculateAllVessels( const unsigned& taskCode, MultiValue& myvals, MultiValue& bvals, std::vector<double>& buffer ){
+bool ActionWithVessel::calculateAllVessels( const unsigned& taskCode, MultiValue& myvals, MultiValue& bvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ){
   bool keep=false; 
   for(unsigned j=0;j<functions.size();++j){
       // Calculate returns a bool that tells us if this particular
       // quantity is contributing more than the tolerance
-      if( functions[j]->calculate( taskCode, functions[j]->transformDerivatives(taskCode, myvals, bvals), buffer ) ) keep=true;
+      if( functions[j]->calculate( taskCode, functions[j]->transformDerivatives(taskCode, myvals, bvals), buffer, der_list ) ) keep=true;
       if( !actionIsBridged ) bvals.clearAll(); 
   }
   return keep;
diff --git a/src/vesselbase/ActionWithVessel.h b/src/vesselbase/ActionWithVessel.h
index e57a93ce18aec3e4ac413c30e9cb412c5036292e..f8ea4c7987d243a1f6cc2f448470a74deda3f2a2 100644
--- a/src/vesselbase/ActionWithVessel.h
+++ b/src/vesselbase/ActionWithVessel.h
@@ -68,6 +68,8 @@ private:
   double nl_tolerance;
 /// Pointers to the functions we are using on each value
   std::vector<Vessel*> functions;
+/// A pointer to the object that stores data
+  StoreDataVessel* mydata;
 /// Tempory storage for forces
   std::vector<double> tmpforces;
 /// Ths full list of tasks we have to perform
@@ -110,7 +112,7 @@ protected:
   void resizeFunctions();
 /// This loops over all the vessels calculating them and also 
 /// sets all the element derivatives equal to zero
-  bool calculateAllVessels( const unsigned& taskCode, MultiValue& myvals, MultiValue& bvals, std::vector<double>& buffer );
+  bool calculateAllVessels( const unsigned& taskCode, MultiValue& myvals, MultiValue& bvals, std::vector<double>& buffer, std::vector<unsigned>& der_list );
 /// Retrieve the forces from all the vessels (used in apply)
   bool getForcesFromVessels( std::vector<double>& forcesToApply );
 /// Is the calculation being done in serial
diff --git a/src/vesselbase/Between.cpp b/src/vesselbase/Between.cpp
index 2a3dabea53d84ca2f8aed576de6b59c915c36045..af2024b0926eb0f00091bf60077c45042614bf79 100644
--- a/src/vesselbase/Between.cpp
+++ b/src/vesselbase/Between.cpp
@@ -67,7 +67,7 @@ std::string Between::function_description(){
   return "the number of values " + hist.description();
 }
 
-double Between::calcTransform( const double& val, double& dv ){
+double Between::calcTransform( const double& val, double& dv ) const {
   double f = hist.calculate(val, dv); return f; 
 }
 
diff --git a/src/vesselbase/Between.h b/src/vesselbase/Between.h
index 593f1993d17378f5331fde6bc429b0f15eb32560..02f6eb059bc85a1419b37b4a394791107c5eefd3 100644
--- a/src/vesselbase/Between.h
+++ b/src/vesselbase/Between.h
@@ -36,7 +36,7 @@ public:
   static void reserveKeyword( Keywords& keys );
   Between( const VesselOptions& da );
   std::string function_description();
-  double calcTransform( const double& val, double& dv );
+  double calcTransform( const double& val, double& dv ) const ;
   double getCutoff();
 };
 
diff --git a/src/vesselbase/BridgeVessel.cpp b/src/vesselbase/BridgeVessel.cpp
index dceec770437191c5a221de5ccb4349c59eaf7fcb..41cb4eb84c86347b264cb391bd7c5f6fc1388b1d 100644
--- a/src/vesselbase/BridgeVessel.cpp
+++ b/src/vesselbase/BridgeVessel.cpp
@@ -28,8 +28,8 @@ namespace vesselbase {
 
 BridgeVessel::BridgeVessel( const VesselOptions& da ):
 Vessel(da),
-inum(0),
-in_normal_calculate(false)
+inum(0)
+// in_normal_calculate(false)
 {
 }
 
@@ -71,13 +71,13 @@ MultiValue& BridgeVessel::transformDerivatives( const unsigned& current, MultiVa
   return outvals;
 }
 
-bool BridgeVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer ){
-  in_normal_calculate=true;
+bool BridgeVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
+  // in_normal_calculate=true;
   if( myvals.get(0)<myOutputAction->getTolerance() ){
       return ( !myOutputAction->contributorsAreUnlocked || myvals.get(0)>=myOutputAction->getNLTolerance() );
   }
-  bool keep=myOutputAction->calculateAllVessels( current, myvals, myvals, buffer );    
-  in_normal_calculate=false;
+  bool keep=myOutputAction->calculateAllVessels( current, myvals, myvals, buffer, der_list );    
+  // in_normal_calculate=false;
   return ( !myOutputAction->contributorsAreUnlocked || keep );
 }
 
diff --git a/src/vesselbase/BridgeVessel.h b/src/vesselbase/BridgeVessel.h
index fc6c914bd3dc06a01a70b4029c504552c393d768..fedaf39ea922c5f8236633baf602450bacc29e22 100644
--- a/src/vesselbase/BridgeVessel.h
+++ b/src/vesselbase/BridgeVessel.h
@@ -40,7 +40,7 @@ it is created in a different Action however.  At the moment this is used for reg
 class BridgeVessel : public Vessel {
 private:
   unsigned inum;
-  bool in_normal_calculate;
+  // bool in_normal_calculate;
   std::vector<double> mynumerical_values;
   ActionWithVessel* myOutputAction;
   ActionWithValue* myOutputValues;
@@ -63,20 +63,13 @@ public:
 /// This transforms the derivatives using the output value
   MultiValue& transformDerivatives( const unsigned& current, MultiValue& invals, MultiValue& outvals );
 /// Actually do the calculation
-  bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer );
+  bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const ;
 /// Finish the calculation
   void finish( const std::vector<double>& buffer );
 /// Calculate numerical derivatives
   void completeNumericalDerivatives();
-/// This is used to tell if the bridge has been called in recompute
-  bool prerequisitsCalculated();
 };
 
-inline
-bool BridgeVessel::prerequisitsCalculated(){
-  return in_normal_calculate;
-}
-
 }
 }
 #endif
diff --git a/src/vesselbase/FunctionVessel.cpp b/src/vesselbase/FunctionVessel.cpp
index 5aad279d770862235bd21d97e04e6c03faac9acf..02e7a6d695258e9851a4f56f79f2c87003dd29ec 100644
--- a/src/vesselbase/FunctionVessel.cpp
+++ b/src/vesselbase/FunctionVessel.cpp
@@ -73,7 +73,7 @@ void FunctionVessel::setNumberOfDerivatives( const unsigned& nder ){
   final_value->resizeDerivatives( nder );
 }
 
-bool FunctionVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer ){
+bool FunctionVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
   double weight=myvals.get(0); 
   plumed_dbg_assert( weight>=getTolerance() );  
 
@@ -96,7 +96,7 @@ bool FunctionVessel::calculate( const unsigned& current, MultiValue& myvals, std
   return true;
 }
 
-double FunctionVessel::calcTransform( const double& , double& ){ 
+double FunctionVessel::calcTransform( const double& , double& ) const { 
   plumed_error(); return 1.0; 
 }
 
diff --git a/src/vesselbase/FunctionVessel.h b/src/vesselbase/FunctionVessel.h
index 08157928f0ac309db713e91edeb979bad261d34a..9c787c96a816a9919c5f090cdb406afd116cca25 100644
--- a/src/vesselbase/FunctionVessel.h
+++ b/src/vesselbase/FunctionVessel.h
@@ -53,8 +53,6 @@ protected:
   bool usetol;
 /// Set the final value
   void setOutputValue( const double& val );
-/// This does a combination of the product and chain rules
-  void mergeFinalDerivatives( const std::vector<double>& df );
 /// Resize the vector containing the derivatives
   void setNumberOfDerivatives( const unsigned& nder );
 /// Return a pointer to the final value
@@ -71,9 +69,9 @@ public:
 /// The rest of the description of what we are calculating
   virtual std::string function_description()=0;
 /// Do the calcualtion
-  virtual bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer );
+  virtual bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const ;
 /// Do any transformations of the value that are required
-  virtual double calcTransform( const double& val, double& df );
+  virtual double calcTransform( const double& val, double& df ) const ;
 /// Finish the calculation of the quantity
   virtual void finish( const std::vector<double>& buffer );
 /// Finish with any transforms required
diff --git a/src/vesselbase/LessThan.cpp b/src/vesselbase/LessThan.cpp
index 479f3afac78e3c101afd4942af97c5f67b09e2ff..2c3b98ce40843194a21ff817f4346250d185127f 100644
--- a/src/vesselbase/LessThan.cpp
+++ b/src/vesselbase/LessThan.cpp
@@ -54,7 +54,7 @@ std::string LessThan::function_description(){
   return "the number of values less than " + sf.description();
 }
 
-double LessThan::calcTransform( const double& val, double& dv ){
+double LessThan::calcTransform( const double& val, double& dv ) const {
   double f = sf.calculate(val, dv); dv*=val; return f;
 }
 
diff --git a/src/vesselbase/LessThan.h b/src/vesselbase/LessThan.h
index 24bf15148f2c37e5da16ecb28025f2a318467428..7570618ecc0fb61377652afb8f3b93dc2a8b5342 100644
--- a/src/vesselbase/LessThan.h
+++ b/src/vesselbase/LessThan.h
@@ -37,7 +37,7 @@ public:
   static void reserveKeyword( Keywords& keys ); 
   LessThan( const VesselOptions& da );
   std::string function_description();
-  double calcTransform( const double& val, double& dv );
+  double calcTransform( const double& val, double& dv ) const ;
   double getCutoff();
 };
 
diff --git a/src/vesselbase/Max.cpp b/src/vesselbase/Max.cpp
index 13b5a995ce6770dffcadaad33e84569f225a283c..e27dabd79531c6ce712a3410d52e60ff5b67e252 100644
--- a/src/vesselbase/Max.cpp
+++ b/src/vesselbase/Max.cpp
@@ -34,7 +34,7 @@ public:
   static void reserveKeyword( Keywords& keys );
   Max( const VesselOptions& da );
   std::string function_description();
-  double calcTransform( const double& val, double& dv );
+  double calcTransform( const double& val, double& dv ) const ;
   double finalTransform( const double& val, double& dv );
 };
 
@@ -69,7 +69,7 @@ std::string Max::function_description(){
   return "the maximum value. Beta is equal to " + str_beta;
 }
 
-double Max::calcTransform( const double& val, double& dv ){
+double Max::calcTransform( const double& val, double& dv ) const {
   double f = exp(val/beta); dv=f/beta; return f;
 }
 
diff --git a/src/vesselbase/Mean.cpp b/src/vesselbase/Mean.cpp
index 529788360d4ecf1a270300ed77e00cd42d6a48de..764629e7cc916e990e3091c20da5314f1cceee43 100644
--- a/src/vesselbase/Mean.cpp
+++ b/src/vesselbase/Mean.cpp
@@ -32,7 +32,7 @@ public:
   static void reserveKeyword( Keywords& keys );
   Mean( const vesselbase::VesselOptions& da );
   std::string function_description();
-  double calcTransform( const double& val, double& dv );
+  double calcTransform( const double& val, double& dv ) const ;
 };
 
 PLUMED_REGISTER_VESSEL(Mean,"MEAN")
@@ -58,7 +58,7 @@ std::string Mean::function_description(){
   return "the mean value";
 }
 
-double Mean::calcTransform( const double& val, double& dv ){
+double Mean::calcTransform( const double& val, double& dv ) const {
   dv=1.0; return val;
 }
 
diff --git a/src/vesselbase/Min.cpp b/src/vesselbase/Min.cpp
index a14c1f6bb8fc51c021100a3be15ed926105bcac7..8e5ba1a422d62f32b585d4471851598d044e040c 100644
--- a/src/vesselbase/Min.cpp
+++ b/src/vesselbase/Min.cpp
@@ -34,7 +34,7 @@ public:
   static void reserveKeyword( Keywords& keys );
   Min( const VesselOptions& da );
   std::string function_description();
-  double calcTransform( const double& val, double& dv );
+  double calcTransform( const double& val, double& dv ) const ;
   double finalTransform( const double& val, double& dv );
 };
 
@@ -68,7 +68,7 @@ std::string Min::function_description(){
   return "the minimum value. Beta is equal to " + str_beta;
 }
 
-double Min::calcTransform( const double& val, double& dv ){
+double Min::calcTransform( const double& val, double& dv ) const {
   double f = exp(beta/val); dv=f/(val*val);
   return f; 
 }
diff --git a/src/vesselbase/MoreThan.cpp b/src/vesselbase/MoreThan.cpp
index 8012481550d10fc9bfd80a11f75385ceeb151c37..938da5d6b688f59ab604bf2258796f9eeee7f0f2 100644
--- a/src/vesselbase/MoreThan.cpp
+++ b/src/vesselbase/MoreThan.cpp
@@ -36,7 +36,7 @@ public:
   static void reserveKeyword( Keywords& keys );
   MoreThan( const VesselOptions& da );
   std::string function_description();
-  double calcTransform( const double& val, double& dv );
+  double calcTransform( const double& val, double& dv ) const ;
 };
 
 PLUMED_REGISTER_VESSEL(MoreThan,"MORE_THAN")
@@ -68,7 +68,7 @@ std::string MoreThan::function_description(){
   return "the number of values more than " + sf.description();
 }
 
-double MoreThan::calcTransform( const double& val, double& dv ){
+double MoreThan::calcTransform( const double& val, double& dv ) const {
   double f = 1.0 - sf.calculate(val, dv); dv*=-val; return f; 
 }
 
diff --git a/src/vesselbase/ShortcutVessel.h b/src/vesselbase/ShortcutVessel.h
index b1a2d88f2fd8d274a9aa8fb5157120d53f1775ce..0e58d50c5546457b1b1450d86800811b78c1a83f 100644
--- a/src/vesselbase/ShortcutVessel.h
+++ b/src/vesselbase/ShortcutVessel.h
@@ -39,7 +39,7 @@ public:
   ShortcutVessel( const VesselOptions& );
   std::string description(){ return ""; }
   void resize(){ plumed_error(); }
-  bool calculate( const unsigned& taskCode, MultiValue& myvals, std::vector<double>& buffer ){ plumed_error(); }
+  bool calculate( const unsigned& taskCode, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const { plumed_error(); }
   void finish( const std::vector<double>& buffer ){ plumed_error(); }
   bool applyForce( std::vector<double>& forces ){ plumed_error(); }
 };
diff --git a/src/vesselbase/StoreDataVessel.cpp b/src/vesselbase/StoreDataVessel.cpp
index b1b35f484c82eccc8a58083873f9a84b035a7406..b496ce841049ab852eaf2815a52ecb4019bce1db 100644
--- a/src/vesselbase/StoreDataVessel.cpp
+++ b/src/vesselbase/StoreDataVessel.cpp
@@ -49,80 +49,46 @@ bool StoreDataVessel::weightCutoffIsOn() const {
   return hard_cut;
 }
 
-// void StoreDataVessel::completeSetup( const unsigned& dstart, const unsigned& nvec ){
-//   if( (dstart+nvec)>getAction()->getNumberOfQuantities() ) error("this vessel can not be used with this action");
-//   data_start=dstart; vecsize = nvec; // fvec.resize( vecsize ); 
-// }
-
 void StoreDataVessel::resize(){
   plumed_dbg_assert( vecsize>0 );
-//  if( getAction()->derivativesAreRequired() ) final_derivatives.resize( getAction()->getNumberOfDerivatives() );
 
   if( getAction()->lowmem || !getAction()->derivativesAreRequired() ){
      nspace = 1;
      active_der.resize( max_lowmem_stash * ( 1 + getAction()->getNumberOfDerivatives() ) );
-//     local_derivatives.resize( max_lowmem_stash * vecsize * getAction()->getNumberOfDerivatives() );
   } else {
      nspace = 1 + getAction()->maxderivatives;
      active_der.resize( getAction()->getFullNumberOfTasks() * ( 1 + getAction()->maxderivatives ) );
   }
-  active_val.resize( getAction()->getFullNumberOfTasks() );
+  //active_val.resize( getAction()->getFullNumberOfTasks() );
   resizeBuffer( getAction()->getFullNumberOfTasks()*vecsize*nspace );
   local_buffer.resize( getAction()->getFullNumberOfTasks()*vecsize*nspace );
 }
 
-void StoreDataVessel::prepare(){
-// Clear bookeeping arrays  
-  active_der.assign(active_der.size(),0);
-  active_val.assign(active_val.size(),0);
-}
-
-// void StoreDataVessel::getIndexList( const unsigned& ntotal, const unsigned& jstore, const unsigned& maxder, std::vector<unsigned>& indices ){
-// //  getAction()->getIndexList( ntotal, jstore, maxder, indices );
-// }
-// 
-// void StoreDataVessel::setTaskToRecompute( const unsigned& ivec ){
-// // getAction()->current = getAction()->fullTaskList[ivec];
-// // getAction()->task_index = ivec;
-// }
-// 
-// void StoreDataVessel::recompute( const unsigned& ivec, const unsigned& jstore ){
-//   plumed_dbg_assert( getAction()->derivativesAreRequired() && getAction()->lowmem && jstore<max_lowmem_stash );
-//   // Set the task we want to reperform
-//   setTaskToRecompute( ivec );
-//   // Reperform the task
-//   performTask( jstore );
-//   // Store the derivatives
-//   storeDerivativesLowMem( jstore );
-//  
-//   // Clear up afterwards
-//   // getAction()->clearAfterTask();
+// void StoreDataVessel::prepare(){
+// // Clear bookeeping arrays  
+// //  active_der.assign(active_der.size(),0);
+//   //active_val.assign(active_val.size(),0);
 // }
 
-void StoreDataVessel::storeValues( const unsigned& myelem, MultiValue& myvals, std::vector<double>& buffer ){
-  unsigned ibuf = bufstart + myelem * vecsize * nspace; active_val[myelem]=1;
+void StoreDataVessel::storeValues( const unsigned& myelem, MultiValue& myvals, std::vector<double>& buffer ) const {
+  unsigned ibuf = bufstart + myelem * vecsize * nspace; // active_val[myelem]=1;
   for(unsigned icomp=0;icomp<vecsize;++icomp){
       buffer[ibuf] = myvals.get(icomp); ibuf+=nspace;
   }
 }
 
-void StoreDataVessel::storeDerivatives( const unsigned& myelem, MultiValue& myvals, std::vector<double>& buffer ){
+void StoreDataVessel::storeDerivatives( const unsigned& myelem, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
   plumed_dbg_assert( getAction()->derivativesAreRequired() && myelem<getAction()->getFullNumberOfTasks() );
-  active_der[myelem]=myvals.getNumberActive();
-  // getIndexList( act->getFullNumberOfTasks(), myelem, nspace-1, active_der );
+  der_list[myelem]=myvals.getNumberActive();
 
   // Store the values of the components and the derivatives if 
-  // unsigned nder =getAction()->getNumberOfDerivatives();
   for(unsigned icomp=0;icomp<vecsize;++icomp){
      unsigned ibuf = bufstart + myelem * ( vecsize*nspace ) + icomp*nspace + 1;
      unsigned kder = getAction()->getFullNumberOfTasks() + myelem * ( nspace - 1 );
      for(unsigned j=0;j<myvals.getNumberActive();++j){
         unsigned jder=myvals.getActiveIndex(j);
         buffer[ibuf] = myvals.getDerivative( icomp, jder ); 
-        active_der[kder] =  jder;
-        // act->getElementDerivative(nder*icomp + active_der[ kder ]);
-        //addToBufferElement( ibuf, act->getElementDerivative(nder*icomp + active_der[ kder ]) );
-        ibuf++; kder++;
+        der_list[kder] = jder; ibuf++; kder++;
      }
   }
 }
@@ -199,130 +165,51 @@ void StoreDataVessel::retrieveDerivatives( const unsigned& myelem, const bool& n
   } 
 }
 
+bool StoreDataVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const {
 
-// 
-// void StoreDataVessel::storeDerivativesLowMem( const unsigned& jstore ){
-//   plumed_dbg_assert( getAction()->derivativesAreRequired() );
-//   // Store the indexes that have derivatives
-//   ActionWithVessel* act=getAction();
-//   unsigned nder = act->getNumberOfDerivatives();
-//   getIndexList( max_lowmem_stash, jstore, nder, active_der );
-// 
-//   // Stash the derivatives
-//   for(unsigned icomp=data_start;icomp<data_start + vecsize;++icomp){
-//      unsigned ibuf = jstore * vecsize * nder + (icomp-data_start)*nder;
-//      unsigned kder = max_lowmem_stash + jstore*nder;
-//      for(unsigned jder=0;jder<active_der[jstore];++jder){
-//         // local_derivatives[ibuf] = act->getElementDerivative( nder*icomp + active_der[ kder ] );
-//         ibuf++; kder++;
-//      }
-//   }
-// }
-
-bool StoreDataVessel::calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer ){
-
-   if( !hard_cut ){
-      storeValues( current, myvals, buffer );
-      if( !(getAction()->lowmem) && getAction()->derivativesAreRequired() ) storeDerivatives( current, myvals, buffer );
-   } else if( myvals.get(0)>wtol ){
-      storeValues( current, myvals, buffer );
-      if( !(getAction()->lowmem) && getAction()->derivativesAreRequired() ) storeDerivatives( current, myvals, buffer );
-   } 
-
-//  unsigned myelem = getAction()->getCurrentPositionInTaskList();
-  // Normalize vector if it is required
-//  unsigned myelem; finishTask( myelem );
+  if( !hard_cut ){
+     storeValues( current, myvals, buffer );
+     if( !(getAction()->lowmem) && getAction()->derivativesAreRequired() ) storeDerivatives( current, myvals, buffer, der_list );
+  } else if( myvals.get(0)>wtol ){
+     storeValues( current, myvals, buffer );
+     if( !(getAction()->lowmem) && getAction()->derivativesAreRequired() ) storeDerivatives( current, myvals, buffer, der_list );
+  } 
 
-  // Store the values 
-//  if( !hard_cut ){
-//     storeValues( myelem, buffer );
-//     // Store the derivatives if we are not using low memory
-//     if( !(getAction()->lowmem) && getAction()->derivativesAreRequired() ) storeDerivativesHighMem( myelem, buffer );  
-//  } else if( getAction()->getElementValue(0)>wtol ){ 
-//     storeValues( myelem, buffer );
-//     // Store the derivatives if we are not using low memory
-//     if( !(getAction()->lowmem) && getAction()->derivativesAreRequired() ) storeDerivativesHighMem( myelem, buffer );
-//  }
-//
   return true;
 }
 
+// void StoreDataVessel::buildIndexStores( const unsigned& current, MultiValue& myvals, std::vector<unsigned>& val_index, std::vector<unsigned>& der_index ) const {
+//   if( !hard_cut ){
+//      val_index[current]=1;
+//      if( !(getAction()->lowmem) && getAction()->derivativesAreRequired() ){
+//          der_index[current]=myvals.getNumberActive();
+//          unsigned kder = getAction()->getFullNumberOfTasks() + current * ( nspace - 1 );
+//          for(unsigned j=0;j<myvals.getNumberActive();++j){ der_index[kder] = myvals.getActiveIndex(j); kder++; }
+//      }
+//   } else if( myvals.get(0)>wtol ){
+//      val_index[current]=1;
+//      if( !(getAction()->lowmem) && getAction()->derivativesAreRequired() ){
+//          der_index[current]=myvals.getNumberActive();
+//          unsigned kder = getAction()->getFullNumberOfTasks() + current * ( nspace - 1 );
+//          for(unsigned j=0;j<myvals.getNumberActive();++j){ der_index[kder] = myvals.getActiveIndex(j); kder++; }
+//      } 
+//   }  
+// }
+
 void StoreDataVessel::finish( const std::vector<double>& buffer ){
   // Store the buffer locally
   for(unsigned i=0;i<local_buffer.size();++i) local_buffer[i]=buffer[bufstart+i];
-
-  // Update what variables we need
-  comm.Sum( &active_val[0], active_val.size() );
-
-  // Get the active derivatives
-  if(!getAction()->lowmem && getAction()->derivativesAreRequired() ) comm.Sum( &active_der[0], active_der.size() );
 }
 
-// double StoreDataVessel::chainRule( const unsigned& ival, const unsigned& ider, const std::vector<double>& df ){
-//   plumed_dbg_assert( getAction()->derivativesAreRequired() && df.size()==vecsize );
-//   // Clear final derivatives array
-//   final_derivatives.assign( final_derivatives.size(), 0.0 );
-// 
-//   double dfout = 0.0;
-//   if(getAction()->lowmem){
-//      plumed_dbg_assert( ival<max_lowmem_stash );
-//      unsigned maxder = getAction()->getNumberOfDerivatives();
-//      unsigned ibuf=ival*(vecsize*maxder) + ider;
-//      for(unsigned icomp=0;icomp<vecsize;++icomp){
-//          dfout+=df[icomp]*local_derivatives[ibuf];
-//          ibuf+=maxder;
-//      }  
-//   } else {
-//      plumed_dbg_assert( ival<getAction()->getFullNumberOfTasks() );
-//      unsigned ibuf=ival*(vecsize*nspace) + 1 + ider;
-//      for(unsigned icomp=0;icomp<vecsize;++icomp){
-//          dfout+=df[icomp]*local_buffer[ibuf]; //getBufferElement(ibuf);
-//          ibuf+=nspace;
-//      }
-//   }
-//   return dfout;
-// }
 
-// void StoreDataVessel::chainRule( const unsigned& ival, const std::vector<double>& df ){
-//   plumed_dbg_assert( getAction()->derivativesAreRequired() && df.size()==vecsize );
-//   // Clear final derivatives array
-//   final_derivatives.assign( final_derivatives.size(), 0.0 );
-// 
-//   if(getAction()->lowmem){
-//      plumed_dbg_assert( ival<max_lowmem_stash );
-//      unsigned maxder = getAction()->getNumberOfDerivatives();
-//      for(unsigned ider=0;ider<active_der[ival];++ider){
-//         final_derivatives[ider]=0.0;
-//         unsigned ibuf=ival*(vecsize*maxder) + ider; 
-//         for(unsigned jcomp=0;jcomp<vecsize;++jcomp){
-//             final_derivatives[ider]+=df[jcomp]*local_derivatives[ibuf]; 
-//             ibuf+=maxder;
-//         }
-//      }
-//   } else {
-//      plumed_dbg_assert( ival<getAction()->getFullNumberOfTasks() );
-//      for(unsigned ider=0;ider<active_der[ival];++ider){
-//          final_derivatives[ider]=0.0;
-//          unsigned ibuf=ival*(vecsize*nspace) + 1 + ider; 
-//          for(unsigned jcomp=0;jcomp<vecsize;++jcomp){
-//              final_derivatives[ider]+=df[jcomp]*local_buffer[ibuf];
-//              //final_derivatives[ider]+=df[jcomp]*getBufferElement(ibuf);
-//              ibuf+=nspace;
-//          }
-//      }
+void StoreDataVessel::setActiveValsAndDerivatives( const std::vector<unsigned>& der_index ){
+//   for(unsigned i=0;i<val_index.size();++i){   // Need to think about this bit GAT
+//      if( ) active_val[i]    // =val_index[i];
 //   }
-// }
-
-// void StoreDataVessel::chainRule( const unsigned& ival, const std::vector<double>& df, Value* val ){
-//   plumed_dbg_assert( getAction()->derivativesAreRequired() && val->getNumberOfDerivatives()==final_derivatives.size() );
-//   chainRule( ival, df );
-// 
-//   unsigned kder;
-//   if( getAction()->lowmem ) kder = max_lowmem_stash + ival*getAction()->getNumberOfDerivatives();
-//   else kder = getAction()->getFullNumberOfTasks() + ival*(nspace-1);
-// 
-//   for(unsigned i=0;i<active_der[ival];++i) val->addDerivative( active_der[kder+i] , final_derivatives[i] );
-// }
+  if( !getAction()->lowmem && getAction()->derivativesAreRequired() ){
+     for(unsigned i=0;i<der_index.size();++i) active_der[i]=der_index[i]; 
+  }
+}
 
 }
 }
diff --git a/src/vesselbase/StoreDataVessel.h b/src/vesselbase/StoreDataVessel.h
index 02a0f7f93147938b5fda421c216003b3d8438564..77b537cc5445bbb55a1220a2c84591ac2fa8e47f 100644
--- a/src/vesselbase/StoreDataVessel.h
+++ b/src/vesselbase/StoreDataVessel.h
@@ -51,17 +51,11 @@ private:
 /// The amount of data per vector element
   unsigned nspace;
 /// The currently active values 
-  std::vector<unsigned> active_val;
+//  std::vector<unsigned> active_val;
 /// The active derivative elements
   std::vector<unsigned> active_der;
-/// This is a tempory vector that is used to store data
-//  std::vector<double> fvec;
 /// The buffer
    std::vector<double> local_buffer;
-// /// The local derivatives
-//   std::vector<double> local_derivatives;
-// /// The final derivatives
-//   std::vector<double> final_derivatives;
 protected:
 /// Apply a hard cutoff on the weight
   bool hard_cut;
@@ -77,15 +71,9 @@ protected:
 /// Return value of nspace
   unsigned getNumberOfDerivativeSpacesPerComponent() const ;
 /// Retrieve the values from the underlying ActionWithVessel
-   void storeValues( const unsigned& , MultiValue& , std::vector<double>& );
-/// Set the Task that needs redoing
-//   void setTaskToRecompute( const unsigned& ivec );
-/// Set a component of one of the vectors
-//   void setComponent( const unsigned& , const unsigned& , const double& );
-/// This is the proper chain rule for vectors
-//    double chainRule( const unsigned&, const unsigned&, const std::vector<double>& );
-//  /// Chain rule the vector and output derivatives to a value
-//    void chainRule( const unsigned& , const std::vector<double>&, Value* );
+   void storeValues( const unsigned& , MultiValue& , std::vector<double>& ) const ;
+/// This stores the data we get from the calculation
+  void storeDerivatives( const unsigned& , MultiValue& myvals, std::vector<double>&, std::vector<unsigned>& ) const ;
 /// Get the ibuf'th local derivative value
   double getLocalDerivative( const unsigned& ibuf );
 /// Set the ibuf'th local derivative value
@@ -106,40 +94,25 @@ public:
 /// Do all resizing of data
   virtual void resize();
 /// Clear certain data before start of main loop
-  virtual void prepare();
+//  virtual void prepare();
 ///
   virtual std::string description(){ return ""; }
 /// Get the number of derivatives for the ith value
   unsigned getNumberOfDerivatives( const unsigned& );
-/// Get one of the stored indexes
-//  unsigned getStoredIndex( const unsigned& , const unsigned& );
-/// Get a component of the stored vector
-//  double getComponent( const unsigned& , const unsigned& );
-/// Recalculate a vector - used in lowmem mode
-//  virtual void recompute( const unsigned& , const unsigned& );
-/// This reperforms the task in the underlying action
-//  virtual void performTask( const unsigned& );
-/// This reperforms the task
-//  virtual void finishTask( const unsigned& ){};
-/// Chain rule and store output in local array called final_derivatives
-/// with vectors this does chain rule for dot products
-//   void chainRule( const unsigned& , const std::vector<double>& );
-/// Get the ider'th final derivative value
-//  double getFinalDerivative( const unsigned& ider ) const ;
+/// Get the size of the derivative list
+  unsigned getSizeOfDerivativeList() const ;
 /// This stores the data when not using lowmem
-  bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer );
-/// This stores the data we get from the calculation
-//   void storeDerivativesLowMem( const unsigned& );
-/// This stores the data we get from the calculation
-  void storeDerivatives( const unsigned& , MultiValue& myvals, std::vector<double>& );
+  bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_index ) const ;
+/// Build index stores
+//  void buildIndexStores( const unsigned& current, MultiValue& myvals, std::vector<unsigned>& val_index, std::vector<unsigned>& der_index ) const ;
 /// Final step in gathering data
   virtual void finish( const std::vector<double>& buffer );
 /// Is a particular stored value active at the present time
   bool storedValueIsActive( const unsigned& iatom ); 
+/// Set the active values
+  void setActiveValsAndDerivatives( const std::vector<unsigned>& der_index );
 /// Activate indexes (this is used at end of chain rule)
   virtual void activateIndices( ActionWithVessel* ){}
-/// Get the list of indices that we are storing data for
-//  virtual void getIndexList( const unsigned& , const unsigned& , const unsigned& , std::vector<unsigned>& );
 /// Forces on vectors should always be applied elsewhere
   virtual bool applyForce(std::vector<double>&){ return false; }
 };
@@ -154,68 +127,21 @@ bool StoreDataVessel::usingLowMem(){
   return getAction()->lowmem;
 }
 
-// inline
-// void StoreDataVessel::performTask( const unsigned& ivec ){
-// //  if( usingLowMem() ) getAction()->performTask();
-// }
-
-// inline
-// double StoreDataVessel::getComponent( const unsigned& ival, const unsigned& jcomp ){
-//   plumed_dbg_assert( ival<getAction()->getFullNumberOfTasks() && jcomp<vecsize );
-//   return 0.0;    // GAT Broken
-// //  return getBufferElement( ival*(vecsize*nspace) + jcomp*nspace ); 
-// }
-// 
-// inline
-// void StoreDataVessel::setComponent( const unsigned& ival, const unsigned& jcomp, const double& val ){
-//   plumed_dbg_assert( ival<getAction()->getFullNumberOfTasks() && jcomp<vecsize );
-// //  setBufferElement( ival*(vecsize*nspace) + jcomp*nspace, val );
-// }
-
 inline
 unsigned StoreDataVessel::getNumberOfDerivativeSpacesPerComponent() const {
   return nspace;
 }
 
-// inline
-// unsigned StoreDataVessel::getNumberOfDerivatives( const unsigned& ival ){
-//   plumed_dbg_assert( ival<getAction()->getFullNumberOfTasks() );
-//   return active_der[ival];
-// }
-
-// inline
-// unsigned StoreDataVessel::getStoredIndex( const unsigned& ival, const unsigned& jindex ){
-//   plumed_dbg_assert( ival<getAction()->getFullNumberOfTasks() && jindex<active_der[ival] );
-// 
-//   unsigned kder;
-//   if( getAction()->lowmem ) kder = max_lowmem_stash + ival*getAction()->getNumberOfDerivatives();
-//   else kder = getAction()->getFullNumberOfTasks() + ival*(nspace-1); 
-// 
-//   return active_der[kder + jindex];
-// }
-
-// inline
-// double StoreDataVessel::getLocalDerivative( const unsigned& ibuf ){
-//   plumed_dbg_assert( getAction()->lowmem && ibuf<local_derivatives.size() );
-//   return local_derivatives[ibuf];
-// }
-// 
-// inline
-// double StoreDataVessel::getFinalDerivative( const unsigned& ider ) const {
-//   plumed_dbg_assert( ider<final_derivatives.size() );
-//   return final_derivatives[ider];
-// }
-// 
-// inline
-// void StoreDataVessel::setLocalDerivative( const unsigned& ibuf, const double& val ){
-//   plumed_dbg_assert( getAction()->lowmem && ibuf<local_derivatives.size() );
-//   local_derivatives[ibuf]=val;
-// }
-
 inline
 bool StoreDataVessel::storedValueIsActive( const unsigned& iatom ){
-  plumed_dbg_assert( iatom<active_val.size() );
-  return (active_val[iatom]==1);
+  plumed_dbg_assert( iatom<getAction()->getFullNumberOfTasks() );
+  if( !hard_cut ) return true; 
+  return local_buffer[iatom*vecsize*nspace]>wtol;   // (active_val[iatom]==1);
+}
+
+inline
+unsigned StoreDataVessel::getSizeOfDerivativeList() const {
+  return active_der.size();
 }
 
 }
diff --git a/src/vesselbase/Sum.cpp b/src/vesselbase/Sum.cpp
index 08ce670bece4c78d4c9c4a034d6c8711156062c9..485844765c20b02c4a47c938d4cc7231a78bf9d5 100644
--- a/src/vesselbase/Sum.cpp
+++ b/src/vesselbase/Sum.cpp
@@ -31,7 +31,7 @@ public:
   static void reserveKeyword( Keywords& keys );
   Sum( const VesselOptions& da );
   std::string function_description();
-  double calcTransform( const double& val, double& dv );
+  double calcTransform( const double& val, double& dv ) const ;
 };
 
 PLUMED_REGISTER_VESSEL(Sum,"SUM")
@@ -54,7 +54,7 @@ std::string Sum::function_description(){
   return "the sum of all the values"; 
 }
 
-double Sum::calcTransform( const double& val, double& dv ){
+double Sum::calcTransform( const double& val, double& dv ) const {
   dv=1.0; return val;
 }
 
diff --git a/src/vesselbase/Vessel.h b/src/vesselbase/Vessel.h
index 8ff117e6ac33981467b31d4737a403e22ad16928..304669dc9e3b0b11f7ecd175282226ee173d4cf4 100644
--- a/src/vesselbase/Vessel.h
+++ b/src/vesselbase/Vessel.h
@@ -113,7 +113,7 @@ protected:
 /// This returns the whole input line (it is used for less_than/more_than/between)
   std::string getAllInput(); 
 /// Return a pointer to the action we are working in
-  ActionWithVessel* getAction();
+  ActionWithVessel* getAction() const ;
 /// Return the value of the tolerance
   double getTolerance() const ;
 /// Return the value of the neighbor list tolerance
@@ -146,7 +146,7 @@ public:
 /// This is replaced in bridges so we can transform the derivatives
   virtual MultiValue& transformDerivatives( const unsigned& current, MultiValue& myvals, MultiValue& bvals );
 /// Calculate the part of the vessel that is done in the loop
-  virtual bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer )=0;
+  virtual bool calculate( const unsigned& current, MultiValue& myvals, std::vector<double>& buffer, std::vector<unsigned>& der_list ) const = 0;
 /// Complete the calculation once the loop is finished
   virtual void finish( const std::vector<double>& )=0;
 /// Reset the size of the buffers
@@ -233,7 +233,7 @@ double Vessel::getNLTolerance() const {
 }
 
 inline
-ActionWithVessel* Vessel::getAction(){
+ActionWithVessel* Vessel::getAction() const {
   return action;
 }