diff --git a/src/adjmat/TopologyMatrix.cpp b/src/adjmat/TopologyMatrix.cpp
index 9baad49c01..2c6028f4dc 100644
--- a/src/adjmat/TopologyMatrix.cpp
+++ b/src/adjmat/TopologyMatrix.cpp
@@ -127,7 +127,7 @@ double TopologyMatrix::calculateWeight( const Vector& pos1, const Vector& pos2,
   // Now run through all sea atoms
   HistogramBead bead; bead.isNotPeriodic(); bead.setKernelType( kerneltype );
   Vector g1derivf,g2derivf,lderivf; Tensor vir; double binlength = maxbins * binw_mat;
-  MultiValue tvals; tvals.resize( maxbins, myvals.getNumberOfDerivatives() );
+  MultiValue tvals; tvals.resize( maxbins, myvals.getNumberOfDerivatives(), 0 );
   for(unsigned i=0; i<natoms; ++i) {
     // Position of sea atom (this will be the origin)
     Vector d2 = getPosition(i,myvals);
diff --git a/src/colvar/MultiColvarTemplate.h b/src/colvar/MultiColvarTemplate.h
index bd11210de5..d4cfe2393b 100644
--- a/src/colvar/MultiColvarTemplate.h
+++ b/src/colvar/MultiColvarTemplate.h
@@ -42,6 +42,7 @@ class MultiColvarTemplate : public ActionWithVector {
   static void registerKeywords(Keywords&);
   explicit MultiColvarTemplate(const ActionOptions&);
   unsigned getNumberOfDerivatives() override ;
+  unsigned getNumberOfAtomsPerTask() const 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 performTask( const unsigned&, MultiValue& ) const override ;
@@ -127,11 +128,15 @@ void MultiColvarTemplate<T>::addComponentWithDerivatives( const std::string& nam
   std::vector<unsigned> s(1); s[0]=ablocks[0].size(); addComponent( name, s );
 }
 
+template <class T>
+unsigned MultiColvarTemplate<T>::getNumberOfAtomsPerTask() const { 
+  return ablocks.size();
+}
+
 template <class T>
 void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue& myvals ) const {
   // Retrieve the positions
   std::vector<Vector> & fpositions( myvals.getFirstAtomVector() );
-  if( fpositions.size()!=ablocks.size() ) fpositions.resize( ablocks.size() );
   for(unsigned i=0; i<ablocks.size(); ++i) fpositions[i] = getPosition( ablocks[i][task_index] );
   // If we are using pbc make whole
   if( usepbc ) {
@@ -145,19 +150,13 @@ void MultiColvarTemplate<T>::performTask( const unsigned& task_index, MultiValue
     }
   } else if( fpositions.size()==1 ) fpositions[0]=delta(Vector(0.0,0.0,0.0),getPosition( ablocks[0][task_index] ) );
   // Retrieve the masses and charges
-  myvals.resizeTemporyVector(2);
   std::vector<double> & mass( myvals.getTemporyVector(0) );
   std::vector<double> & charge( myvals.getTemporyVector(1) );
-  if( mass.size()!=ablocks.size() ) { mass.resize(ablocks.size()); charge.resize(ablocks.size()); }
   for(unsigned i=0; i<ablocks.size(); ++i) { mass[i]=getMass( ablocks[i][task_index] ); charge[i]=getCharge( ablocks[i][task_index] ); }
   // Make some space to store various things
   std::vector<double> values( getNumberOfComponents() );
   std::vector<Tensor> & virial( myvals.getFirstAtomVirialVector() );
   std::vector<std::vector<Vector> > & derivs( myvals.getFirstAtomDerivativeVector() );
-  if( derivs.size()!=values.size() ) { derivs.resize( values.size() ); virial.resize( values.size() ); }
-  for(unsigned i=0; i<derivs.size(); ++i) {
-    if( derivs[i].size()<ablocks.size() ) derivs[i].resize( ablocks.size() );
-  }
   // Calculate the CVs using the method in the Colvar
   T::calculateCV( mode, mass, charge, fpositions, values, derivs, virial, this );
   for(unsigned i=0; i<values.size(); ++i) myvals.setValue( i, values[i] );
diff --git a/src/core/ActionWithVector.cpp b/src/core/ActionWithVector.cpp
index 8eed7e7771..a1abec5c06 100644
--- a/src/core/ActionWithVector.cpp
+++ b/src/core/ActionWithVector.cpp
@@ -222,7 +222,9 @@ void ActionWithVector::runAllTasks() {
     std::vector<double> omp_buffer;
     const unsigned t=OpenMP::getThreadNum();
     if( nt>1 ) omp_buffer.resize( bufsize, 0.0 );
-    if( myvals[t].getNumberOfValues()!=getNumberOfComponents() || myvals[t].getNumberOfDerivatives()!=nderivatives ) myvals[t].resize( getNumberOfComponents(), nderivatives );
+    if( myvals[t].getNumberOfValues()!=getNumberOfComponents() || myvals[t].getNumberOfDerivatives()!=nderivatives || myvals[t].getAtomVector().size()!=getNumberOfAtomsPerTask() ) { 
+        myvals[t].resize( getNumberOfComponents(), nderivatives, getNumberOfAtomsPerTask() );
+    }
     myvals[t].clearAll();
 
     #pragma omp for nowait
@@ -379,7 +381,9 @@ bool ActionWithVector::checkForForces() {
       if( omp_forces[t].size()!=forcesForApply.size() ) omp_forces[t].resize( forcesForApply.size(), 0.0 );
       else omp_forces[t].assign( forcesForApply.size(), 0.0 );
     }
-    if( myvals[t].getNumberOfValues()!=getNumberOfComponents() || myvals[t].getNumberOfDerivatives()!=nderiv ) myvals[t].resize( getNumberOfComponents(), nderiv );
+    if( myvals[t].getNumberOfValues()!=getNumberOfComponents() || myvals[t].getNumberOfDerivatives()!=nderiv || myvals[t].getAtomVector().size()!=getNumberOfAtomsPerTask() ) {
+        myvals[t].resize( getNumberOfComponents(), nderiv, getNumberOfAtomsPerTask() );
+    }
     myvals[t].clearAll();
 
     #pragma omp for nowait
diff --git a/src/core/ActionWithVector.h b/src/core/ActionWithVector.h
index fa6e63c0c7..cf0bb6902c 100644
--- a/src/core/ActionWithVector.h
+++ b/src/core/ActionWithVector.h
@@ -84,6 +84,8 @@ class ActionWithVector:
 /// Check if a mask has been set
   int getNumberOfMasks() const ;
   void calculateNumericalDerivatives(ActionWithValue* av) override;
+/// This is for resizing the task list
+  virtual unsigned getNumberOfAtomsPerTask() const { return 0; }
 /// Turn off the calculation of the derivatives during the forward pass through a calculation
   bool doNotCalculateDerivatives() const override ;
 /// Get the list of tasks that are active
diff --git a/src/refdist/MatrixProductDiagonal.cpp b/src/refdist/MatrixProductDiagonal.cpp
index 1679ac53df..f9b2e5b2a8 100644
--- a/src/refdist/MatrixProductDiagonal.cpp
+++ b/src/refdist/MatrixProductDiagonal.cpp
@@ -127,7 +127,7 @@ void MatrixProductDiagonal::performTask( const unsigned& task_index, MultiValue&
 void MatrixProductDiagonal::calculate() {
   if( getPntrToArgument(1)->getRank()==1 ) {
     unsigned nder = getNumberOfDerivatives();
-    MultiValue myvals; myvals.resize( 1, nder ); performTask( 0, myvals );
+    MultiValue myvals; myvals.resize( 1, nder, 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 24e1153a96..680ed0b3ef 100644
--- a/src/tools/MultiValue.cpp
+++ b/src/tools/MultiValue.cpp
@@ -24,12 +24,17 @@
 
 namespace PLMD {
 
-void MultiValue::resize( const size_t& nvals, const size_t& nder ) {
+void MultiValue::resize( const size_t& nvals, const size_t& nder, const size_t& natoms ) {
   if( values.size()==nvals && nderivatives>nder ) return;
   values.resize(nvals); nderivatives=nder; derivatives.resize( nvals*nder );
   hasderiv.resize(nvals*nder,false); nactive.resize(nvals); active_list.resize(nvals*nder);
   matrix_force_stash.resize(nder,0);
   matrix_row_nderivatives=0; matrix_row_derivative_indices.resize(nder); atLeastOneSet=false;
+  if( natoms>0 ) {
+      tmp_vectors.resize(2); tmp_vectors[0].resize(natoms); tmp_vectors[1].resize(natoms);
+      tmp_atom_der.resize(nvals); tmp_atom_virial.resize(nvals); tmp_atoms.resize(natoms);
+      for(unsigned i=0; i<nvals; ++i) tmp_atom_der[i].resize(natoms); 
+  }
 }
 
 void MultiValue::clearAll() {
diff --git a/src/tools/MultiValue.h b/src/tools/MultiValue.h
index 73df0cbb92..1991e119f4 100644
--- a/src/tools/MultiValue.h
+++ b/src/tools/MultiValue.h
@@ -63,7 +63,7 @@ class MultiValue {
   std::vector<std::vector<double> > tmp_vectors;
 public:
   MultiValue() : task_index(0), task2_index(0), nderivatives(0), atLeastOneSet(false), vector_call(false), nindices(0), nsplit(0), matrix_row_nderivatives(0) {}
-  void resize( const std::size_t& nvals, const std::size_t& nder );
+  void resize( const std::size_t& nvals, const std::size_t& nder, const std::size_t& natoms );
 /// Set the task index prior to the loop
   void setTaskIndex( const std::size_t& tindex );
 ///
@@ -91,7 +91,6 @@ class MultiValue {
   std::vector<std::vector<Vector> >& getFirstAtomDerivativeVector();
   const std::vector<std::vector<Vector> >& getConstFirstAtomDerivativeVector() const ;
   std::vector<Tensor>& getFirstAtomVirialVector();
-  void resizeTemporyVector(const unsigned& n );
   std::vector<double>& getTemporyVector(const unsigned& ind );
 ///
   bool inVectorCall() const ;
@@ -322,11 +321,6 @@ double MultiValue::getStashedMatrixForce( const unsigned& jind ) const {
   return matrix_force_stash[jind];
 }
 
-inline
-void MultiValue::resizeTemporyVector(const unsigned& n ) {
-  if( n>tmp_vectors.size() ) tmp_vectors.resize(n);
-}
-
 inline
 std::vector<double>& MultiValue::getTemporyVector(const unsigned& ind ) {
   plumed_dbg_assert( ind<tmp_vectors.size() );