Small fixes that correct the calculation forces when there is a funct…

…ion of a matrix in the apply loop

src/core/ActionWithMatrix.cpp

@@ -162,9 +162,9 @@ void ActionWithMatrix::runTask( const std::string& controller, const unsigned& c
       unsigned matindex = myval->getPositionInMatrixStash(), col_stash_index = colno;
       if( colno>=myval->getShape()[0] ) col_stash_index = colno - myval->getShape()[0];
       if( myval->forcesWereAdded() ) {
-        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 );
+        unsigned sind = myval->getPositionInStream(); double fforce;
+        if( ncols<myval->getShape()[1] ) fforce = myval->getForce( myvals.getTaskIndex()*ncols + myval->matrix_bookeeping[current*(1+ncols)] );
+        else 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 ) );
         }

src/function/FunctionOfMatrix.h

@@ -233,9 +233,9 @@ unsigned FunctionOfMatrix<T>::getNumberOfColumns() const {
 template <class T>
 void FunctionOfMatrix<T>::setupForTask( const unsigned& task_index, std::vector<unsigned>& indices, MultiValue& myvals ) const {
   for(unsigned i=0; i<getNumberOfArguments(); ++i) plumed_assert( getPntrToArgument(i)->getRank()==2 );
-  unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(0)->getShape()[1];
+  unsigned start_n = getPntrToArgument(0)->getShape()[0], size_v = getPntrToArgument(0)->getRowLength(task_index);
   if( indices.size()!=size_v+1 ) indices.resize( size_v+1 );
-  for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + i;
+  for(unsigned i=0; i<size_v; ++i) indices[i+1] = start_n + getPntrToArgument(0)->getRowIndex(task_index, i);
   myvals.setSplitIndex( size_v + 1 );
 }
 
@@ -337,7 +337,7 @@ void FunctionOfMatrix<T>::performTask( const std::string& controller, const unsi
       if( getPntrToArgument(j)->getRank()==2 ) {
         for(int i=0; i<getNumberOfComponents(); ++i) {
           unsigned ostrn=getConstPntrToComponent(i)->getPositionInStream();
-          unsigned myind = base + getPntrToArgument(j)->getShape()[1]*index1 + ind2;
+          unsigned myind = base + getPntrToArgument(j)->getNumberOfColumns()*index1 + ind2;
           myvals.addDerivative( ostrn, myind, derivatives(i,j) );
           myvals.updateIndex( ostrn, myind );
         }
@@ -348,7 +348,7 @@ void FunctionOfMatrix<T>::performTask( const std::string& controller, const unsi
           myvals.updateIndex( ostrn, base );
         }
       }
-      base += getPntrToArgument(j)->getNumberOfValues();
+      base += getPntrToArgument(j)->getNumberOfStoredValues();
     }
   }
 }
@@ -413,9 +413,9 @@ void FunctionOfMatrix<T>::runEndOfRowJobs( const unsigned& ind, const std::vecto
       if( mat_indices.size()<nderivatives ) mat_indices.resize( nderivatives ); unsigned matderbase = 0;
       for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
         if( getPntrToArgument(i)->getRank()==0 ) continue ;
-        unsigned ss = getPntrToArgument(i)->getShape()[1]; unsigned tbase = matderbase + ss*myvals.getTaskIndex();
+        unsigned ss = getPntrToArgument(i)->getNumberOfColumns(); unsigned tbase = matderbase + ss*myvals.getTaskIndex();
         for(unsigned k=0; k<ss; ++k) mat_indices[ntot_mat + k] = tbase + k;
-        ntot_mat += ss; matderbase += getPntrToArgument(i)->getNumberOfValues();
+        ntot_mat += ss; matderbase += getPntrToArgument(i)->getNumberOfStoredValues();
       }
       myvals.setNumberOfMatrixRowDerivatives( nmat, ntot_mat );
     }