Small speedup when calculating matrices without using derivatives

Derivatives do not need to be cleared in these cases so we avoid clearing them now.

src/adjmat/AdjacencyMatrixBase.cpp

@@ -52,6 +52,7 @@ AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
   neighbour_list_updated(false),
   linkcells(comm),
   threecells(comm),
+  maxcol(0),
   natoms_per_list(0)
 {
   std::vector<unsigned> shape(2); std::vector<AtomNumber> t; parseAtomList("GROUP", t );
@@ -291,9 +292,9 @@ void AdjacencyMatrixBase::setupForTask( const unsigned& current, std::vector<uns
 void AdjacencyMatrixBase::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
   Vector zero; zero.zero(); plumed_dbg_assert( index2<myvals.getAtomVector().size() );
   double weight = calculateWeight( zero, myvals.getAtomVector()[index2], myvals.getNumberOfIndices()-myvals.getSplitIndex(), myvals );
-  if( fabs(weight)<epsilon ) return;
-
   unsigned w_ind = getConstPntrToComponent(0)->getPositionInStream(); myvals.setValue( w_ind, weight );
+  if( fabs(weight)<epsilon ) { myvals.setValue( w_ind, 0 ); return; }
+
   if( !doNotCalculateDerivatives() ) {
     // Update dynamic list indices for central atom
     myvals.updateIndex( w_ind, 3*index1+0 ); myvals.updateIndex( w_ind, 3*index1+1 ); myvals.updateIndex( w_ind, 3*index1+2 );

src/core/ActionWithMatrix.cpp

@@ -33,7 +33,8 @@ ActionWithMatrix::ActionWithMatrix(const ActionOptions&ao):
   ActionWithVector(ao),
   next_action_in_chain(NULL),
   matrix_to_do_before(NULL),
-  matrix_to_do_after(NULL)
+  matrix_to_do_after(NULL),
+  clearOnEachCycle(true)
 {
 }
 
@@ -62,6 +63,14 @@ void ActionWithMatrix::setupStreamedComponents( const std::string& headstr, unsi
     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;
+  for(int i=0; i<getNumberOfComponents(); ++i) {
+    const Value* myval=getConstPntrToComponent(i);
+    if( myval->getRank()==2 && !myval->hasDerivatives() ) { clearOnEachCycle = true; break; }
+  }
+  // Turn off clearing of derivatives if we have only the values of adjacency matrices
+  if( doNotCalculateDerivatives() && isAdjacencyMatrix() ) clearOnEachCycle = false;
 }
 
 void ActionWithMatrix::finishChainBuild( ActionWithVector* act ) {
@@ -222,7 +231,7 @@ void ActionWithMatrix::gatherForcesOnStoredValue( const Value* myval, const unsi
 }
 
 void ActionWithMatrix::clearMatrixElements( MultiValue& myvals ) const {
-  if( isActive() ) {
+  if( isActive() && clearOnEachCycle ) {
     for(int i=0; i<getNumberOfComponents(); ++i) {
       const Value* myval=getConstPntrToComponent(i);
       if( myval->getRank()==2 && !myval->hasDerivatives() ) myvals.clearDerivatives( myval->getPositionInStream() );

src/core/ActionWithMatrix.h

@@ -44,6 +44,8 @@ class ActionWithMatrix : public ActionWithVector {
 /// 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:
+/// This turns off derivative clearing for contact matrix if we are not storing derivatives
+  bool clearOnEachCycle;
 /// Does the matrix chain continue on from this action
   bool matrixChainContinues() const ;
 /// This returns the jelem th element of argument ic