Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    plumed is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "ActionWithVector.h"
      23             : #include "ActionWithMatrix.h"
      24             : #include "PlumedMain.h"
      25             : #include "ActionSet.h"
      26             : #include "tools/OpenMP.h"
      27             : #include "tools/Communicator.h"
      28             : #include <cstring>
      29             : #include <algorithm>
      30             : 
      31             : namespace PLMD {
      32             : 
      33        5420 : void ActionWithVector::registerKeywords( Keywords& keys ) {
      34        5420 :   Action::registerKeywords( keys );
      35        5420 :   ActionAtomistic::registerKeywords( keys );
      36        5420 :   ActionWithValue::registerKeywords( keys );
      37        5420 :   keys.remove("NUMERICAL_DERIVATIVES");
      38        5420 :   ActionWithArguments::registerKeywords( keys );
      39        5420 : }
      40             : 
      41        3953 : ActionWithVector::ActionWithVector(const ActionOptions&ao):
      42             :   Action(ao),
      43             :   ActionAtomistic(ao),
      44             :   ActionWithValue(ao),
      45        3953 :   ActionWithArguments(ao) {
      46        3953 :   if( !keywords.exists("MASKED_INPUT_ALLOWED") ) {
      47        4104 :     for(unsigned i=0; i<getNumberOfArguments(); ++i) {
      48        2328 :       ActionWithVector* av = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
      49        2328 :       if( av && av->getNumberOfMasks()>=0 ) {
      50          70 :         nmask=0;
      51             :       }
      52             :     }
      53             :   }
      54             : 
      55        3953 :   if( keywords.exists("MASK") ) {
      56        3023 :     std::vector<Value*> mask;
      57        6046 :     parseArgumentList("MASK",mask);
      58        3023 :     if( mask.size()>0 ) {
      59         245 :       if( nmask>=0 && getNumberOfArguments()==1 ) {
      60          58 :         ActionWithVector* av=dynamic_cast<ActionWithVector*>( getPntrToArgument(0)->getPntrToAction() );
      61          58 :         plumed_massert( av, "input should be a vector from ActionWithVector" );
      62          58 :         unsigned j=0, nargs = av->getNumberOfArguments();
      63         118 :         for(unsigned i=nargs-av->nmask; i<nargs; ++i) {
      64          60 :           if( av->getPntrToArgument(i)!=mask[j] ) {
      65           0 :             error("the masks in subsequent actions do not match");
      66             :           }
      67          60 :           j++;
      68             :         }
      69             :       }
      70         245 :       if( getNumberOfArguments()>0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")==std::string::npos && getPntrToArgument(0)->hasDerivatives() ) {
      71           0 :         error("input for mask should be vector or matrix");
      72         245 :       } else if( mask[0]->getRank()==2 ) {
      73          99 :         if( mask.size()>1 ) {
      74           0 :           error("MASK should only have one argument");
      75             :         }
      76          99 :         log.printf("  only computing elements of matrix that correspond to non-zero elements of matrix %s \n", mask[0]->getName().c_str() );
      77         146 :       } else if( mask[0]->getRank()==1 ) {
      78         146 :         log.printf("  only computing elements of vector that correspond to non-zero elements of vectors %s", mask[0]->getName().c_str() );
      79         150 :         for(unsigned i=1; i<mask.size(); ++i) {
      80           4 :           if( mask[i]->getRank()!=1 ) {
      81           0 :             log.printf("\n");
      82           0 :             error("input to mask should be vector");
      83             :           }
      84           4 :           log.printf(", %s", mask[i]->getName().c_str() );
      85             :         }
      86         146 :         log.printf("\n");
      87             :       }
      88         245 :       std::vector<Value*> allargs( getArguments() );
      89         245 :       nmask=mask.size();
      90         494 :       for(unsigned i=0; i<mask.size(); ++i) {
      91         249 :         allargs.push_back( mask[i] );
      92             :       }
      93         245 :       requestArguments( allargs );
      94             :     }
      95             :   }
      96        3953 : }
      97             : 
      98      218559 : void ActionWithVector::lockRequests() {
      99             :   ActionAtomistic::lockRequests();
     100             :   ActionWithArguments::lockRequests();
     101      218559 : }
     102             : 
     103      218559 : void ActionWithVector::unlockRequests() {
     104             :   ActionAtomistic::unlockRequests();
     105             :   ActionWithArguments::unlockRequests();
     106      218559 : }
     107             : 
     108           0 : void ActionWithVector::calculateNumericalDerivatives(ActionWithValue* av) {
     109           0 :   plumed_merror("cannot calculate numerical derivative for action " + getName() + " with label " + getLabel() );
     110             : }
     111             : 
     112      219708 : void ActionWithVector::prepare() {
     113      219708 :   active_tasks.resize(0);
     114      219708 : }
     115             : 
     116    16536426 : int ActionWithVector::checkTaskIsActive( const unsigned& itask ) const {
     117    16536426 :   unsigned nargs = getNumberOfArguments();
     118    16536426 :   if( nargs==0 ) {
     119             :     return 1;
     120    13399857 :   } else if( nmask>0 ) {
     121     7316417 :     for(unsigned j=nargs-nmask; j<nargs; ++j) {
     122             :       Value* myarg = getPntrToArgument(j);
     123     3857927 :       if( myarg->getRank()==1 && !myarg->hasDerivatives() ) {
     124     3595352 :         if( fabs(myarg->get(itask))>0.0 ) {
     125             :           return 1;
     126             :         }
     127      262575 :       } else if( myarg->getRank()==2 && !myarg->hasDerivatives() ) {
     128      262575 :         if( myarg->getRowLength(itask)>0 ) {
     129             :           return 1;
     130             :         }
     131             :       } else {
     132           0 :         plumed_merror("only matrices and vectors should be used as masks");
     133             :       }
     134             :     }
     135             :   } else {
     136    11531424 :     for(unsigned i=0; i<nargs; ++i) {
     137             :       Value* myarg = getPntrToArgument(i);
     138    11178958 :       if( !myarg->isDerivativeZeroWhenValueIsZero() ) {
     139             :         return 1;
     140             :       }
     141             : 
     142     8140729 :       if( myarg->getRank()==0 ) {
     143             :         return 1;
     144     8140729 :       } else if( myarg->getRank()==1 && !myarg->hasDerivatives() ) {
     145     7375841 :         if( fabs(myarg->get(itask))>0.0 ) {
     146             :           return 1;
     147             :         }
     148      764888 :       } else if( myarg->getRank()==2 && !myarg->hasDerivatives() ) {
     149        3287 :         const unsigned ncol = myarg->getRowLength(itask);
     150        3287 :         const unsigned base = itask*myarg->getNumberOfColumns();
     151        3329 :         for(unsigned k=0; k<ncol; ++k) {
     152        3329 :           if( fabs(myarg->get(base+k,false))>0.0 ) {
     153             :             return 1;
     154             :           }
     155             :         }
     156             :       } else if( myarg->getRank()>0 ) {
     157             :         return 1;
     158             :       } else {
     159             :         plumed_merror("should not be in action " + getName() );
     160             :       }
     161             :     }
     162             :   }
     163             :   return -1;
     164             : }
     165             : 
     166      336067 : std::vector<unsigned>& ActionWithVector::getListOfActiveTasks() {
     167      336067 :   if( active_tasks.size()>0 ) {
     168      156664 :     return active_tasks;
     169             :   }
     170      179403 :   unsigned ntasks=0;
     171      179403 :   getNumberOfTasks( ntasks );
     172             : 
     173      179403 :   active_tasks.resize(0);
     174      179403 :   active_tasks.reserve(ntasks);
     175    21138395 :   for(unsigned i=0; i<ntasks; ++i) {
     176    20958992 :     if( checkTaskIsActive(i)>0 ) {
     177             : //no resize are triggered, since we have reserved the number of tasks
     178    15937578 :       active_tasks.push_back(i);
     179             :     }
     180             :   }
     181      179403 :   return active_tasks;
     182             : }
     183             : 
     184      137260 : void ActionWithVector::getAtomicInputData( std::vector<double>& inputdata ) const {
     185             :   plumed_dbg_assert( getNumberOfAtoms()==0 );
     186      137260 : }
     187             : 
     188           0 : void ActionWithVector::getAtomicInputData( std::vector<float>& inputdata ) const {
     189             :   plumed_dbg_assert( getNumberOfAtoms()==0 );
     190           0 : }
     191             : 
     192       58705 : void ActionWithVector::getArgumentInputData( std::vector<double>& inputdata ) const {
     193       58705 :   unsigned nargs = getNumberOfArguments();
     194       58705 :   unsigned nmasks=getNumberOfMasks();
     195             :   // getNumberOfMasks(); returns nmask, that it is an int
     196             :   // nmasks cant be <0 (it is unsigned), so I check nmask for that
     197       58705 :   if( nargs>=nmasks && nmask>0 ) {
     198          76 :     nargs = nargs - nmasks;
     199             :   }
     200             : 
     201             :   std::size_t total_args = 0;
     202      179877 :   for(unsigned i=0; i<nargs; ++i) {
     203      121172 :     total_args += getPntrToArgument(i)->getNumberOfStoredValues();
     204             :   }
     205             : 
     206       58705 :   if( inputdata.size()!=total_args ) {
     207         939 :     inputdata.resize( total_args );
     208             :   }
     209             : 
     210             :   total_args = 0;
     211      179877 :   for(unsigned i=0; i<nargs; ++i) {
     212             :     Value* myarg = getPntrToArgument(i);
     213      121172 :     total_args+= myarg->assignValues(View<double> {inputdata.data()+total_args,inputdata.size()-total_args});
     214             :   }
     215       58705 : }
     216             : 
     217           0 : void ActionWithVector::getArgumentInputData( std::vector<float>& inputdata ) const {
     218           0 :   unsigned nargs = getNumberOfArguments();
     219             :   int nmasks=getNumberOfMasks();
     220           0 :   if( nargs>=static_cast<unsigned>(nmasks) && nmasks>0 ) {
     221           0 :     nargs = nargs - nmasks;
     222             :   }
     223             : 
     224             :   std::size_t total_args = 0;
     225           0 :   for(unsigned i=0; i<nargs; ++i) {
     226           0 :     total_args += getPntrToArgument(i)->getNumberOfStoredValues();
     227             :   }
     228             : 
     229           0 :   if( inputdata.size()!=total_args ) {
     230           0 :     inputdata.resize( total_args );
     231             :   }
     232             : 
     233             :   total_args = 0;
     234           0 :   for(unsigned i=0; i<nargs; ++i) {
     235             :     Value* myarg = getPntrToArgument(i);
     236           0 :     total_args+= myarg->assignValues(View<float> {inputdata.data()+total_args,inputdata.size()-total_args});
     237             :   }
     238           0 : }
     239             : 
     240      138013 : void ActionWithVector::copyArgumentBookeeping( ArgumentsBookkeeping& argumentsMap ) const {
     241      138013 :   argumentsMap.nargs = getNumberOfArguments();
     242      138013 :   argumentsMap.ranks.resize( argumentsMap.nargs );
     243      138013 :   argumentsMap.shapestarts.resize( argumentsMap.nargs );
     244      138013 :   argumentsMap.argstarts.resize( argumentsMap.nargs );
     245             :   std::size_t s = 0;
     246             :   std::size_t ts = 0;
     247      391049 :   for(unsigned i=0; i<argumentsMap.nargs; ++i) {
     248             :     Value* myarg = getPntrToArgument(i);
     249      253036 :     argumentsMap.shapestarts[i] = ts;
     250      253036 :     argumentsMap.ranks[i] = myarg->getRank();
     251      253036 :     ts += argumentsMap.ranks[i];
     252      253036 :     argumentsMap.argstarts[i] = s;
     253      253036 :     s += myarg->getNumberOfStoredValues();
     254             :   }
     255      138013 :   argumentsMap.shapedata.resize( ts );
     256             :   ts = 0;
     257      138013 :   argumentsMap.ncols.resize( argumentsMap.nargs );
     258      138013 :   argumentsMap.bookstarts.resize( argumentsMap.nargs );
     259      138013 :   argumentsMap.booksizes.resize( argumentsMap.nargs );
     260             :   // Matrices whose bookkeeping was copied from the same value (e.g. the
     261             :   // components of a function of a matrix) share one copy of it: packing and
     262             :   // transferring (and, in the GPU task manager, comparing) 26 copies of
     263             :   // the bookkeeping of a contact matrix was a visible part of a step
     264      138013 :   std::vector<int> sharedWith( argumentsMap.nargs, -1 );
     265             :   std::size_t nbook = 0;
     266      391049 :   for(unsigned i=0; i<argumentsMap.nargs; ++i) {
     267             :     Value* myarg = getPntrToArgument(i);
     268      484533 :     for(unsigned j=0; j<argumentsMap.ranks[i]; ++j) {
     269      231497 :       argumentsMap.shapedata[ts] = myarg->getShape()[j];
     270      231497 :       ++ts;
     271             :     }
     272      253036 :     if( argumentsMap.ranks[i]==1 ) {
     273      161528 :       argumentsMap.ncols[i] = 1;
     274      161528 :       argumentsMap.booksizes[i] = 2*myarg->getShape()[0];
     275       91508 :     } else if( argumentsMap.ranks[i]==2 ) {
     276       34806 :       argumentsMap.ncols[i] = myarg->getNumberOfColumns();
     277       34806 :       argumentsMap.booksizes[i] = myarg->matrix_bookeeping.size();
     278       46855 :       for(unsigned k=0; k<i; ++k) {
     279       15868 :         if( argumentsMap.ranks[k]==2 && argumentsMap.booksizes[k]==argumentsMap.booksizes[i]
     280       21327 :             && getPntrToArgument(k)->getBookeepingRoot()==myarg->getBookeepingRoot() ) {
     281        3819 :           sharedWith[i] = k;
     282        3819 :           break;
     283             :         }
     284             :       }
     285             :     }
     286      253036 :     if( sharedWith[i]>=0 ) {
     287        3819 :       argumentsMap.bookstarts[i] = argumentsMap.bookstarts[sharedWith[i]];
     288             :     } else {
     289      249217 :       argumentsMap.bookstarts[i] = nbook;
     290      249217 :       nbook += argumentsMap.booksizes[i];
     291             :     }
     292             :   }
     293      138013 :   argumentsMap.bookeeping.resize( nbook );
     294      391049 :   for(unsigned i=0; i<argumentsMap.nargs; ++i) {
     295      253036 :     if( sharedWith[i]>=0 ) {
     296        3819 :       continue;
     297             :     }
     298             :     Value* myarg = getPntrToArgument(i);
     299      249217 :     std::size_t* dst = argumentsMap.bookeeping.data() + argumentsMap.bookstarts[i];
     300      249217 :     if( argumentsMap.ranks[i]==1 ) {
     301    32793443 :       for(std::size_t j=0; j<myarg->getShape()[0]; ++j) {
     302    32631915 :         dst[2*j] = 1;
     303    32631915 :         dst[2*j+1] = 0;
     304             :       }
     305       87689 :     } else if( argumentsMap.ranks[i]==2 ) {
     306             :       // the bookkeeping of the value is stored as unsigned, this one as size_t
     307             :       std::copy( myarg->matrix_bookeeping.begin(), myarg->matrix_bookeeping.end(), dst );
     308             :     }
     309             :   }
     310      138013 : }
     311             : 
     312      182306 : void ActionWithVector::getInputData( std::vector<double>& inputdata, ArgumentsBookkeeping& argumentsMap ) const {
     313      182306 :   getAtomicInputData( inputdata );
     314      182306 :   if( getNumberOfArguments()==0 ) {
     315             :     return;
     316             :   }
     317             :   // The masks are not part of the input data.  If the only arguments are masks
     318             :   // the buffer holds the atomic data only and must not be touched: the default
     319             :   // getArgumentInputData resizes it to the size of the (non-mask) arguments,
     320             :   // which would discard the atomic data (and leave a null pointer on the GPU)
     321      137780 :   if( nmask<=0 || getNumberOfArguments()>static_cast<unsigned>(nmask) ) {
     322      137315 :     getArgumentInputData( inputdata );
     323             :   }
     324      137780 :   copyArgumentBookeeping( argumentsMap );
     325             : }
     326             : 
     327           0 : void ActionWithVector::getInputData( std::vector<float>& inputdata, ArgumentsBookkeeping& argumentsMap ) const {
     328           0 :   getAtomicInputData( inputdata );
     329           0 :   if( getNumberOfArguments()==0 ) {
     330             :     return;
     331             :   }
     332             :   // The masks are not part of the input data.  If the only arguments are masks
     333             :   // the buffer holds the atomic data only and must not be touched: the default
     334             :   // getArgumentInputData resizes it to the size of the (non-mask) arguments,
     335             :   // which would discard the atomic data (and leave a null pointer on the GPU)
     336           0 :   if( nmask<=0 || getNumberOfArguments()>static_cast<unsigned>(nmask) ) {
     337           0 :     getArgumentInputData( inputdata );
     338             :   }
     339           0 :   copyArgumentBookeeping( argumentsMap );
     340             : }
     341             : 
     342             : // The stash holds the values of a task contiguously ([task][component]),
     343             : // the Values hold each component contiguously. The copies below write the
     344             : // storage of the Values directly (ActionWithVector is a friend of Value):
     345             : // Value::set per element, with its checks and periodicity test, was the
     346             : // most expensive part of a step for actions with many components (the 26
     347             : // components of the spherical harmonics on a million matrix elements).
     348             : // The components are independent, so the loop over them is parallel.
     349             : template <typename T>
     350      159235 : void ActionWithVector::stashToValuesImpl( const std::vector<unsigned>& partialTaskList, const std::vector<T>& stash ) {
     351             :   const std::size_t ntask = partialTaskList.size();
     352             :   const std::size_t ncomponents = getNumberOfComponents();
     353             : #ifdef _OPENMP
     354             :   // without OpenMP nt would be unused (an error with -Werror)
     355      159235 :   const unsigned nt = OpenMP::getGoodNumThreads( stash );
     356             : #endif
     357      159235 :   #pragma omp parallel for num_threads(nt) schedule(static)
     358             :   for(std::size_t i=0; i<ncomponents; ++i) {
     359             :     Value* myval = copyOutput(i);
     360             :     if( myval->hasDeriv || myval->getRank()==0 ) {
     361             :       for(std::size_t j=0; j<ntask; ++j) {
     362             :         myval->set( partialTaskList[j], stash[partialTaskList[j]*ncomponents+i] );
     363             :       }
     364             :       continue;
     365             :     }
     366             :     myval->value_set=true;
     367             :     double* d = myval->data.data();
     368             :     plumed_dbg_massert( ntask==0 || myval->data.size()>partialTaskList[ntask-1], "failing in " + myval->getName() );
     369             :     if( myval->isPeriodic() ) {
     370             :       for(std::size_t j=0; j<ntask; ++j) {
     371             :         const std::size_t t = partialTaskList[j];
     372             :         d[t] = stash[t*ncomponents+i];
     373             :         myval->applyPeriodicity(t);
     374             :       }
     375             :     } else {
     376             :       for(std::size_t j=0; j<ntask; ++j) {
     377             :         const std::size_t t = partialTaskList[j];
     378             :         d[t] = stash[t*ncomponents+i];
     379             :       }
     380             :     }
     381             :   }
     382      159235 : }
     383             : 
     384             : template <typename T>
     385      134239 : void ActionWithVector::forcesToStashImpl( const std::vector<unsigned>& partialTaskList, std::vector<T>& stash ) const {
     386             :   const std::size_t ntask = partialTaskList.size();
     387             :   const std::size_t ncomponents = getNumberOfComponents();
     388             : #ifdef _OPENMP
     389             :   // without OpenMP nt would be unused (an error with -Werror)
     390      134239 :   const unsigned nt = OpenMP::getGoodNumThreads( stash );
     391             : #endif
     392      134239 :   #pragma omp parallel for num_threads(nt) schedule(static)
     393             :   for(std::size_t i=0; i<ncomponents; ++i) {
     394             :     const double* f = getConstPntrToComponent(i)->inputForce.data();
     395             :     plumed_dbg_massert( ntask==0 || getConstPntrToComponent(i)->inputForce.size()>partialTaskList[ntask-1], "failing in " + getConstPntrToComponent(i)->getName() );
     396             :     for(std::size_t j=0; j<ntask; ++j) {
     397             :       const std::size_t t = partialTaskList[j];
     398             :       stash[t*ncomponents+i] = f[t];
     399             :     }
     400             :   }
     401      134239 : }
     402             : 
     403      159235 : void ActionWithVector::transferStashToValues( const std::vector<unsigned>& partialTaskList, const std::vector<double>& stash ) {
     404      159235 :   stashToValuesImpl( partialTaskList, stash );
     405      159235 : }
     406             : 
     407           0 : void ActionWithVector::transferStashToValues( const std::vector<unsigned>& partialTaskList, const std::vector<float>& stash ) {
     408           0 :   stashToValuesImpl( partialTaskList, stash );
     409           0 : }
     410             : 
     411      134239 : void ActionWithVector::transferForcesToStash( const std::vector<unsigned>& partialTaskList, std::vector<double>& stash ) const {
     412      134239 :   forcesToStashImpl( partialTaskList, stash );
     413      134239 : }
     414             : 
     415           0 : void ActionWithVector::transferForcesToStash( const std::vector<unsigned>& partialTaskList, std::vector<float>& stash ) const {
     416           0 :   forcesToStashImpl( partialTaskList, stash );
     417           0 : }
     418             : 
     419      211361 : void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) {
     420      211361 :   if( ntasks==0 ) {
     421      211361 :     if( getNumberOfArguments()==1 && getNumberOfComponents()==1 && getPntrToComponent(0)->getRank()==0 ) {
     422           0 :       if( !getPntrToArgument(0)->hasDerivatives() && getPntrToArgument(0)->getRank()==2 ) {
     423           0 :         ntasks = getPntrToArgument(0)->getShape()[0];
     424             :       } else {
     425           0 :         ntasks = getPntrToArgument(0)->getNumberOfValues();
     426             :       }
     427             :     } else {
     428      211361 :       plumed_assert( getNumberOfComponents()>0 && getPntrToComponent(0)->getRank()>0 );
     429      211361 :       if( getPntrToComponent(0)->hasDerivatives() ) {
     430       10085 :         ntasks = getPntrToComponent(0)->getNumberOfValues();
     431             :       } else {
     432      201276 :         ntasks = getPntrToComponent(0)->getShape()[0];
     433             :       }
     434             :     }
     435             :   }
     436      459776 :   for(unsigned i=0; i<getNumberOfComponents(); ++i) {
     437      248415 :     if( getPntrToComponent(i)->getRank()==0 ) {
     438           0 :       if( getNumberOfArguments()!=1 ) {
     439           0 :         error("mismatched numbers of tasks in streamed quantities");
     440             :       }
     441           0 :       if( getPntrToArgument(0)->hasDerivatives() && ntasks!=getPntrToArgument(0)->getNumberOfValues() ) {
     442           0 :         error("mismatched numbers of tasks in streamed quantities");
     443           0 :       } else if ( !getPntrToArgument(0)->hasDerivatives() && ntasks!=getPntrToArgument(0)->getShape()[0] ) {
     444           0 :         error("mismatched numbers of tasks in streamed quantities");
     445             :       }
     446      248415 :     } else if( getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getNumberOfValues() ) {
     447           0 :       error("mismatched numbers of tasks in streamed quantities");
     448      248415 :     } else if( !getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getShape()[0] ) {
     449           0 :       error("mismatched numbers of tasks in streamed quantities");
     450             :     }
     451             :   }
     452      211361 : }
     453             : 
     454      227117 : unsigned ActionWithVector::getNumberOfForceDerivatives() const {
     455             :   unsigned nforces=0;
     456      227117 :   unsigned nargs = getNumberOfArguments();
     457      227117 :   unsigned  nmasks = getNumberOfMasks();
     458      227117 :   if( nargs>=nmasks && nmasks>0 ) {
     459        2379 :     nargs = nargs - nmasks;
     460             :   }
     461      227117 :   if( getNumberOfAtoms()>0 ) {
     462       60501 :     nforces += 3*getNumberOfAtoms() + 9;
     463             :   }
     464      545249 :   for(unsigned i=0; i<nargs; ++i) {
     465      318132 :     nforces += getPntrToArgument(i)->getNumberOfStoredValues();
     466             :   }
     467      227117 :   return nforces;
     468             : }
     469             : 
     470      188388 : bool ActionWithVector::checkForForces() {
     471      188388 :   if( getPntrToComponent(0)->getRank()==0 ) {
     472        8676 :     return ActionWithValue::checkForForces();
     473             :   }
     474             : 
     475             :   // Check if there are any forces
     476             :   bool hasforce=false;
     477      227915 :   for(unsigned i=0; i<getNumberOfComponents(); ++i) {
     478      195496 :     if( getConstPntrToComponent(i)->getRank()>0 && getConstPntrToComponent(i)->forcesWereAdded() ) {
     479             :       hasforce=true;
     480             :       break;
     481             :     }
     482             :   }
     483      179712 :   if( !hasforce ) {
     484             :     return false;
     485             :   }
     486      147293 :   applyNonZeroRankForces( forcesForApply );
     487      147293 :   return true;
     488             : }
     489             : 
     490      188388 : void ActionWithVector::apply() {
     491      188388 :   unsigned nf =  getNumberOfForceDerivatives();
     492      188388 :   if( forcesForApply.size()!=nf ) {
     493        3722 :     forcesForApply.resize( nf, 0 );
     494             :   }
     495             : 
     496      188388 :   if( !checkForForces() ) {
     497       36322 :     return;
     498             :   }
     499             :   // Find the top of the chain and add forces
     500      152066 :   unsigned ind=0;
     501      152066 :   addForcesOnArguments( 0, forcesForApply, ind );
     502      152066 :   setForcesOnAtoms( forcesForApply, ind );
     503             : }
     504             : 
     505             : }