diff --git a/src/multicolvar/ActionVolume.cpp b/src/multicolvar/ActionVolume.cpp
index d5e515af46f2935a272cdddad54463ab99132d36..9b65d698f6cc8a28825b75a6682fcc61eb5a3556 100644
--- a/src/multicolvar/ActionVolume.cpp
+++ b/src/multicolvar/ActionVolume.cpp
@@ -91,6 +91,8 @@ lastUpdate(0)
 
   // Now set up the bridging vessel (has to be done this way for internal arrays to be resized properly)
   addDependency(mycolv); myBridgeVessel = mycolv->addBridgingVessel( this );
+  // Setup a dynamic list 
+  resizeLocalArrays();
 }
 
 void ActionVolume::requestAtoms( const std::vector<AtomNumber>& atoms ){
@@ -106,7 +108,8 @@ void ActionVolume::doJobsRequiredBeforeTaskList(){
 }
 
 void ActionVolume::prepare(){
-  if( contributorsAreUnlocked ) lockContributors();
+  bool updatetime=false;
+  if( contributorsAreUnlocked ){ updatetime=true; lockContributors(); }
   if( updateFreq>0 && (getStep()-lastUpdate)>=updateFreq ){
       if( !mycolv->isDensity() ){
           mycolv->taskList.activateAll();
@@ -117,10 +120,19 @@ void ActionVolume::prepare(){
           plumed_massert( mycolv->contributorsAreUnlocked, "contributors are not unlocked in base multicolvar" );
       }
       unlockContributors(); 
+      updatetime=true;
   }
+  if( updatetime ) resizeLocalArrays();
+}
+
+void ActionVolume::resizeLocalArrays(){
+  activeAtoms.clear();
+  for(unsigned i=0;i<mycolv->getNumberOfAtoms();++i) activeAtoms.addIndexToList(i);
+  activeAtoms.deactivateAll();
 }
 
 void ActionVolume::performTask( const unsigned& j ){
+  activeAtoms.deactivateAll(); // Currently no atoms are active so deactivate them all
   Vector catom_pos=mycolv->retrieveCentralAtomPos();
 
   double weight; Vector wdf; 
@@ -131,8 +143,8 @@ void ActionVolume::performTask( const unsigned& j ){
      unsigned nder=getNumberOfDerivatives();
      setElementValue( 1, weight ); setElementValue( 0, 1.0 ); 
      for(unsigned i=0;i<mycolv->atomsWithCatomDer.getNumberActive();++i){
-        unsigned n=mycolv->atomsWithCatomDer[i];
-        unsigned nx=nder + 3*n; 
+        unsigned n=mycolv->atomsWithCatomDer[i], nx=nder + 3*n;
+        activeAtoms.activate(n);
         addElementDerivative( nx+0, mycolv->getCentralAtomDerivative(n, 0, wdf ) );
         addElementDerivative( nx+1, mycolv->getCentralAtomDerivative(n, 1, wdf ) );
         addElementDerivative( nx+2, mycolv->getCentralAtomDerivative(n, 2, wdf ) );
@@ -141,10 +153,11 @@ void ActionVolume::performTask( const unsigned& j ){
      // Copy derivatives of the colvar and the value of the colvar
      double colv=mycolv->getElementValue(0); setElementValue( 0, colv );
      for(unsigned i=0;i<mycolv->atoms_with_derivatives.getNumberActive();++i){
-        unsigned n=3*mycolv->atoms_with_derivatives(i);
-        addElementDerivative( n, mycolv->getElementDerivative(n) ); n++;
-        addElementDerivative( n, mycolv->getElementDerivative(n) ); n++;
-        addElementDerivative( n, mycolv->getElementDerivative(n) ); 
+        unsigned n=mycolv->atoms_with_derivatives[i], nx=3*n;
+        activeAtoms.activate(n);
+        addElementDerivative( nx+0, mycolv->getElementDerivative(nx+0) );
+        addElementDerivative( nx+1, mycolv->getElementDerivative(nx+1) );
+        addElementDerivative( nx+2, mycolv->getElementDerivative(nx+2) ); 
      }
      unsigned nvir=3*mycolv->getNumberOfAtoms();
      for(unsigned i=0;i<9;++i){ 
@@ -158,10 +171,11 @@ void ActionVolume::performTask( const unsigned& j ){
      // Add derivatives of weight if we have a weight
      if( mycolv->weightHasDerivatives ){
         for(unsigned i=0;i<mycolv->atoms_with_derivatives.getNumberActive();++i){
-           unsigned n=3*mycolv->atoms_with_derivatives(i);
-           addElementDerivative( nder+n, weight*mycolv->getElementDerivative(nder+n) ); n++;
-           addElementDerivative( nder+n, weight*mycolv->getElementDerivative(nder+n) ); n++;
-           addElementDerivative( nder+n, weight*mycolv->getElementDerivative(nder+n) );
+           unsigned n=mycolv->atoms_with_derivatives[i], nx=nder + 3*n;
+           activeAtoms.activate(n); 
+           addElementDerivative( nx+0, weight*mycolv->getElementDerivative(nx+0) ); 
+           addElementDerivative( nx+1, weight*mycolv->getElementDerivative(nx+1) ); 
+           addElementDerivative( nx+2, weight*mycolv->getElementDerivative(nx+2) );
         } 
         unsigned nwvir=3*mycolv->getNumberOfAtoms();
         for(unsigned i=0;i<9;++i){
@@ -171,13 +185,60 @@ void ActionVolume::performTask( const unsigned& j ){
 
      // Add derivatives of central atoms
      for(unsigned i=0;i<mycolv->atomsWithCatomDer.getNumberActive();++i){
-         unsigned n=mycolv->atomsWithCatomDer[i];
-         unsigned nx=nder+3*n; 
+         unsigned n=mycolv->atomsWithCatomDer[i], nx=nder+3*n;
+         activeAtoms.activate(n); 
          addElementDerivative( nx+0, ww*mycolv->getCentralAtomDerivative(n, 0, wdf ) );
          addElementDerivative( nx+1, ww*mycolv->getCentralAtomDerivative(n, 1, wdf ) );
          addElementDerivative( nx+2, ww*mycolv->getCentralAtomDerivative(n, 2, wdf ) );
      }
   }
+  activeAtoms.updateActiveMembers();
+}
+
+void ActionVolume::mergeDerivatives( const unsigned& ider, const double& df ){
+  unsigned vstart=getNumberOfDerivatives()*ider;
+  // Merge atom derivatives
+  for(unsigned i=0;i<activeAtoms.getNumberActive();++i){
+     unsigned iatom=3*activeAtoms[i];
+     accumulateDerivative( iatom, df*getElementDerivative(vstart+iatom) ); iatom++;
+     accumulateDerivative( iatom, df*getElementDerivative(vstart+iatom) ); iatom++;
+     accumulateDerivative( iatom, df*getElementDerivative(vstart+iatom) );
+  }
+  // Merge virial derivatives
+  unsigned nvir=3*mycolv->getNumberOfAtoms();
+  for(unsigned j=0;j<9;++j){
+     accumulateDerivative( nvir, df*getElementDerivative(vstart+nvir) ); nvir++;
+  }
+  // Merge local atom derivatives
+  for(unsigned j=0;j<getNumberOfAtoms();++j){
+     accumulateDerivative( nvir, df*getElementDerivative(vstart+nvir) ); nvir++;
+     accumulateDerivative( nvir, df*getElementDerivative(vstart+nvir) ); nvir++;
+     accumulateDerivative( nvir, df*getElementDerivative(vstart+nvir) ); nvir++;
+  }
+  plumed_dbg_assert( nvir==getNumberOfDerivatives() );
+}
+
+void ActionVolume::clearDerivativesAfterTask( const unsigned& ider ){
+  unsigned vstart=getNumberOfDerivatives()*ider;
+  // Clear atom derivatives
+  for(unsigned i=0;i<activeAtoms.getNumberActive();++i){
+     unsigned iatom=vstart+3*activeAtoms[i];
+     setElementDerivative( iatom, 0.0 ); iatom++;
+     setElementDerivative( iatom, 0.0 ); iatom++;
+     setElementDerivative( iatom, 0.0 );
+  }
+  // Clear virial contribution
+  unsigned nvir=vstart+3*mycolv->getNumberOfAtoms();
+  for(unsigned j=0;j<9;++j){
+     setElementDerivative( nvir, 0.0 ); nvir++;
+  }
+  // Clear derivatives of local atoms
+  for(unsigned j=0;j<getNumberOfAtoms();++j){
+     setElementDerivative( nvir, 0.0 ); nvir++;
+     setElementDerivative( nvir, 0.0 ); nvir++;
+     setElementDerivative( nvir, 0.0 ); nvir++;
+  }
+  plumed_dbg_assert( nvir==getNumberOfDerivatives() );
 }
 
 void ActionVolume::calculateNumericalDerivatives( ActionWithValue* a ){
diff --git a/src/multicolvar/ActionVolume.h b/src/multicolvar/ActionVolume.h
index 1f896b90cff8c1ce4d3c8be8c0686e8e59ebb3a5..d1a8ee86ca66e9d1b487f2fd01853d12fe8e3fc8 100644
--- a/src/multicolvar/ActionVolume.h
+++ b/src/multicolvar/ActionVolume.h
@@ -111,6 +111,10 @@ public:
 /// Forces here are applied through the bridge
   void applyBridgeForces( const std::vector<double>& bb );
   void apply(){};
+/// These routines replace the virtual routines in ActionWithVessel for 
+/// code optimization
+  void mergeDerivatives( const unsigned& ider, const double& df );
+  void clearDerivativesAfterTask( const unsigned& ider );
 };
 
 inline
@@ -174,15 +178,15 @@ void ActionVolume::addBoxDerivatives( const Tensor& vir ){
   double pref=mycolv->getElementValue(1);
   if( not_in ) pref*=-1;
   unsigned nstart = getNumberOfDerivatives() + mycolv->getNumberOfDerivatives() - 9;
-  addElementDerivative( nstart + 0, pref*vir(0,0) ); // printf("VIR %f ",vir(0,0) ); 
-  addElementDerivative( nstart + 1, pref*vir(0,1) ); // printf("%f ",vir(0,1) ); 
-  addElementDerivative( nstart + 2, pref*vir(0,2) ); // printf("%f ",vir(0,2) ); 
-  addElementDerivative( nstart + 3, pref*vir(1,0) ); // printf("%f ",vir(1,0) ); 
-  addElementDerivative( nstart + 4, pref*vir(1,1) ); // printf("%f ",vir(1,1) ); 
-  addElementDerivative( nstart + 5, pref*vir(1,2) ); // printf("%f ",vir(1,2) ); 
-  addElementDerivative( nstart + 6, pref*vir(2,0) ); // printf("%f ",vir(2,0) ); 
-  addElementDerivative( nstart + 7, pref*vir(2,1) ); // printf("%f ",vir(2,1) ); 
-  addElementDerivative( nstart + 8, pref*vir(2,2) ); // printf("%f \n",vir(2,2) ); 
+  addElementDerivative( nstart + 0, pref*vir(0,0) );  
+  addElementDerivative( nstart + 1, pref*vir(0,1) ); 
+  addElementDerivative( nstart + 2, pref*vir(0,2) ); 
+  addElementDerivative( nstart + 3, pref*vir(1,0) ); 
+  addElementDerivative( nstart + 4, pref*vir(1,1) ); 
+  addElementDerivative( nstart + 5, pref*vir(1,2) ); 
+  addElementDerivative( nstart + 6, pref*vir(2,0) ); 
+  addElementDerivative( nstart + 7, pref*vir(2,1) ); 
+  addElementDerivative( nstart + 8, pref*vir(2,2) ); 
 }
 
 }
diff --git a/src/vesselbase/BridgeVessel.cpp b/src/vesselbase/BridgeVessel.cpp
index 22a5fb8839eb52c814b26c9678b5226af610cbe6..1f2b5ede383dd0104e4cdaa301ab131345d74656 100644
--- a/src/vesselbase/BridgeVessel.cpp
+++ b/src/vesselbase/BridgeVessel.cpp
@@ -57,8 +57,7 @@ void BridgeVessel::prepare(){
 }
 
 bool BridgeVessel::calculate(){
-  unsigned j;
-  myOutputAction->performTask(j);
+  myOutputAction->performTask( getAction()->current );
   if( myOutputAction->thisval[1]<myOutputAction->getTolerance() ){
       myOutputAction->clearAfterTask();
       return ( !myOutputAction->contributorsAreUnlocked || myOutputAction->thisval[1]>=myOutputAction->getNLTolerance() );