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 "PlumedMain.h"
24 : #include "ActionSet.h"
25 : #include "tools/OpenMP.h"
26 : #include "tools/Communicator.h"
27 :
28 : namespace PLMD {
29 :
30 5614 : void ActionWithVector::registerKeywords( Keywords& keys ) {
31 5614 : Action::registerKeywords( keys );
32 5614 : ActionAtomistic::registerKeywords( keys );
33 5614 : ActionWithValue::registerKeywords( keys );
34 5614 : keys.remove("NUMERICAL_DERIVATIVES");
35 5614 : ActionWithArguments::registerKeywords( keys );
36 11228 : keys.addFlag("SERIAL",false,"do the calculation in serial. Do not parallelize");
37 5614 : }
38 :
39 5001 : ActionWithVector::ActionWithVector(const ActionOptions&ao):
40 : Action(ao),
41 : ActionAtomistic(ao),
42 : ActionWithValue(ao),
43 : ActionWithArguments(ao),
44 5001 : serial(false),
45 5001 : action_to_do_before(NULL),
46 5001 : action_to_do_after(NULL),
47 5001 : never_reduce_tasks(false),
48 5001 : reduce_tasks(false),
49 5001 : atomsWereRetrieved(false),
50 5001 : done_in_chain(false) {
51 10002 : if( keywords.exists("SERIAL") ) {
52 4956 : parseFlag("SERIAL",serial);
53 : }
54 5001 : warning("you will likely obtain superior performance for this action if you update to version 2.11 or higher");
55 5001 : }
56 :
57 5001 : ActionWithVector::~ActionWithVector() {
58 5001 : if( action_to_do_before ) {
59 2485 : action_to_do_before->action_to_do_after=NULL;
60 : }
61 5001 : }
62 :
63 288992 : void ActionWithVector::lockRequests() {
64 : ActionAtomistic::lockRequests();
65 : ActionWithArguments::lockRequests();
66 288992 : }
67 :
68 288992 : void ActionWithVector::unlockRequests() {
69 : ActionAtomistic::unlockRequests();
70 : ActionWithArguments::unlockRequests();
71 288992 : }
72 :
73 0 : void ActionWithVector::calculateNumericalDerivatives(ActionWithValue* av) {
74 0 : plumed_merror("cannot calculate numerical derivative for action " + getName() + " with label " + getLabel() );
75 : }
76 :
77 473345 : void ActionWithVector::clearDerivatives( const bool& force ) {
78 473345 : if( !force && actionInChain() ) {
79 : return;
80 : }
81 354357 : ActionWithValue::clearDerivatives();
82 354357 : if( action_to_do_after ) {
83 184956 : action_to_do_after->clearDerivatives( true );
84 : }
85 : }
86 :
87 473345 : void ActionWithVector::clearInputForces( const bool& force ) {
88 473345 : if( !force && actionInChain() ) {
89 : return;
90 : }
91 354357 : ActionWithValue::clearInputForces();
92 354357 : if( action_to_do_after ) {
93 184956 : action_to_do_after->clearInputForces( true );
94 : }
95 : }
96 :
97 46 : const ActionWithVector* ActionWithVector::getFirstActionInChain() const {
98 46 : if( !actionInChain() ) {
99 : return this;
100 : }
101 46 : return action_to_do_before->getFirstActionInChain();
102 : }
103 :
104 209474253 : ActionWithVector* ActionWithVector::getFirstActionInChain() {
105 209474253 : if( !actionInChain() ) {
106 : return this;
107 : }
108 198584362 : return action_to_do_before->getFirstActionInChain();
109 : }
110 :
111 472099 : void ActionWithVector::retrieveAtoms( const bool& force ) {
112 472099 : if( !force && actionInChain() || atomsWereRetrieved ) {
113 : return;
114 : }
115 352781 : ActionAtomistic::retrieveAtoms();
116 352781 : atomsWereRetrieved = !actionInChain();
117 352781 : if( action_to_do_after ) {
118 184942 : action_to_do_after->retrieveAtoms( true );
119 : }
120 : }
121 :
122 46 : bool ActionWithVector::hasStoredArguments() const {
123 46 : std::string headstr=getFirstActionInChain()->getLabel();
124 100 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
125 68 : if( !getPntrToArgument(i)->ignoreStoredValue(headstr) ) {
126 : return true;
127 : }
128 : }
129 : return false;
130 : }
131 :
132 346 : bool ActionWithVector::argumentDependsOn( const std::string& headstr, ActionWithVector* faction, Value* thearg ) {
133 830 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
134 490 : if( this!=faction && thearg==getPntrToArgument(i) ) {
135 : return true;
136 : }
137 487 : ActionWithVector* av = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
138 487 : if( av && (av->getFirstActionInChain())->getLabel()==headstr ) {
139 180 : if( av->argumentDependsOn( headstr, faction, thearg ) ) {
140 : return true;
141 : };
142 : }
143 : }
144 : return false;
145 : }
146 :
147 3452 : unsigned ActionWithVector::buildArgumentStore( const unsigned& argstart ) {
148 : // Don't use chains for grids
149 10729 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
150 7667 : if( getPntrToArgument(i)->isConstant() ) {
151 911 : continue;
152 : }
153 6756 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction());
154 6756 : if( !av || getPntrToArgument(i)->getRank()>0 && getPntrToArgument(i)->hasDerivatives() ) {
155 390 : done_in_chain=false;
156 390 : break;
157 : }
158 : }
159 3452 : if( done_in_chain ) {
160 : std::vector<std::string> alabels;
161 : std::vector<ActionWithVector*> f_actions;
162 8582 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
163 : bool found=false;
164 6163 : std::string mylab = (getPntrToArgument(i)->getPntrToAction())->getLabel();
165 12731 : for(unsigned j=0; j<alabels.size(); ++j) {
166 8717 : if( alabels[j]==mylab ) {
167 : found=true;
168 : break;
169 : }
170 : }
171 6163 : if( !found ) {
172 4014 : alabels.push_back( mylab );
173 : }
174 :
175 : // If this is calculated in setup we never need to add to chain
176 6163 : if( getPntrToArgument(i)->isConstant() ) {
177 : continue;
178 : }
179 : // Find the chain we need to add this to from the arguments
180 5744 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction());
181 5744 : plumed_assert( av );
182 : found=false;
183 5744 : ActionWithVector* myact = av->getFirstActionInChain();
184 5744 : if( getPntrToArgument(i)->getRank()==1 && getPntrToArgument(i)->storedata ) {
185 5601 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
186 5282 : ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getPntrToArgument(j)->getPntrToAction() );
187 5282 : if( !aarg || i==j ) {
188 344 : continue;
189 : }
190 14180 : for(unsigned k=0; k<aarg->getNumberOfArguments(); ++k) {
191 9261 : if( aarg->getPntrToArgument(k)==getPntrToArgument(i) ) {
192 19 : done_in_chain=false;
193 19 : return reallyBuildArgumentStore( argstart );
194 : }
195 : }
196 : }
197 : }
198 6053 : for(unsigned j=0; j<f_actions.size(); ++j) {
199 3287 : if( f_actions[j]==myact ) {
200 : found=true;
201 : break;
202 : }
203 : }
204 5725 : if( !found ) {
205 2766 : if( f_actions.size()>0 ) {
206 298 : if( f_actions[0]->checkForDependency(myact) ) {
207 13 : getPntrToArgument(i)->buildDataStore();
208 : }
209 298 : if( myact->checkForDependency(f_actions[0]) ) {
210 0 : error("cannot deal with arguments in this order. Try swapping argument order");
211 : }
212 : }
213 2766 : if( !getPntrToArgument(i)->storedata && getPntrToArgument(i)->getRank()>0 ) {
214 2391 : f_actions.push_back( myact );
215 : }
216 : }
217 : }
218 : // Now make sure that everything we need is in the chain
219 2419 : if( f_actions.size()>0 ) {
220 : // Check everything for later f_actions is done before f_actions[0]
221 2377 : for(unsigned i=1; i<f_actions.size(); ++i) {
222 84 : ActionWithArguments* aarg = dynamic_cast<ActionWithArguments*>( f_actions[i] );
223 84 : if( !aarg || aarg->getNumberOfArguments()==0 ) {
224 30 : continue;
225 : }
226 175 : for(unsigned j=0; j<aarg->getNumberOfArguments(); ++j) {
227 130 : if( (aarg->getPntrToArgument(j))->isConstant() ) {
228 9 : continue ;
229 : }
230 : bool found=false;
231 121 : std::string dep_argname = (aarg->getPntrToArgument(j))->getPntrToAction()->getLabel();
232 17436 : for(const auto & pp : plumed.getActionSet()) {
233 : Action* p(pp.get());
234 : // Check if this is the dependency
235 17436 : if( p->getLabel()==dep_argname ) {
236 : found=true;
237 : break;
238 : }
239 : // Check if this is the first of the arguments that will appear in this chain
240 17324 : else if( p->getLabel()==f_actions[0]->getLabel() ) {
241 : break;
242 : }
243 : }
244 121 : if( !found ) {
245 9 : done_in_chain=false;
246 : break;
247 : }
248 : }
249 : // Stop trying to add things in the chain if we cannot
250 54 : if( !done_in_chain ) {
251 9 : return reallyBuildArgumentStore( argstart );
252 : }
253 : }
254 2293 : std::vector<std::string> empty(1);
255 2293 : empty[0] = f_actions[0]->getLabel();
256 2368 : for(unsigned i=1; i<f_actions.size(); ++i) {
257 75 : f_actions[0]->addActionToChain( empty, f_actions[i] );
258 : }
259 2293 : }
260 : // Now add this argument to the chain
261 : bool added=false;
262 2459 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
263 : // Add this function to jobs to do in recursive loop in previous action
264 2459 : if( getPntrToArgument(i)->getRank()>0 && !getPntrToArgument(i)->isConstant() ) {
265 2459 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
266 2459 : if( av && av->addActionToChain( alabels, this ) ) {
267 : added=true;
268 : break;
269 : }
270 : }
271 : }
272 2410 : plumed_massert(added, "could not add action " + getLabel() + " to chain of any of its arguments");
273 : // And get the number of derivatives
274 2410 : ActionWithVector* head=getFirstActionInChain();
275 2410 : unsigned nder=0;
276 2410 : arg_deriv_starts.resize( getNumberOfArguments() );
277 8532 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
278 : // Ensures that we ignore the grid (first argument) if we are evaluating a function on a grid
279 6122 : if( i==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) {
280 1 : continue ;
281 : }
282 6121 : ActionWithVector* iaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction());
283 6121 : if( actionInChain() && !getPntrToArgument(i)->ignoreStoredValue(head->getLabel()) ) {
284 677 : arg_deriv_starts[i] = 0;
285 677 : head->getNumberOfStreamedDerivatives( arg_deriv_starts[i], getPntrToArgument(i) );
286 5444 : } else if( iaction && iaction->isInSubChain(nder) ) {
287 2183 : arg_deriv_starts[i] = nder;
288 : // Add the total number of derivatives that we have by this point in the chain to nder
289 : if( iaction ) {
290 2183 : nder=0;
291 2183 : head->getNumberOfStreamedDerivatives( nder, getPntrToArgument(i) );
292 : }
293 : } else {
294 : // Check if we have already found this action
295 : int k=-1;
296 3261 : if( iaction ) {
297 3261 : ActionWithVector* ider_action=iaction->getActionWithDerivatives( iaction );
298 3482 : for(unsigned j=0; j<i; ++j) {
299 1633 : if( j==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) {
300 1 : continue ;
301 : }
302 1632 : ActionWithVector* jaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(j)->getPntrToAction());
303 1632 : if( jaction->getActionWithDerivatives(jaction)==ider_action || jaction->checkForDependency(ider_action) ) {
304 1412 : k=j;
305 1412 : break;
306 : }
307 : }
308 3261 : if( k>=0 ) {
309 1412 : arg_deriv_starts[i] = arg_deriv_starts[k];
310 1412 : continue;
311 : }
312 : }
313 :
314 1849 : if( i>0 ) {
315 : // This is a fudge so that inputs like this work:
316 : // c: CONTACT_MATRIX ATOMS=1-100
317 : // d: MATRIX_PRODUCT ARG=mat1,mat2
318 : // e: CUSTOM ARG=c,d
319 : // f: MATRIX_PRODUCT ARG=mat3,mat4
320 : // g: CUSTOM ARG=c,f
321 : // See symfunc rt-nbonds-q6 for an example
322 : // In this example when we set arg_deriv_starts[1] for f in g nder=number of derivatives of c
323 : // mder is equal to the number of derivatives by the time you get to f minus the number of derivatives for c
324 136 : unsigned mder=0;
325 136 : ActionWithVector* jaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(i-1)->getPntrToAction());
326 136 : if( jaction->action_to_do_after && !(jaction->action_to_do_after)->getNumberOfStoredValues( getPntrToArgument(i-1), mder, i, getArguments() ) ) {
327 0 : mder=0;
328 : }
329 136 : if( mder>0 ) {
330 1 : nder = nder + mder;
331 : }
332 : }
333 :
334 1849 : arg_deriv_starts[i] = nder;
335 : // Add the total number of derivatives that we have by this point in the chain to nder
336 1849 : if( iaction ) {
337 1849 : nder=0;
338 1849 : if( (getPntrToArgument(i)->getPntrToAction())->getName().find("DIFFERENCE")!=std::string::npos ) {
339 0 : ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getPntrToArgument(i)->getPntrToAction() );
340 0 : plumed_assert( aarg && aarg->getNumberOfArguments()==2 );
341 0 : head->getNumberOfStreamedDerivatives( nder, aarg->getPntrToArgument(0) );
342 0 : nder += (aarg->getPntrToArgument(1))->getNumberOfValues();
343 : } else {
344 1849 : head->getNumberOfStreamedDerivatives( nder, getPntrToArgument(i) );
345 : }
346 : }
347 : }
348 : }
349 2410 : nder=0;
350 2410 : head->getNumberOfStreamedDerivatives( nder, NULL );
351 2410 : return nder;
352 2438 : }
353 1014 : return reallyBuildArgumentStore( argstart );
354 : }
355 :
356 1042 : unsigned ActionWithVector::reallyBuildArgumentStore( const unsigned& argstart ) {
357 2846 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
358 1804 : if( getPntrToArgument(i)->getRank()>0 ) {
359 1558 : getPntrToArgument(i)->buildDataStore();
360 : }
361 : }
362 : unsigned nder=0;
363 1042 : arg_deriv_starts.resize( getNumberOfArguments() );
364 2846 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
365 1804 : arg_deriv_starts[i] = nder;
366 1804 : nder += getPntrToArgument(i)->getNumberOfValues();
367 : }
368 1042 : return nder;
369 : }
370 :
371 17841 : ActionWithVector* ActionWithVector::getActionWithDerivatives( ActionWithVector* depaction ) {
372 17841 : if( depaction==this || depaction->checkForDependency(this) ) {
373 11247 : if( getNumberOfAtoms()>0 ) {
374 4893 : return this;
375 : }
376 6851 : std::string c=getFirstActionInChain()->getLabel();
377 42713 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
378 36359 : if( !getPntrToArgument(i)->ignoreStoredValue(c) && !getPntrToArgument(i)->isConstant() ) {
379 : return this;
380 : }
381 : }
382 : }
383 12948 : plumed_assert( action_to_do_before );
384 12948 : return action_to_do_before->getActionWithDerivatives(depaction);
385 : }
386 :
387 10495 : bool ActionWithVector::addActionToChain( const std::vector<std::string>& alabels, ActionWithVector* act ) {
388 10495 : if( action_to_do_after ) {
389 7961 : bool state=action_to_do_after->addActionToChain( alabels, act );
390 7961 : return state;
391 : }
392 : // Check action is not already in chain
393 : std::vector<std::string> mylabels;
394 2534 : getFirstActionInChain()->getAllActionLabelsInChain( mylabels );
395 19308 : for(unsigned i=0; i<mylabels.size(); ++i) {
396 16774 : if( act->getLabel()==mylabels[i] ) {
397 : return true;
398 : }
399 : }
400 :
401 : // Check that everything that is required has been calculated
402 6648 : for(unsigned i=0; i<alabels.size(); ++i) {
403 : bool found=false;
404 19801 : for(unsigned j=0; j<mylabels.size(); ++j) {
405 19069 : if( alabels[i]==mylabels[j] ) {
406 : found=true;
407 : break;
408 : }
409 : }
410 4163 : if( !found ) {
411 732 : ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( alabels[i] );
412 732 : plumed_massert( av, "could not cast " + alabels[i] );
413 : bool storingall=true;
414 1526 : for(int j=0; j<av->getNumberOfComponents(); ++j) {
415 794 : if( !(av->getPntrToComponent(j))->storedata ) {
416 : storingall=false;
417 : }
418 : }
419 732 : if( !storingall ) {
420 : return false;
421 : }
422 : }
423 : }
424 : // This checks that there is nothing that will cause problems in the chain
425 2485 : mylabels.resize(0);
426 2485 : getFirstActionInChain()->getAllActionLabelsInChain( mylabels );
427 18806 : for(unsigned i=0; i<mylabels.size(); ++i) {
428 16321 : ActionWithVector* av1=plumed.getActionSet().selectWithLabel<ActionWithVector*>( mylabels[i] );
429 153887 : for(unsigned j=0; j<i; ++j) {
430 137566 : ActionWithVector* av2=plumed.getActionSet().selectWithLabel<ActionWithVector*>( mylabels[j] );
431 137566 : if( !av1->canBeAfterInChain( av2 ) ) {
432 0 : error("must calculate " + mylabels[j] + " before " + mylabels[i] );
433 : }
434 : }
435 : }
436 2485 : action_to_do_after=act;
437 2485 : act->action_to_do_before=this;
438 2485 : updateTaskListReductionStatus();
439 2485 : ActionWithVector* head = getFirstActionInChain();
440 2485 : head->broadcastThatTasksAreReduced( head );
441 2485 : head->finishChainBuild( act );
442 : return true;
443 2534 : }
444 :
445 2640 : void ActionWithVector::updateTaskListReductionStatus() {
446 2640 : ActionWithVector* head = getFirstActionInChain();
447 : std::vector<ActionWithVector*> task_reducing_actions;
448 2640 : head->canReduceTasks( task_reducing_actions );
449 2640 : if( task_reducing_actions.size()>0 ) {
450 319 : head->reduce_tasks=true;
451 : }
452 2640 : }
453 :
454 3841714 : void ActionWithVector::broadcastThatTasksAreReduced( ActionWithVector* aselect ) {
455 3841714 : std::string c=getFirstActionInChain()->getLabel();
456 18080207 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
457 14238493 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
458 71083 : ActionWithVector* av = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
459 71083 : if( av ) {
460 : bool found=false;
461 12049 : ActionWithVector* av_head = av->getFirstActionInChain();
462 17754 : for(unsigned i=0; i<av_head->task_control_list.size(); ++i) {
463 17217 : if( aselect==av_head->task_control_list[i] ) {
464 : found=true;
465 : break;
466 : }
467 : }
468 12049 : if( !found ) {
469 537 : av_head->task_control_list.insert( av_head->task_control_list.begin(), aselect );
470 : }
471 :
472 12049 : av_head->reduce_tasks=true;
473 12049 : av_head->updateTaskReductionFlag( av_head->reduce_tasks );
474 : }
475 : }
476 : }
477 3841714 : if( action_to_do_after ) {
478 3609739 : action_to_do_after->broadcastThatTasksAreReduced( aselect );
479 : }
480 3841714 : }
481 :
482 229490 : void ActionWithVector::updateTaskReductionFlag( bool& head_reduce_tasks ) {
483 229490 : if( actionInChain() ) {
484 217441 : plumed_assert( task_control_list.size()==0 );
485 : } else {
486 30600 : for(unsigned i=0; i<task_control_list.size(); ++i) {
487 18551 : if( !(task_control_list[i]->getFirstActionInChain())->reduce_tasks ) {
488 17359 : head_reduce_tasks=false;
489 : }
490 : }
491 : }
492 229490 : broadcastThatTasksAreReduced( getFirstActionInChain() );
493 229490 : if( action_to_do_after ) {
494 217441 : action_to_do_after->updateTaskReductionFlag( head_reduce_tasks );
495 : }
496 229490 : }
497 :
498 19013 : void ActionWithVector::canReduceTasks( std::vector<ActionWithVector*>& task_reducing_actions ) {
499 19013 : areAllTasksRequired( task_reducing_actions );
500 19013 : if( action_to_do_after ) {
501 16373 : action_to_do_after->canReduceTasks( task_reducing_actions );
502 : }
503 19013 : }
504 :
505 3874 : void ActionWithVector::finishChainBuild( ActionWithVector* act ) {
506 3874 : if( action_to_do_after ) {
507 3187 : action_to_do_after->finishChainBuild( act );
508 : }
509 3874 : }
510 :
511 33175 : void ActionWithVector::getAllActionLabelsInChain( std::vector<std::string>& mylabels ) const {
512 : bool found = false ;
513 313793 : for(unsigned i=0; i<mylabels.size(); ++i) {
514 280618 : if( getLabel()==mylabels[i] ) {
515 : found=true;
516 : }
517 : }
518 33175 : if( !found ) {
519 33175 : mylabels.push_back( getLabel() );
520 : }
521 33175 : if( action_to_do_after ) {
522 28145 : action_to_do_after->getAllActionLabelsInChain( mylabels );
523 : }
524 33175 : }
525 :
526 6014383 : void ActionWithVector::taskIsActive( const unsigned& current, int& flag ) const {
527 6014383 : if( isActive() ) {
528 6014383 : flag = checkTaskStatus( current, flag );
529 : }
530 6014383 : if( flag<=0 && action_to_do_after ) {
531 3732005 : action_to_do_after->taskIsActive( current, flag );
532 : }
533 6014383 : }
534 :
535 1400 : void ActionWithVector::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
536 1445 : for(unsigned i=0; i<task_control_list.size(); ++i ) {
537 45 : task_control_list[i]->getAdditionalTasksRequired( action, atasks );
538 : }
539 1400 : }
540 :
541 289998 : void ActionWithVector::prepare() {
542 289998 : active_tasks.resize(0);
543 289998 : atomsWereRetrieved=false;
544 289998 : }
545 :
546 147420 : std::vector<unsigned>& ActionWithVector::getListOfActiveTasks( ActionWithVector* action ) {
547 147420 : if( active_tasks.size()>0 ) {
548 11971 : return active_tasks;
549 : }
550 135449 : unsigned ntasks=0;
551 135449 : getNumberOfTasks( ntasks );
552 :
553 135449 : unsigned stride=comm.Get_size();
554 135449 : unsigned rank=comm.Get_rank();
555 135449 : if(serial) {
556 : stride=1;
557 : rank=0;
558 : }
559 :
560 : // Get number of threads for OpenMP
561 135449 : unsigned nt=OpenMP::getNumThreads();
562 135449 : if( nt*stride*10>ntasks ) {
563 58134 : nt=ntasks/stride/10;
564 : }
565 : if( nt==0 ) {
566 : nt=1;
567 : }
568 :
569 135449 : if( !never_reduce_tasks && reduce_tasks ) {
570 1450 : if( task_control_list.size()>0 ) {
571 : // Get the list of tasks that are active in the action that uses the output of this action
572 220 : for(unsigned i=0; i<task_control_list.size(); ++i) {
573 121 : task_control_list[i]->retrieveAtoms();
574 121 : active_tasks = task_control_list[i]->getListOfActiveTasks( action );
575 : }
576 : // Now work out else we need from here to calculate the later action
577 99 : getAdditionalTasksRequired( action, active_tasks );
578 : } else {
579 1351 : std::vector<int> taskFlags( ntasks, -1 );
580 :
581 1351 : #pragma omp parallel num_threads(nt)
582 : {
583 : #pragma omp for nowait
584 : for(unsigned i=rank; i<ntasks; i+=stride ) {
585 : taskIsActive( i, taskFlags[i] );
586 : }
587 : }
588 2482527 : for(unsigned i=0; i<ntasks; ++i) {
589 2481176 : taskFlags[i] = std::abs( taskFlags[i] );
590 : }
591 1351 : if( !serial ) {
592 1351 : comm.Sum( taskFlags );
593 : }
594 :
595 : unsigned nt=0;
596 2482527 : for(unsigned i=0; i<ntasks; ++i) {
597 2481176 : if( taskFlags[i]>=stride ) {
598 261189 : nt++;
599 : }
600 : }
601 1351 : active_tasks.resize(nt);
602 : nt=0;
603 2482527 : for(unsigned i=0; i<ntasks; ++i) {
604 2481176 : if( taskFlags[i]>=stride ) {
605 261189 : active_tasks[nt]=i;
606 261189 : nt++;
607 : }
608 : }
609 1351 : getAdditionalTasksRequired( this, active_tasks );
610 : }
611 : } else {
612 133999 : active_tasks.resize( ntasks );
613 5870454 : for(unsigned i=0; i<ntasks; ++i) {
614 5736455 : active_tasks[i]=i;
615 : }
616 : }
617 135449 : return active_tasks;
618 : }
619 :
620 135604 : void ActionWithVector::runAllTasks() {
621 : // Skip this if this is done elsewhere
622 135604 : if( action_to_do_before ) {
623 35 : return;
624 : }
625 :
626 135569 : unsigned stride=comm.Get_size();
627 135569 : unsigned rank=comm.Get_rank();
628 135569 : if(serial) {
629 : stride=1;
630 : rank=0;
631 : }
632 :
633 : // Get the list of active tasks
634 135569 : std::vector<unsigned> & partialTaskList( getListOfActiveTasks( this ) );
635 135569 : unsigned nactive_tasks=partialTaskList.size();
636 :
637 : // Get number of threads for OpenMP
638 135569 : unsigned nt=OpenMP::getNumThreads();
639 135569 : if( nt*stride*10>nactive_tasks ) {
640 58325 : nt=nactive_tasks/stride/10;
641 : }
642 135569 : if( nt==0 ) {
643 20063 : nt=1;
644 : }
645 :
646 : // Now do all preparations required to run all the tasks
647 : // prepareForTaskLoop();
648 :
649 : // Get the total number of streamed quantities that we need
650 135569 : unsigned nquants=0, nmatrices=0, maxcol=0, nbooks=0;
651 135569 : getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol, nbooks );
652 : // Get size for buffer
653 135569 : unsigned bufsize=0;
654 135569 : getSizeOfBuffer( nactive_tasks, bufsize );
655 135569 : if( buffer.size()!=bufsize ) {
656 2347 : buffer.resize( bufsize );
657 : }
658 : // Clear buffer
659 135569 : buffer.assign( buffer.size(), 0.0 );
660 :
661 : // Recover the number of derivatives we require
662 135569 : unsigned nderivatives = 0;
663 135569 : bool gridsInStream=checkForGrids(nderivatives);
664 135569 : if( !doNotCalculateDerivatives() && !gridsInStream ) {
665 105507 : getNumberOfStreamedDerivatives( nderivatives, NULL );
666 : }
667 :
668 135569 : #pragma omp parallel num_threads(nt)
669 : {
670 : std::vector<double> omp_buffer;
671 : if( nt>1 ) {
672 : omp_buffer.resize( bufsize, 0.0 );
673 : }
674 : MultiValue myvals( nquants, nderivatives, nmatrices, maxcol, nbooks );
675 : myvals.clearAll();
676 :
677 : #pragma omp for nowait
678 : for(unsigned i=rank; i<nactive_tasks; i+=stride) {
679 : // Calculate the stuff in the loop for this action
680 : runTask( partialTaskList[i], myvals );
681 :
682 : // Now transfer the data to the actions that accumulate values from the calculated quantities
683 : if( nt>1 ) {
684 : gatherAccumulators( partialTaskList[i], myvals, omp_buffer );
685 : } else {
686 : gatherAccumulators( partialTaskList[i], myvals, buffer );
687 : }
688 :
689 : // Clear the value
690 : myvals.clearAll();
691 : }
692 : #pragma omp critical
693 : gatherThreads( nt, bufsize, omp_buffer, buffer, myvals );
694 : }
695 :
696 : // MPI Gather everything
697 135569 : if( !serial && buffer.size()>0 ) {
698 135562 : gatherProcesses( buffer );
699 : }
700 135569 : finishComputations( buffer );
701 : }
702 :
703 212813 : void ActionWithVector::gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) {
704 212813 : if( nt>1 )
705 98748328 : for(unsigned i=0; i<bufsize; ++i) {
706 98593840 : buffer[i]+=omp_buffer[i];
707 : }
708 212813 : }
709 :
710 135562 : void ActionWithVector::gatherProcesses( std::vector<double>& buffer ) {
711 135562 : comm.Sum( buffer );
712 135562 : }
713 :
714 320544 : bool ActionWithVector::checkForGrids( unsigned& nder ) const {
715 659569 : for(int i=0; i<getNumberOfComponents(); ++i) {
716 347786 : if( getConstPntrToComponent(i)->getRank()>0 && getConstPntrToComponent(i)->hasDerivatives() ) {
717 8761 : nder=getConstPntrToComponent(i)->getNumberOfGridDerivatives();
718 8761 : return true;
719 : }
720 : }
721 311783 : if( action_to_do_after ) {
722 184975 : return action_to_do_after->checkForGrids(nder);
723 : }
724 : return false;
725 : }
726 :
727 317763 : void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) {
728 317763 : if( ntasks==0 ) {
729 132784 : if( getNumberOfArguments()==1 && getNumberOfComponents()==1 && getPntrToComponent(0)->getRank()==0 ) {
730 2293 : if( !getPntrToArgument(0)->hasDerivatives() && getPntrToArgument(0)->getRank()==2 ) {
731 19 : ntasks = getPntrToArgument(0)->getShape()[0];
732 : } else {
733 2274 : ntasks = getPntrToArgument(0)->getNumberOfValues();
734 : }
735 : } else {
736 130491 : plumed_assert( getNumberOfComponents()>0 && getPntrToComponent(0)->getRank()>0 );
737 130491 : if( getPntrToComponent(0)->hasDerivatives() ) {
738 6077 : ntasks = getPntrToComponent(0)->getNumberOfValues();
739 : } else {
740 124414 : ntasks = getPntrToComponent(0)->getShape()[0];
741 : }
742 : }
743 : }
744 662764 : for(int i=0; i<getNumberOfComponents(); ++i) {
745 345001 : if( getPntrToComponent(i)->getRank()==0 ) {
746 66066 : if( getNumberOfArguments()!=1 ) {
747 0 : error("mismatched numbers of tasks in streamed quantities");
748 : }
749 66066 : if( getPntrToArgument(0)->hasDerivatives() && ntasks!=getPntrToArgument(0)->getNumberOfValues() ) {
750 0 : error("mismatched numbers of tasks in streamed quantities");
751 66066 : } else if ( !getPntrToArgument(0)->hasDerivatives() && ntasks!=getPntrToArgument(0)->getShape()[0] ) {
752 0 : error("mismatched numbers of tasks in streamed quantities");
753 : }
754 278935 : } else if( getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getNumberOfValues() ) {
755 0 : error("mismatched numbers of tasks in streamed quantities");
756 278935 : } else if( !getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getShape()[0] ) {
757 0 : error("mismatched numbers of tasks in streamed quantities");
758 : }
759 : }
760 317763 : if( action_to_do_after ) {
761 184979 : action_to_do_after->getNumberOfTasks( ntasks );
762 : }
763 317763 : }
764 :
765 424565 : void ActionWithVector::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
766 1081840 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
767 657275 : if( !getPntrToArgument(i)->ignoreStoredValue(headstr) ) {
768 415109 : getPntrToArgument(i)->streampos=nquants;
769 415109 : nquants++;
770 : }
771 : }
772 891463 : for(int i=0; i<getNumberOfComponents(); ++i) {
773 466898 : getPntrToComponent(i)->streampos=nquants;
774 466898 : nquants++;
775 : }
776 424565 : }
777 :
778 424565 : void ActionWithVector::getNumberOfStreamedQuantities( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
779 424565 : setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
780 424565 : if( action_to_do_after ) {
781 213997 : action_to_do_after->getNumberOfStreamedQuantities( headstr, nquants, nmat, maxcol, nbookeeping );
782 : }
783 424565 : }
784 :
785 320544 : void ActionWithVector::getSizeOfBuffer( const unsigned& nactive_tasks, unsigned& bufsize ) {
786 668330 : for(int i=0; i<getNumberOfComponents(); ++i) {
787 347786 : getPntrToComponent(i)->bufstart=bufsize;
788 347786 : bufsize += getPntrToComponent(i)->data.size();
789 : }
790 320544 : if( action_to_do_after ) {
791 184975 : action_to_do_after->getSizeOfBuffer( nactive_tasks, bufsize );
792 : }
793 320544 : }
794 :
795 422811 : void ActionWithVector::getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ) {
796 422811 : std::string c=getFirstActionInChain()->getLabel();
797 1140863 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
798 718729 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
799 385417 : if( getPntrToArgument(i)==stopat ) {
800 : return;
801 : }
802 384740 : nderivatives += getPntrToArgument(i)->getNumberOfValues();
803 : }
804 : }
805 422134 : if( getNumberOfAtoms()>0 ) {
806 53871 : nderivatives += 3*getNumberOfAtoms() + 9;
807 : }
808 : // Don't do the whole chain if we have been told to stop early
809 422134 : if( stopat && stopat->getPntrToAction()==this ) {
810 : return;
811 : }
812 :
813 417928 : if( action_to_do_after ) {
814 235012 : action_to_do_after->getNumberOfStreamedDerivatives( nderivatives, stopat );
815 : }
816 : }
817 :
818 136 : bool ActionWithVector::getNumberOfStoredValues( Value* startat, unsigned& nvals, const unsigned& astart, const std::vector<Value*>& stopat ) {
819 181 : for(unsigned j=astart; j<stopat.size(); ++j) {
820 151 : if( stopat[j] && (stopat[j]->getPntrToAction()==this || (stopat[j]->getPntrToAction())->checkForDependency(this)) ) {
821 106 : return true;
822 : }
823 : }
824 :
825 30 : std::string c=getFirstActionInChain()->getLabel();
826 77 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
827 47 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
828 4 : for(unsigned j=astart; j<stopat.size(); ++j) {
829 2 : if( getPntrToArgument(i)==stopat[j] ) {
830 : return true;
831 : }
832 : }
833 2 : nvals += getPntrToArgument(i)->getNumberOfValues();
834 : }
835 : }
836 30 : if( startat->getPntrToAction()!=this && getNumberOfAtoms()>0 ) {
837 : return false;
838 : }
839 :
840 30 : if( action_to_do_after ) {
841 30 : return action_to_do_after->getNumberOfStoredValues( startat, nvals, astart, stopat );
842 : }
843 : return false;
844 : }
845 :
846 13555748 : void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) const {
847 13555748 : if( isActive() ) {
848 12021026 : myvals.setTaskIndex(current);
849 12021026 : myvals.vector_call=true;
850 12021026 : performTask( current, myvals );
851 : }
852 13555748 : if( action_to_do_after ) {
853 6301873 : action_to_do_after->runTask( current, myvals );
854 : }
855 13555748 : }
856 :
857 11442823 : void ActionWithVector::gatherAccumulators( const unsigned& taskCode, const MultiValue& myvals, std::vector<double>& buffer ) const {
858 11442823 : if( isActive() ) {
859 26584722 : for(int i=0; i<getNumberOfComponents(); ++i) {
860 16520398 : unsigned bufstart = getConstPntrToComponent(i)->bufstart;
861 : // This looks after storing of scalars that are summed from vectors/matrices
862 16520398 : if( getConstPntrToComponent(i)->getRank()==0 ) {
863 : plumed_dbg_massert( bufstart<buffer.size(), "problem in " + getLabel() );
864 1336426 : unsigned sind = getConstPntrToComponent(i)->streampos;
865 1336426 : buffer[bufstart] += myvals.get(sind);
866 1336426 : if( getConstPntrToComponent(i)->hasDerivatives() ) {
867 26205534 : for(unsigned k=0; k<myvals.getNumberActive(sind); ++k) {
868 24869108 : unsigned kindex = myvals.getActiveIndex(sind,k);
869 : plumed_dbg_massert( bufstart+1+kindex<buffer.size(), "problem in " + getLabel() );
870 24869108 : buffer[bufstart + 1 + kindex] += myvals.getDerivative(sind,kindex);
871 : }
872 : }
873 : // This looks after storing of vectors
874 15183972 : } else if( getConstPntrToComponent(i)->storedata ) {
875 5738668 : gatherStoredValue( i, taskCode, myvals, bufstart, buffer );
876 : }
877 : }
878 : }
879 11442823 : if( action_to_do_after ) {
880 5614243 : action_to_do_after->gatherAccumulators( taskCode, myvals, buffer );
881 : }
882 11442823 : }
883 :
884 4071938 : void ActionWithVector::gatherStoredValue( const unsigned& valindex, const unsigned& taskCode, const MultiValue& myvals,
885 : const unsigned& bufstart, std::vector<double>& buffer ) const {
886 : plumed_dbg_assert( getConstPntrToComponent(valindex)->getRank()==1 && !getConstPntrToComponent(valindex)->hasDeriv );
887 4071938 : unsigned vindex = getConstPntrToComponent(valindex)->bufstart + taskCode;
888 : plumed_dbg_massert( vindex<buffer.size(), "failing in " + getLabel() );
889 4071938 : buffer[vindex] += myvals.get(getConstPntrToComponent(valindex)->streampos);
890 4071938 : }
891 :
892 320544 : void ActionWithVector::finishComputations( const std::vector<double>& buf ) {
893 320544 : if( isActive() ) {
894 536394 : for(int i=0; i<getNumberOfComponents(); ++i) {
895 : // This gathers vectors and grids at the end of the calculation
896 281818 : unsigned bufstart = getPntrToComponent(i)->bufstart;
897 281818 : getPntrToComponent(i)->data.assign( getPntrToComponent(i)->data.size(), 0 );
898 281818 : if( (getPntrToComponent(i)->getRank()>0 && getPntrToComponent(i)->hasDerivatives()) || getPntrToComponent(i)->storedata ) {
899 136650 : unsigned sz_v = getPntrToComponent(i)->data.size();
900 47392382 : for(unsigned j=0; j<sz_v; ++j) {
901 : plumed_dbg_assert( bufstart+j<buf.size() );
902 47255732 : getPntrToComponent(i)->add( j, buf[bufstart+j] );
903 : }
904 : // Make sure single values are set
905 145168 : } else if( getPntrToComponent(i)->getRank()==0 ) {
906 43628 : getPntrToComponent(i)->set( buf[bufstart] );
907 : }
908 : // This gathers derivatives of scalars
909 281818 : if( !doNotCalculateDerivatives() && getPntrToComponent(i)->hasDeriv && getPntrToComponent(i)->getRank()==0 ) {
910 10216614 : for(unsigned j=0; j<getPntrToComponent(i)->getNumberOfDerivatives(); ++j) {
911 10180276 : getPntrToComponent(i)->setDerivative( j, buf[bufstart+1+j] );
912 : }
913 : }
914 : }
915 : }
916 320544 : if( action_to_do_after ) {
917 184975 : action_to_do_after->finishComputations( buf );
918 : }
919 320544 : }
920 :
921 294478 : bool ActionWithVector::checkChainForNonScalarForces() const {
922 528408 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
923 308929 : if( getConstPntrToComponent(i)->getRank()>0 && getConstPntrToComponent(i)->forcesWereAdded() ) {
924 : return true;
925 : }
926 : }
927 219479 : if( action_to_do_after ) {
928 167738 : return action_to_do_after->checkChainForNonScalarForces();
929 : }
930 : return false;
931 : }
932 :
933 257303 : bool ActionWithVector::checkForForces() {
934 257303 : if( getPntrToComponent(0)->getRank()==0 ) {
935 53033 : return ActionWithValue::checkForForces();
936 204270 : } else if( actionInChain() ) {
937 : return false;
938 : }
939 :
940 : // Check if there are any forces
941 126740 : if( !checkChainForNonScalarForces() ) {
942 : return false;
943 : }
944 :
945 : // Setup MPI parallel loop
946 74999 : unsigned stride=comm.Get_size();
947 74999 : unsigned rank=comm.Get_rank();
948 74999 : if(serial) {
949 : stride=1;
950 : rank=0;
951 : }
952 :
953 : // Get the number of tasks
954 : std::vector<unsigned> force_tasks;
955 74999 : getForceTasks( force_tasks );
956 74999 : Tools::removeDuplicates(force_tasks);
957 74999 : unsigned nf_tasks=force_tasks.size();
958 :
959 : // Get number of threads for OpenMP
960 74999 : unsigned nt=OpenMP::getNumThreads();
961 74999 : if( nt*stride*10>nf_tasks ) {
962 42369 : nt=nf_tasks/stride/10;
963 : }
964 : if( nt==0 ) {
965 : nt=1;
966 : }
967 :
968 : // Now determine how big the multivalue needs to be
969 74999 : unsigned nquants=0, nmatrices=0, maxcol=0, nbooks=0;
970 74999 : getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol, nbooks );
971 : // Recover the number of derivatives we require (this should be equal to the number of forces)
972 74999 : unsigned nderiv=0;
973 74999 : getNumberOfStreamedDerivatives( nderiv, NULL );
974 74999 : if( forcesForApply.size()!=nderiv ) {
975 381 : forcesForApply.resize( nderiv );
976 : }
977 : // Clear force buffer
978 74999 : forcesForApply.assign( forcesForApply.size(), 0.0 );
979 :
980 74999 : #pragma omp parallel num_threads(nt)
981 : {
982 : std::vector<double> omp_forces;
983 : if( nt>1 ) {
984 : omp_forces.resize( forcesForApply.size(), 0.0 );
985 : }
986 : MultiValue myvals( nquants, nderiv, nmatrices, maxcol, nbooks );
987 : myvals.clearAll();
988 :
989 : #pragma omp for nowait
990 : for(unsigned i=rank; i<nf_tasks; i+=stride) {
991 : runTask( force_tasks[i], myvals );
992 :
993 : // Now get the forces
994 : if( nt>1 ) {
995 : gatherForces( force_tasks[i], myvals, omp_forces );
996 : } else {
997 : gatherForces( force_tasks[i], myvals, forcesForApply );
998 : }
999 :
1000 : myvals.clearAll();
1001 : }
1002 : #pragma omp critical
1003 : if(nt>1)
1004 : for(unsigned i=0; i<forcesForApply.size(); ++i) {
1005 : forcesForApply[i]+=omp_forces[i];
1006 : }
1007 : }
1008 : // MPI Gather on forces
1009 74999 : if( !serial ) {
1010 74999 : comm.Sum( forcesForApply );
1011 : }
1012 : return true;
1013 : }
1014 :
1015 2050310 : bool ActionWithVector::checkComponentsForForce() const {
1016 2793497 : for(unsigned i=0; i<values.size(); ++i) {
1017 2557587 : if( getConstPntrToComponent(i)->forcesWereAdded() ) {
1018 : return true;
1019 : }
1020 : }
1021 : return false;
1022 : }
1023 :
1024 5921624 : bool ActionWithVector::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
1025 : plumed_dbg_assert( !(myval->getRank()==2 && !myval->hasDerivatives()) );
1026 5921624 : return fabs(myval->getForce(itask))>epsilon;
1027 : }
1028 :
1029 90929 : void ActionWithVector::updateForceTasksFromValue( const Value* myval, std::vector<unsigned>& force_tasks ) const {
1030 90929 : if( myval->getRank()>0 && myval->forcesWereAdded() ) {
1031 85640 : unsigned nt = myval->getNumberOfValues();
1032 85640 : if( !myval->hasDerivatives() ) {
1033 85171 : nt = myval->getShape()[0];
1034 : }
1035 6510085 : for(unsigned i=0; i<nt; ++i) {
1036 6424445 : if( checkForTaskForce(i, myval) ) {
1037 2079412 : force_tasks.push_back( i );
1038 : }
1039 : }
1040 : }
1041 90929 : }
1042 :
1043 104021 : void ActionWithVector::getForceTasks( std::vector<unsigned>& force_tasks ) const {
1044 104021 : if( isActive() && checkComponentsForForce() ) {
1045 171997 : for(unsigned k=0; k<values.size(); ++k) {
1046 91539 : updateForceTasksFromValue( getConstPntrToComponent(k), force_tasks );
1047 : }
1048 : }
1049 104021 : if( action_to_do_after ) {
1050 29022 : action_to_do_after->getForceTasks( force_tasks );
1051 : }
1052 104021 : }
1053 :
1054 1570301 : void ActionWithVector::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
1055 : plumed_dbg_assert( myval->storedata && !(myval->getRank()==2 && !myval->hasDerivatives()) );
1056 1570301 : double fforce = myval->getForce(itask);
1057 : unsigned sspos = myval->getPositionInStream();
1058 22066949 : for(unsigned j=0; j<myvals.getNumberActive(sspos); ++j) {
1059 20496648 : unsigned jder=myvals.getActiveIndex(sspos, j);
1060 : plumed_dbg_assert( jder<forces.size() );
1061 20496648 : forces[jder] += fforce*myvals.getDerivative( sspos, jder );
1062 : }
1063 1570301 : }
1064 :
1065 2112925 : void ActionWithVector::gatherForces( const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
1066 2112925 : if( isActive() && checkComponentsForForce() ) {
1067 3599634 : for(unsigned k=0; k<getNumberOfComponents(); ++k) {
1068 1865692 : const Value* myval=getConstPntrToComponent(k);
1069 1865692 : if( myval->getRank()>0 && myval->forcesWereAdded() ) {
1070 1699333 : gatherForcesOnStoredValue( myval, itask, myvals, forces );
1071 : }
1072 : }
1073 : }
1074 2112925 : if( action_to_do_after ) {
1075 687630 : action_to_do_after->gatherForces( itask, myvals, forces );
1076 : }
1077 2112925 : }
1078 :
1079 257303 : void ActionWithVector::apply() {
1080 257303 : if( !checkForForces() ) {
1081 139205 : return;
1082 : }
1083 : // Find the top of the chain and add forces
1084 118098 : unsigned ind=0;
1085 118098 : getFirstActionInChain()->addForcesToInput( getForcesToApply(), ind );
1086 : }
1087 :
1088 185112 : void ActionWithVector::addForcesToInput( const std::vector<double>& forcesToApply, unsigned& ind ) {
1089 185112 : if( ind>=forcesToApply.size() ) {
1090 : return;
1091 : }
1092 135354 : addForcesOnArguments( 0, forcesToApply, ind, getFirstActionInChain()->getLabel() );
1093 135354 : setForcesOnAtoms( forcesToApply, ind );
1094 135354 : if( action_to_do_after ) {
1095 67014 : action_to_do_after->addForcesToInput( forcesToApply, ind );
1096 : }
1097 : }
1098 :
1099 : }
|
