Optimised adding of forces in AdjacencyMatrixBase by using the fact t…

…hat we already worked out that some of the elements of the matrix are zero and that the forces on these elements are zero

src/adjmat/AdjacencyMatrixBase.cpp

@@ -267,7 +267,15 @@ void AdjacencyMatrixBase::setupForTask( const unsigned& current, std::vector<uns
   if( indices.size()!=(1+ablocks.size()+threeblocks.size()) ) indices.resize( 1+ablocks.size()+threeblocks.size() );
 
   // Now get the positions
-  unsigned natoms=retrieveNeighbours( current, indices );
+  unsigned natoms; const Value* myval = getConstPntrToComponent(0);
+  if( myval->forcesWereAdded() ) {
+    natoms = 0; indices[natoms] = current; natoms++;
+    unsigned nelements = myval->getRowLength(current), startr = current*myval->getNumberOfColumns();
+    for(unsigned j=0; j<nelements; ++j ) {
+      unsigned jind = myval->getRowIndex( current, j );
+      if( ablocks[jind]!=indices[0] && fabs( myval->getForce( startr + j ) )>epsilon ) { indices[natoms] = ablocks[jind]; natoms++; }
+    }
+  } else natoms=retrieveNeighbours( current, indices );
   unsigned ntwo_atoms=natoms; myvals.setSplitIndex( ntwo_atoms );
 
   // Now retrieve everything for the third atoms

src/valtools/VStack.cpp

@@ -140,7 +140,7 @@ void VStack::performTask( const std::string& controller, const unsigned& index1,
 }
 
 void VStack::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
-  unsigned matind = myval->getPositionInMatrixStash(); const std::vector<unsigned>& mat_indices( myvals.getMatrixRowDerivativeIndices( matind ) );
+  unsigned matind = myval->getPositionInMatrixStash();
   for(unsigned i=0; i<forces.size(); ++i) forces[i] += myvals.getStashedMatrixForce( matind, i );
 }