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 5426 : void ActionWithVector::registerKeywords( Keywords& keys ) {
31 5426 : Action::registerKeywords( keys );
32 5426 : ActionAtomistic::registerKeywords( keys );
33 5426 : ActionWithValue::registerKeywords( keys );
34 5426 : keys.remove("NUMERICAL_DERIVATIVES");
35 5426 : ActionWithArguments::registerKeywords( keys );
36 5426 : keys.addFlag("SERIAL",false,"do the calculation in serial. Do not parallelize");
37 5426 : }
38 :
39 4937 : ActionWithVector::ActionWithVector(const ActionOptions&ao):
40 : Action(ao),
41 : ActionAtomistic(ao),
42 : ActionWithValue(ao),
43 : ActionWithArguments(ao),
44 4937 : serial(false),
45 4937 : action_to_do_before(NULL),
46 4937 : action_to_do_after(NULL),
47 4937 : never_reduce_tasks(false),
48 4937 : reduce_tasks(false),
49 4937 : atomsWereRetrieved(false),
50 4937 : done_in_chain(false) {
51 4937 : if( keywords.exists("SERIAL") ) {
52 4892 : parseFlag("SERIAL",serial);
53 : }
54 4937 : }
55 :
56 4937 : ActionWithVector::~ActionWithVector() {
57 4937 : if( action_to_do_before ) {
58 2538 : action_to_do_before->action_to_do_after=NULL;
59 : }
60 4937 : }
61 :
62 273650 : void ActionWithVector::lockRequests() {
63 : ActionAtomistic::lockRequests();
64 : ActionWithArguments::lockRequests();
65 273650 : }
66 :
67 273650 : void ActionWithVector::unlockRequests() {
68 : ActionAtomistic::unlockRequests();
69 : ActionWithArguments::unlockRequests();
70 273650 : }
71 :
72 0 : void ActionWithVector::calculateNumericalDerivatives(ActionWithValue* av) {
73 0 : plumed_merror("cannot calculate numerical derivative for action " + getName() + " with label " + getLabel() );
74 : }
75 :
76 457848 : void ActionWithVector::clearDerivatives( const bool& force ) {
77 457848 : if( !force && actionInChain() ) {
78 : return;
79 : }
80 339079 : ActionWithValue::clearDerivatives();
81 339079 : if( action_to_do_after ) {
82 184785 : action_to_do_after->clearDerivatives( true );
83 : }
84 : }
85 :
86 457848 : void ActionWithVector::clearInputForces( const bool& force ) {
87 457848 : if( !force && actionInChain() ) {
88 : return;
89 : }
90 339079 : ActionWithValue::clearInputForces();
91 339079 : if( action_to_do_after ) {
92 184785 : action_to_do_after->clearInputForces( true );
93 : }
94 : }
95 :
96 46 : const ActionWithVector* ActionWithVector::getFirstActionInChain() const {
97 46 : if( !actionInChain() ) {
98 : return this;
99 : }
100 46 : return action_to_do_before->getFirstActionInChain();
101 : }
102 :
103 209536157 : ActionWithVector* ActionWithVector::getFirstActionInChain() {
104 209536157 : if( !actionInChain() ) {
105 : return this;
106 : }
107 198647495 : return action_to_do_before->getFirstActionInChain();
108 : }
109 :
110 456611 : void ActionWithVector::retrieveAtoms( const bool& force ) {
111 456611 : if( !force && actionInChain() || atomsWereRetrieved ) {
112 : return;
113 : }
114 337503 : ActionAtomistic::retrieveAtoms();
115 337503 : atomsWereRetrieved = !actionInChain();
116 337503 : if( action_to_do_after ) {
117 184771 : action_to_do_after->retrieveAtoms( true );
118 : }
119 : }
120 :
121 46 : bool ActionWithVector::hasStoredArguments() const {
122 46 : std::string headstr=getFirstActionInChain()->getLabel();
123 100 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
124 68 : if( !getPntrToArgument(i)->ignoreStoredValue(headstr) ) {
125 : return true;
126 : }
127 : }
128 : return false;
129 : }
130 :
131 344 : bool ActionWithVector::argumentDependsOn( const std::string& headstr, ActionWithVector* faction, Value* thearg ) {
132 826 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
133 488 : if( this!=faction && thearg==getPntrToArgument(i) ) {
134 : return true;
135 : }
136 485 : ActionWithVector* av = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
137 485 : if( av && (av->getFirstActionInChain())->getLabel()==headstr ) {
138 180 : if( av->argumentDependsOn( headstr, faction, thearg ) ) {
139 : return true;
140 : };
141 : }
142 : }
143 : return false;
144 : }
145 :
146 3416 : unsigned ActionWithVector::buildArgumentStore( const unsigned& argstart ) {
147 : // Don't use chains for grids
148 10660 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
149 7626 : if( getPntrToArgument(i)->isConstant() ) {
150 901 : continue;
151 : }
152 6725 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction());
153 6725 : if( !av || getPntrToArgument(i)->getRank()>0 && getPntrToArgument(i)->hasDerivatives() ) {
154 382 : done_in_chain=false;
155 382 : break;
156 : }
157 : }
158 3416 : if( done_in_chain ) {
159 : std::vector<std::string> alabels;
160 : std::vector<ActionWithVector*> f_actions;
161 8656 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
162 : bool found=false;
163 6218 : std::string mylab = (getPntrToArgument(i)->getPntrToAction())->getLabel();
164 12821 : for(unsigned j=0; j<alabels.size(); ++j) {
165 8752 : if( alabels[j]==mylab ) {
166 : found=true;
167 : break;
168 : }
169 : }
170 6218 : if( !found ) {
171 4069 : alabels.push_back( mylab );
172 : }
173 :
174 : // If this is calculated in setup we never need to add to chain
175 6218 : if( getPntrToArgument(i)->isConstant() ) {
176 : continue;
177 : }
178 : // Find the chain we need to add this to from the arguments
179 5801 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction());
180 5801 : plumed_assert( av );
181 : found=false;
182 5801 : ActionWithVector* myact = av->getFirstActionInChain();
183 5801 : if( getPntrToArgument(i)->getRank()==1 && getPntrToArgument(i)->storedata ) {
184 5540 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
185 5252 : ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getPntrToArgument(j)->getPntrToAction() );
186 5252 : if( !aarg || i==j ) {
187 313 : continue;
188 : }
189 14182 : for(unsigned k=0; k<aarg->getNumberOfArguments(); ++k) {
190 9263 : if( aarg->getPntrToArgument(k)==getPntrToArgument(i) ) {
191 20 : done_in_chain=false;
192 20 : return reallyBuildArgumentStore( argstart );
193 : }
194 : }
195 : }
196 : }
197 6147 : for(unsigned j=0; j<f_actions.size(); ++j) {
198 3323 : if( f_actions[j]==myact ) {
199 : found=true;
200 : break;
201 : }
202 : }
203 5781 : if( !found ) {
204 2824 : if( f_actions.size()>0 ) {
205 336 : if( f_actions[0]->checkForDependency(myact) ) {
206 13 : getPntrToArgument(i)->buildDataStore();
207 : }
208 336 : if( myact->checkForDependency(f_actions[0]) ) {
209 0 : error("cannot deal with arguments in this order. Try swapping argument order");
210 : }
211 : }
212 2824 : if( !getPntrToArgument(i)->storedata && getPntrToArgument(i)->getRank()>0 ) {
213 2480 : f_actions.push_back( myact );
214 : }
215 : }
216 : }
217 : // Now make sure that everything we need is in the chain
218 2438 : if( f_actions.size()>0 ) {
219 : // Check everything for later f_actions is done before f_actions[0]
220 2463 : for(unsigned i=1; i<f_actions.size(); ++i) {
221 122 : ActionWithArguments* aarg = dynamic_cast<ActionWithArguments*>( f_actions[i] );
222 122 : if( !aarg || aarg->getNumberOfArguments()==0 ) {
223 66 : continue;
224 : }
225 179 : for(unsigned j=0; j<aarg->getNumberOfArguments(); ++j) {
226 134 : if( (aarg->getPntrToArgument(j))->isConstant() ) {
227 11 : continue ;
228 : }
229 : bool found=false;
230 123 : std::string dep_argname = (aarg->getPntrToArgument(j))->getPntrToAction()->getLabel();
231 17497 : for(const auto & pp : plumed.getActionSet()) {
232 : Action* p(pp.get());
233 : // Check if this is the dependency
234 17497 : if( p->getLabel()==dep_argname ) {
235 : found=true;
236 : break;
237 : }
238 : // Check if this is the first of the arguments that will appear in this chain
239 17385 : else if( p->getLabel()==f_actions[0]->getLabel() ) {
240 : break;
241 : }
242 : }
243 123 : if( !found ) {
244 11 : done_in_chain=false;
245 : break;
246 : }
247 : }
248 : // Stop trying to add things in the chain if we cannot
249 56 : if( !done_in_chain ) {
250 11 : return reallyBuildArgumentStore( argstart );
251 : }
252 : }
253 2341 : std::vector<std::string> empty(1);
254 2341 : empty[0] = f_actions[0]->getLabel();
255 2452 : for(unsigned i=1; i<f_actions.size(); ++i) {
256 111 : f_actions[0]->addActionToChain( empty, f_actions[i] );
257 : }
258 2341 : }
259 : // Now add this argument to the chain
260 : bool added=false;
261 2476 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
262 : // Add this function to jobs to do in recursive loop in previous action
263 2476 : if( getPntrToArgument(i)->getRank()>0 && !getPntrToArgument(i)->isConstant() ) {
264 2476 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
265 2476 : if( av && av->addActionToChain( alabels, this ) ) {
266 : added=true;
267 : break;
268 : }
269 : }
270 : }
271 2427 : plumed_massert(added, "could not add action " + getLabel() + " to chain of any of its arguments");
272 : // And get the number of derivatives
273 2427 : ActionWithVector* head=getFirstActionInChain();
274 2427 : unsigned nder=0;
275 2427 : arg_deriv_starts.resize( getNumberOfArguments() );
276 8598 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
277 : // Ensures that we ignore the grid (first argument) if we are evaluating a function on a grid
278 6171 : if( i==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) {
279 1 : continue ;
280 : }
281 6170 : ActionWithVector* iaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction());
282 6170 : if( actionInChain() && !getPntrToArgument(i)->ignoreStoredValue(head->getLabel()) ) {
283 675 : arg_deriv_starts[i] = 0;
284 675 : head->getNumberOfStreamedDerivatives( arg_deriv_starts[i], getPntrToArgument(i) );
285 5495 : } else if( iaction && iaction->isInSubChain(nder) ) {
286 2255 : arg_deriv_starts[i] = nder;
287 : // Add the total number of derivatives that we have by this point in the chain to nder
288 : if( iaction ) {
289 2255 : nder=0;
290 2255 : head->getNumberOfStreamedDerivatives( nder, getPntrToArgument(i) );
291 : }
292 : } else {
293 : // Check if we have already found this action
294 : int k=-1;
295 3240 : if( iaction ) {
296 3240 : ActionWithVector* ider_action=iaction->getActionWithDerivatives( iaction );
297 3461 : for(unsigned j=0; j<i; ++j) {
298 1631 : if( j==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) {
299 1 : continue ;
300 : }
301 1630 : ActionWithVector* jaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(j)->getPntrToAction());
302 1630 : if( jaction->getActionWithDerivatives(jaction)==ider_action || jaction->checkForDependency(ider_action) ) {
303 1410 : k=j;
304 1410 : break;
305 : }
306 : }
307 3240 : if( k>=0 ) {
308 1410 : arg_deriv_starts[i] = arg_deriv_starts[k];
309 1410 : continue;
310 : }
311 : }
312 :
313 1830 : if( i>0 ) {
314 : // This is a fudge so that inputs like this work:
315 : // c: CONTACT_MATRIX ATOMS=1-100
316 : // d: MATRIX_PRODUCT ARG=mat1,mat2
317 : // e: CUSTOM ARG=c,d
318 : // f: MATRIX_PRODUCT ARG=mat3,mat4
319 : // g: CUSTOM ARG=c,f
320 : // See symfunc rt-nbonds-q6 for an example
321 : // In this example when we set arg_deriv_starts[1] for f in g nder=number of derivatives of c
322 : // mder is equal to the number of derivatives by the time you get to f minus the number of derivatives for c
323 136 : unsigned mder=0;
324 136 : ActionWithVector* jaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(i-1)->getPntrToAction());
325 136 : if( jaction->action_to_do_after && !(jaction->action_to_do_after)->getNumberOfStoredValues( getPntrToArgument(i-1), mder, i, getArguments() ) ) {
326 0 : mder=0;
327 : }
328 136 : if( mder>0 ) {
329 1 : nder = nder + mder;
330 : }
331 : }
332 :
333 1830 : arg_deriv_starts[i] = nder;
334 : // Add the total number of derivatives that we have by this point in the chain to nder
335 1830 : if( iaction ) {
336 1830 : nder=0;
337 1830 : if( (getPntrToArgument(i)->getPntrToAction())->getName().find("DIFFERENCE")!=std::string::npos ) {
338 0 : ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getPntrToArgument(i)->getPntrToAction() );
339 0 : plumed_assert( aarg && aarg->getNumberOfArguments()==2 );
340 0 : head->getNumberOfStreamedDerivatives( nder, aarg->getPntrToArgument(0) );
341 0 : nder += (aarg->getPntrToArgument(1))->getNumberOfValues();
342 : } else {
343 1830 : head->getNumberOfStreamedDerivatives( nder, getPntrToArgument(i) );
344 : }
345 : }
346 : }
347 : }
348 2427 : nder=0;
349 2427 : head->getNumberOfStreamedDerivatives( nder, NULL );
350 2427 : return nder;
351 2458 : }
352 958 : return reallyBuildArgumentStore( argstart );
353 : }
354 :
355 989 : unsigned ActionWithVector::reallyBuildArgumentStore( const unsigned& argstart ) {
356 2703 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
357 1714 : if( getPntrToArgument(i)->getRank()>0 ) {
358 1468 : getPntrToArgument(i)->buildDataStore();
359 : }
360 : }
361 : unsigned nder=0;
362 989 : arg_deriv_starts.resize( getNumberOfArguments() );
363 2703 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
364 1714 : arg_deriv_starts[i] = nder;
365 1714 : nder += getPntrToArgument(i)->getNumberOfValues();
366 : }
367 989 : return nder;
368 : }
369 :
370 17724 : ActionWithVector* ActionWithVector::getActionWithDerivatives( ActionWithVector* depaction ) {
371 17724 : if( depaction==this || depaction->checkForDependency(this) ) {
372 11208 : if( getNumberOfAtoms()>0 ) {
373 4870 : return this;
374 : }
375 6798 : std::string c=getFirstActionInChain()->getLabel();
376 42678 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
377 36340 : if( !getPntrToArgument(i)->ignoreStoredValue(c) && !getPntrToArgument(i)->isConstant() ) {
378 : return this;
379 : }
380 : }
381 : }
382 12854 : plumed_assert( action_to_do_before );
383 12854 : return action_to_do_before->getActionWithDerivatives(depaction);
384 : }
385 :
386 11777 : bool ActionWithVector::addActionToChain( const std::vector<std::string>& alabels, ActionWithVector* act ) {
387 11777 : if( action_to_do_after ) {
388 9190 : bool state=action_to_do_after->addActionToChain( alabels, act );
389 9190 : return state;
390 : }
391 : // Check action is not already in chain
392 : std::vector<std::string> mylabels;
393 2587 : getFirstActionInChain()->getAllActionLabelsInChain( mylabels );
394 21253 : for(unsigned i=0; i<mylabels.size(); ++i) {
395 18666 : if( act->getLabel()==mylabels[i] ) {
396 : return true;
397 : }
398 : }
399 :
400 : // Check that everything that is required has been calculated
401 6786 : for(unsigned i=0; i<alabels.size(); ++i) {
402 : bool found=false;
403 20794 : for(unsigned j=0; j<mylabels.size(); ++j) {
404 20064 : if( alabels[i]==mylabels[j] ) {
405 : found=true;
406 : break;
407 : }
408 : }
409 4248 : if( !found ) {
410 730 : ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( alabels[i] );
411 730 : plumed_massert( av, "could not cast " + alabels[i] );
412 : bool storingall=true;
413 1522 : for(int j=0; j<av->getNumberOfComponents(); ++j) {
414 792 : if( !(av->getPntrToComponent(j))->storedata ) {
415 : storingall=false;
416 : }
417 : }
418 730 : if( !storingall ) {
419 : return false;
420 : }
421 : }
422 : }
423 : // This checks that there is nothing that will cause problems in the chain
424 2538 : mylabels.resize(0);
425 2538 : getFirstActionInChain()->getAllActionLabelsInChain( mylabels );
426 20751 : for(unsigned i=0; i<mylabels.size(); ++i) {
427 18213 : ActionWithVector* av1=plumed.getActionSet().selectWithLabel<ActionWithVector*>( mylabels[i] );
428 176690 : for(unsigned j=0; j<i; ++j) {
429 158477 : ActionWithVector* av2=plumed.getActionSet().selectWithLabel<ActionWithVector*>( mylabels[j] );
430 158477 : if( !av1->canBeAfterInChain( av2 ) ) {
431 0 : error("must calculate " + mylabels[j] + " before " + mylabels[i] );
432 : }
433 : }
434 : }
435 2538 : action_to_do_after=act;
436 2538 : act->action_to_do_before=this;
437 2538 : updateTaskListReductionStatus();
438 2538 : ActionWithVector* head = getFirstActionInChain();
439 2538 : head->broadcastThatTasksAreReduced( head );
440 2538 : head->finishChainBuild( act );
441 : return true;
442 2587 : }
443 :
444 2689 : void ActionWithVector::updateTaskListReductionStatus() {
445 2689 : ActionWithVector* head = getFirstActionInChain();
446 : std::vector<ActionWithVector*> task_reducing_actions;
447 2689 : head->canReduceTasks( task_reducing_actions );
448 2689 : if( task_reducing_actions.size()>0 ) {
449 417 : head->reduce_tasks=true;
450 : }
451 2689 : }
452 :
453 3843560 : void ActionWithVector::broadcastThatTasksAreReduced( ActionWithVector* aselect ) {
454 3843560 : std::string c=getFirstActionInChain()->getLabel();
455 18085469 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
456 14241909 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
457 71210 : ActionWithVector* av = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
458 71210 : if( av ) {
459 : bool found=false;
460 12111 : ActionWithVector* av_head = av->getFirstActionInChain();
461 17731 : for(unsigned i=0; i<av_head->task_control_list.size(); ++i) {
462 17268 : if( aselect==av_head->task_control_list[i] ) {
463 : found=true;
464 : break;
465 : }
466 : }
467 12111 : if( !found ) {
468 463 : av_head->task_control_list.insert( av_head->task_control_list.begin(), aselect );
469 : }
470 :
471 12111 : av_head->reduce_tasks=true;
472 12111 : av_head->updateTaskReductionFlag( av_head->reduce_tasks );
473 : }
474 : }
475 : }
476 3843560 : if( action_to_do_after ) {
477 3611558 : action_to_do_after->broadcastThatTasksAreReduced( aselect );
478 : }
479 3843560 : }
480 :
481 229464 : void ActionWithVector::updateTaskReductionFlag( bool& head_reduce_tasks ) {
482 229464 : if( actionInChain() ) {
483 217353 : plumed_assert( task_control_list.size()==0 );
484 : } else {
485 30767 : for(unsigned i=0; i<task_control_list.size(); ++i) {
486 18656 : if( !(task_control_list[i]->getFirstActionInChain())->reduce_tasks ) {
487 16928 : head_reduce_tasks=false;
488 : }
489 : }
490 : }
491 229464 : broadcastThatTasksAreReduced( getFirstActionInChain() );
492 229464 : if( action_to_do_after ) {
493 217353 : action_to_do_after->updateTaskReductionFlag( head_reduce_tasks );
494 : }
495 229464 : }
496 :
497 20958 : void ActionWithVector::canReduceTasks( std::vector<ActionWithVector*>& task_reducing_actions ) {
498 20958 : areAllTasksRequired( task_reducing_actions );
499 20958 : if( action_to_do_after ) {
500 18269 : action_to_do_after->canReduceTasks( task_reducing_actions );
501 : }
502 20958 : }
503 :
504 3586 : void ActionWithVector::finishChainBuild( ActionWithVector* act ) {
505 3586 : if( action_to_do_after ) {
506 2983 : action_to_do_after->finishChainBuild( act );
507 : }
508 3586 : }
509 :
510 36959 : void ActionWithVector::getAllActionLabelsInChain( std::vector<std::string>& mylabels ) const {
511 : bool found = false ;
512 359399 : for(unsigned i=0; i<mylabels.size(); ++i) {
513 322440 : if( getLabel()==mylabels[i] ) {
514 : found=true;
515 : }
516 : }
517 36959 : if( !found ) {
518 36959 : mylabels.push_back( getLabel() );
519 : }
520 36959 : if( action_to_do_after ) {
521 31823 : action_to_do_after->getAllActionLabelsInChain( mylabels );
522 : }
523 36959 : }
524 :
525 7349456 : void ActionWithVector::taskIsActive( const unsigned& current, int& flag ) const {
526 7349456 : if( isActive() ) {
527 7081984 : flag = checkTaskStatus( current, flag );
528 : }
529 7349456 : if( flag<=0 && action_to_do_after ) {
530 4891432 : action_to_do_after->taskIsActive( current, flag );
531 : }
532 7349456 : }
533 :
534 1345 : void ActionWithVector::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
535 1377 : for(unsigned i=0; i<task_control_list.size(); ++i ) {
536 32 : task_control_list[i]->getAdditionalTasksRequired( action, atasks );
537 : }
538 1345 : }
539 :
540 274656 : void ActionWithVector::prepare() {
541 274656 : active_tasks.resize(0);
542 274656 : atomsWereRetrieved=false;
543 274656 : }
544 :
545 144885 : std::vector<unsigned>& ActionWithVector::getListOfActiveTasks( ActionWithVector* action ) {
546 144885 : if( active_tasks.size()>0 ) {
547 11990 : return active_tasks;
548 : }
549 132895 : unsigned ntasks=0;
550 132895 : getNumberOfTasks( ntasks );
551 :
552 132895 : unsigned stride=comm.Get_size();
553 132895 : unsigned rank=comm.Get_rank();
554 132895 : if(serial) {
555 : stride=1;
556 : rank=0;
557 : }
558 :
559 : // Get number of threads for OpenMP
560 132895 : unsigned nt=OpenMP::getNumThreads();
561 132895 : if( nt*stride*10>ntasks ) {
562 55787 : nt=ntasks/stride/10;
563 : }
564 : if( nt==0 ) {
565 : nt=1;
566 : }
567 :
568 132895 : if( !never_reduce_tasks && reduce_tasks ) {
569 1529 : if( task_control_list.size()>0 ) {
570 : // Get the list of tasks that are active in the action that uses the output of this action
571 238 : for(unsigned i=0; i<task_control_list.size(); ++i) {
572 130 : task_control_list[i]->retrieveAtoms();
573 130 : active_tasks = task_control_list[i]->getListOfActiveTasks( action );
574 : }
575 : // Now work out else we need from here to calculate the later action
576 108 : getAdditionalTasksRequired( action, active_tasks );
577 : } else {
578 1421 : std::vector<int> taskFlags( ntasks, -1 );
579 :
580 1421 : #pragma omp parallel num_threads(nt)
581 : {
582 : #pragma omp for nowait
583 : for(unsigned i=rank; i<ntasks; i+=stride ) {
584 : taskIsActive( i, taskFlags[i] );
585 : }
586 : }
587 2658057 : for(unsigned i=0; i<ntasks; ++i) {
588 2656636 : taskFlags[i] = std::abs( taskFlags[i] );
589 : }
590 1421 : if( !serial ) {
591 1421 : comm.Sum( taskFlags );
592 : }
593 :
594 : unsigned nt=0;
595 2658057 : for(unsigned i=0; i<ntasks; ++i) {
596 2656636 : if( taskFlags[i]>=stride ) {
597 264019 : nt++;
598 : }
599 : }
600 1421 : active_tasks.resize(nt);
601 : nt=0;
602 2658057 : for(unsigned i=0; i<ntasks; ++i) {
603 2656636 : if( taskFlags[i]>=stride ) {
604 264019 : active_tasks[nt]=i;
605 264019 : nt++;
606 : }
607 : }
608 1421 : getAdditionalTasksRequired( this, active_tasks );
609 : }
610 : } else {
611 131366 : active_tasks.resize( ntasks );
612 5116683 : for(unsigned i=0; i<ntasks; ++i) {
613 4985317 : active_tasks[i]=i;
614 : }
615 : }
616 132895 : return active_tasks;
617 : }
618 :
619 133040 : void ActionWithVector::runAllTasks() {
620 : // Skip this if this is done elsewhere
621 133040 : if( action_to_do_before ) {
622 35 : return;
623 : }
624 :
625 133005 : unsigned stride=comm.Get_size();
626 133005 : unsigned rank=comm.Get_rank();
627 133005 : if(serial) {
628 : stride=1;
629 : rank=0;
630 : }
631 :
632 : // Get the list of active tasks
633 133005 : std::vector<unsigned> & partialTaskList( getListOfActiveTasks( this ) );
634 133005 : unsigned nactive_tasks=partialTaskList.size();
635 :
636 : // Get number of threads for OpenMP
637 133005 : unsigned nt=OpenMP::getNumThreads();
638 133005 : if( nt*stride*10>nactive_tasks ) {
639 56005 : nt=nactive_tasks/stride/10;
640 : }
641 133005 : if( nt==0 ) {
642 18865 : nt=1;
643 : }
644 :
645 : // Now do all preparations required to run all the tasks
646 : // prepareForTaskLoop();
647 :
648 : // Get the total number of streamed quantities that we need
649 133005 : unsigned nquants=0, nmatrices=0, maxcol=0, nbooks=0;
650 133005 : getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol, nbooks );
651 : // Get size for buffer
652 133005 : unsigned bufsize=0;
653 133005 : getSizeOfBuffer( nactive_tasks, bufsize );
654 133005 : if( buffer.size()!=bufsize ) {
655 2246 : buffer.resize( bufsize );
656 : }
657 : // Clear buffer
658 133005 : buffer.assign( buffer.size(), 0.0 );
659 :
660 : // Recover the number of derivatives we require
661 133005 : unsigned nderivatives = 0;
662 133005 : bool gridsInStream=checkForGrids(nderivatives);
663 133005 : if( !doNotCalculateDerivatives() && !gridsInStream ) {
664 105266 : getNumberOfStreamedDerivatives( nderivatives, NULL );
665 : }
666 :
667 133005 : #pragma omp parallel num_threads(nt)
668 : {
669 : std::vector<double> omp_buffer;
670 : if( nt>1 ) {
671 : omp_buffer.resize( bufsize, 0.0 );
672 : }
673 : MultiValue myvals( nquants, nderivatives, nmatrices, maxcol, nbooks );
674 : myvals.clearAll();
675 :
676 : #pragma omp for nowait
677 : for(unsigned i=rank; i<nactive_tasks; i+=stride) {
678 : // Calculate the stuff in the loop for this action
679 : runTask( partialTaskList[i], myvals );
680 :
681 : // Now transfer the data to the actions that accumulate values from the calculated quantities
682 : if( nt>1 ) {
683 : gatherAccumulators( partialTaskList[i], myvals, omp_buffer );
684 : } else {
685 : gatherAccumulators( partialTaskList[i], myvals, buffer );
686 : }
687 :
688 : // Clear the value
689 : myvals.clearAll();
690 : }
691 : #pragma omp critical
692 : gatherThreads( nt, bufsize, omp_buffer, buffer, myvals );
693 : }
694 :
695 : // MPI Gather everything
696 133005 : if( !serial && buffer.size()>0 ) {
697 132998 : gatherProcesses( buffer );
698 : }
699 133005 : finishComputations( buffer );
700 : }
701 :
702 210005 : void ActionWithVector::gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) {
703 210005 : if( nt>1 )
704 100237334 : for(unsigned i=0; i<bufsize; ++i) {
705 100083334 : buffer[i]+=omp_buffer[i];
706 : }
707 210005 : }
708 :
709 132998 : void ActionWithVector::gatherProcesses( std::vector<double>& buffer ) {
710 132998 : comm.Sum( buffer );
711 132998 : }
712 :
713 317809 : bool ActionWithVector::checkForGrids( unsigned& nder ) const {
714 650851 : for(int i=0; i<getNumberOfComponents(); ++i) {
715 341775 : if( getConstPntrToComponent(i)->getRank()>0 && getConstPntrToComponent(i)->hasDerivatives() ) {
716 8733 : nder=getConstPntrToComponent(i)->getNumberOfGridDerivatives();
717 8733 : return true;
718 : }
719 : }
720 309076 : if( action_to_do_after ) {
721 184804 : return action_to_do_after->checkForGrids(nder);
722 : }
723 : return false;
724 : }
725 :
726 315090 : void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) {
727 315090 : if( ntasks==0 ) {
728 130242 : if( getNumberOfArguments()==1 && getNumberOfComponents()==1 && getPntrToComponent(0)->getRank()==0 ) {
729 2285 : if( !getPntrToArgument(0)->hasDerivatives() && getPntrToArgument(0)->getRank()==2 ) {
730 17 : ntasks = getPntrToArgument(0)->getShape()[0];
731 : } else {
732 2268 : ntasks = getPntrToArgument(0)->getNumberOfValues();
733 : }
734 : } else {
735 127957 : plumed_assert( getNumberOfComponents()>0 && getPntrToComponent(0)->getRank()>0 );
736 127957 : if( getPntrToComponent(0)->hasDerivatives() ) {
737 6061 : ntasks = getPntrToComponent(0)->getNumberOfValues();
738 : } else {
739 121896 : ntasks = getPntrToComponent(0)->getShape()[0];
740 : }
741 : }
742 : }
743 654142 : for(int i=0; i<getNumberOfComponents(); ++i) {
744 339052 : if( getPntrToComponent(i)->getRank()==0 ) {
745 65973 : if( getNumberOfArguments()!=1 ) {
746 0 : error("mismatched numbers of tasks in streamed quantities");
747 : }
748 65973 : if( getPntrToArgument(0)->hasDerivatives() && ntasks!=getPntrToArgument(0)->getNumberOfValues() ) {
749 0 : error("mismatched numbers of tasks in streamed quantities");
750 65973 : } else if ( !getPntrToArgument(0)->hasDerivatives() && ntasks!=getPntrToArgument(0)->getShape()[0] ) {
751 0 : error("mismatched numbers of tasks in streamed quantities");
752 : }
753 273079 : } else if( getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getNumberOfValues() ) {
754 0 : error("mismatched numbers of tasks in streamed quantities");
755 273079 : } else if( !getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getShape()[0] ) {
756 0 : error("mismatched numbers of tasks in streamed quantities");
757 : }
758 : }
759 315090 : if( action_to_do_after ) {
760 184848 : action_to_do_after->getNumberOfTasks( ntasks );
761 : }
762 315090 : }
763 :
764 421404 : void ActionWithVector::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
765 1075453 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
766 654049 : if( !getPntrToArgument(i)->ignoreStoredValue(headstr) ) {
767 412404 : getPntrToArgument(i)->streampos=nquants;
768 412404 : nquants++;
769 : }
770 : }
771 881865 : for(int i=0; i<getNumberOfComponents(); ++i) {
772 460461 : getPntrToComponent(i)->streampos=nquants;
773 460461 : nquants++;
774 : }
775 421404 : }
776 :
777 421404 : void ActionWithVector::getNumberOfStreamedQuantities( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
778 421404 : setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
779 421404 : if( action_to_do_after ) {
780 213562 : action_to_do_after->getNumberOfStreamedQuantities( headstr, nquants, nmat, maxcol, nbookeeping );
781 : }
782 421404 : }
783 :
784 317809 : void ActionWithVector::getSizeOfBuffer( const unsigned& nactive_tasks, unsigned& bufsize ) {
785 659584 : for(int i=0; i<getNumberOfComponents(); ++i) {
786 341775 : getPntrToComponent(i)->bufstart=bufsize;
787 341775 : bufsize += getPntrToComponent(i)->data.size();
788 : }
789 317809 : if( action_to_do_after ) {
790 184804 : action_to_do_after->getSizeOfBuffer( nactive_tasks, bufsize );
791 : }
792 317809 : }
793 :
794 424606 : void ActionWithVector::getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ) {
795 424606 : std::string c=getFirstActionInChain()->getLabel();
796 1145263 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
797 721332 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
798 385865 : if( getPntrToArgument(i)==stopat ) {
799 : return;
800 : }
801 385190 : nderivatives += getPntrToArgument(i)->getNumberOfValues();
802 : }
803 : }
804 423931 : if( getNumberOfAtoms()>0 ) {
805 54800 : nderivatives += 3*getNumberOfAtoms() + 9;
806 : }
807 : // Don't do the whole chain if we have been told to stop early
808 423931 : if( stopat && stopat->getPntrToAction()==this ) {
809 : return;
810 : }
811 :
812 419742 : if( action_to_do_after ) {
813 237212 : action_to_do_after->getNumberOfStreamedDerivatives( nderivatives, stopat );
814 : }
815 : }
816 :
817 136 : bool ActionWithVector::getNumberOfStoredValues( Value* startat, unsigned& nvals, const unsigned& astart, const std::vector<Value*>& stopat ) {
818 181 : for(unsigned j=astart; j<stopat.size(); ++j) {
819 151 : if( stopat[j] && (stopat[j]->getPntrToAction()==this || (stopat[j]->getPntrToAction())->checkForDependency(this)) ) {
820 106 : return true;
821 : }
822 : }
823 :
824 30 : std::string c=getFirstActionInChain()->getLabel();
825 77 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
826 47 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
827 4 : for(unsigned j=astart; j<stopat.size(); ++j) {
828 2 : if( getPntrToArgument(i)==stopat[j] ) {
829 : return true;
830 : }
831 : }
832 2 : nvals += getPntrToArgument(i)->getNumberOfValues();
833 : }
834 : }
835 30 : if( startat->getPntrToAction()!=this && getNumberOfAtoms()>0 ) {
836 : return false;
837 : }
838 :
839 30 : if( action_to_do_after ) {
840 30 : return action_to_do_after->getNumberOfStoredValues( startat, nvals, astart, stopat );
841 : }
842 : return false;
843 : }
844 :
845 11395074 : void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) const {
846 11395074 : if( isActive() ) {
847 10279025 : myvals.setTaskIndex(current);
848 10279025 : myvals.vector_call=true;
849 10279025 : performTask( current, myvals );
850 : }
851 11395074 : if( action_to_do_after ) {
852 4886746 : action_to_do_after->runTask( current, myvals );
853 : }
854 11395074 : }
855 :
856 9293892 : void ActionWithVector::gatherAccumulators( const unsigned& taskCode, const MultiValue& myvals, std::vector<double>& buffer ) const {
857 9293892 : if( isActive() ) {
858 18680617 : for(int i=0; i<getNumberOfComponents(); ++i) {
859 10347771 : unsigned bufstart = getConstPntrToComponent(i)->bufstart;
860 : // This looks after storing of scalars that are summed from vectors/matrices
861 10347771 : if( getConstPntrToComponent(i)->getRank()==0 ) {
862 : plumed_dbg_massert( bufstart<buffer.size(), "problem in " + getLabel() );
863 1161228 : unsigned sind = getConstPntrToComponent(i)->streampos;
864 1161228 : buffer[bufstart] += myvals.get(sind);
865 1161228 : if( getConstPntrToComponent(i)->hasDerivatives() ) {
866 28422118 : for(unsigned k=0; k<myvals.getNumberActive(sind); ++k) {
867 27260890 : unsigned kindex = myvals.getActiveIndex(sind,k);
868 : plumed_dbg_massert( bufstart+1+kindex<buffer.size(), "problem in " + getLabel() );
869 27260890 : buffer[bufstart + 1 + kindex] += myvals.getDerivative(sind,kindex);
870 : }
871 : }
872 : // This looks after storing of vectors
873 9186543 : } else if( getConstPntrToComponent(i)->storedata ) {
874 5062271 : gatherStoredValue( i, taskCode, myvals, bufstart, buffer );
875 : }
876 : }
877 : }
878 9293892 : if( action_to_do_after ) {
879 4206354 : action_to_do_after->gatherAccumulators( taskCode, myvals, buffer );
880 : }
881 9293892 : }
882 :
883 3571815 : void ActionWithVector::gatherStoredValue( const unsigned& valindex, const unsigned& taskCode, const MultiValue& myvals,
884 : const unsigned& bufstart, std::vector<double>& buffer ) const {
885 : plumed_dbg_assert( getConstPntrToComponent(valindex)->getRank()==1 && !getConstPntrToComponent(valindex)->hasDeriv );
886 3571815 : unsigned vindex = getConstPntrToComponent(valindex)->bufstart + taskCode;
887 : plumed_dbg_massert( vindex<buffer.size(), "failing in " + getLabel() );
888 3571815 : buffer[vindex] += myvals.get(getConstPntrToComponent(valindex)->streampos);
889 3571815 : }
890 :
891 317809 : void ActionWithVector::finishComputations( const std::vector<double>& buf ) {
892 317809 : if( isActive() ) {
893 527552 : for(int i=0; i<getNumberOfComponents(); ++i) {
894 : // This gathers vectors and grids at the end of the calculation
895 275759 : unsigned bufstart = getPntrToComponent(i)->bufstart;
896 275759 : getPntrToComponent(i)->data.assign( getPntrToComponent(i)->data.size(), 0 );
897 275759 : if( (getPntrToComponent(i)->getRank()>0 && getPntrToComponent(i)->hasDerivatives()) || getPntrToComponent(i)->storedata ) {
898 130918 : unsigned sz_v = getPntrToComponent(i)->data.size();
899 46583551 : for(unsigned j=0; j<sz_v; ++j) {
900 : plumed_dbg_assert( bufstart+j<buf.size() );
901 46452633 : getPntrToComponent(i)->add( j, buf[bufstart+j] );
902 : }
903 : // Make sure single values are set
904 144841 : } else if( getPntrToComponent(i)->getRank()==0 ) {
905 43495 : getPntrToComponent(i)->set( buf[bufstart] );
906 : }
907 : // This gathers derivatives of scalars
908 275759 : if( !doNotCalculateDerivatives() && getPntrToComponent(i)->hasDeriv && getPntrToComponent(i)->getRank()==0 ) {
909 16307287 : for(unsigned j=0; j<getPntrToComponent(i)->getNumberOfDerivatives(); ++j) {
910 16271018 : getPntrToComponent(i)->setDerivative( j, buf[bufstart+1+j] );
911 : }
912 : }
913 : }
914 : }
915 317809 : if( action_to_do_after ) {
916 184804 : action_to_do_after->finishComputations( buf );
917 : }
918 317809 : }
919 :
920 293010 : bool ActionWithVector::checkChainForNonScalarForces() const {
921 523450 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
922 305277 : if( getConstPntrToComponent(i)->getRank()>0 && getConstPntrToComponent(i)->forcesWereAdded() ) {
923 : return true;
924 : }
925 : }
926 218173 : if( action_to_do_after ) {
927 167710 : return action_to_do_after->checkChainForNonScalarForces();
928 : }
929 : return false;
930 : }
931 :
932 254548 : bool ActionWithVector::checkForForces() {
933 254548 : if( getPntrToComponent(0)->getRank()==0 ) {
934 51808 : return ActionWithValue::checkForForces();
935 202740 : } else if( actionInChain() ) {
936 : return false;
937 : }
938 :
939 : // Check if there are any forces
940 125300 : if( !checkChainForNonScalarForces() ) {
941 : return false;
942 : }
943 :
944 : // Setup MPI parallel loop
945 74837 : unsigned stride=comm.Get_size();
946 74837 : unsigned rank=comm.Get_rank();
947 74837 : if(serial) {
948 : stride=1;
949 : rank=0;
950 : }
951 :
952 : // Get the number of tasks
953 : std::vector<unsigned> force_tasks;
954 74837 : getForceTasks( force_tasks );
955 74837 : Tools::removeDuplicates(force_tasks);
956 74837 : unsigned nf_tasks=force_tasks.size();
957 :
958 : // Get number of threads for OpenMP
959 74837 : unsigned nt=OpenMP::getNumThreads();
960 74837 : if( nt*stride*10>nf_tasks ) {
961 42255 : nt=nf_tasks/stride/10;
962 : }
963 : if( nt==0 ) {
964 : nt=1;
965 : }
966 :
967 : // Now determine how big the multivalue needs to be
968 74837 : unsigned nquants=0, nmatrices=0, maxcol=0, nbooks=0;
969 74837 : getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol, nbooks );
970 : // Recover the number of derivatives we require (this should be equal to the number of forces)
971 74837 : unsigned nderiv=0;
972 74837 : getNumberOfStreamedDerivatives( nderiv, NULL );
973 74837 : if( forcesForApply.size()!=nderiv ) {
974 351 : forcesForApply.resize( nderiv );
975 : }
976 : // Clear force buffer
977 74837 : forcesForApply.assign( forcesForApply.size(), 0.0 );
978 :
979 74837 : #pragma omp parallel num_threads(nt)
980 : {
981 : std::vector<double> omp_forces;
982 : if( nt>1 ) {
983 : omp_forces.resize( forcesForApply.size(), 0.0 );
984 : }
985 : MultiValue myvals( nquants, nderiv, nmatrices, maxcol, nbooks );
986 : myvals.clearAll();
987 :
988 : #pragma omp for nowait
989 : for(unsigned i=rank; i<nf_tasks; i+=stride) {
990 : runTask( force_tasks[i], myvals );
991 :
992 : // Now get the forces
993 : if( nt>1 ) {
994 : gatherForces( force_tasks[i], myvals, omp_forces );
995 : } else {
996 : gatherForces( force_tasks[i], myvals, forcesForApply );
997 : }
998 :
999 : myvals.clearAll();
1000 : }
1001 : #pragma omp critical
1002 : if(nt>1)
1003 : for(unsigned i=0; i<forcesForApply.size(); ++i) {
1004 : forcesForApply[i]+=omp_forces[i];
1005 : }
1006 : }
1007 : // MPI Gather on forces
1008 74837 : if( !serial ) {
1009 74837 : comm.Sum( forcesForApply );
1010 : }
1011 : return true;
1012 : }
1013 :
1014 2039456 : bool ActionWithVector::checkComponentsForForce() const {
1015 2777915 : for(unsigned i=0; i<values.size(); ++i) {
1016 2546733 : if( getConstPntrToComponent(i)->forcesWereAdded() ) {
1017 : return true;
1018 : }
1019 : }
1020 : return false;
1021 : }
1022 :
1023 5561312 : bool ActionWithVector::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
1024 : plumed_dbg_assert( !(myval->getRank()==2 && !myval->hasDerivatives()) );
1025 5561312 : return fabs(myval->getForce(itask))>epsilon;
1026 : }
1027 :
1028 90719 : void ActionWithVector::updateForceTasksFromValue( const Value* myval, std::vector<unsigned>& force_tasks ) const {
1029 90719 : if( myval->getRank()>0 && myval->forcesWereAdded() ) {
1030 85478 : unsigned nt = myval->getNumberOfValues();
1031 85478 : if( !myval->hasDerivatives() ) {
1032 85009 : nt = myval->getShape()[0];
1033 : }
1034 6149611 : for(unsigned i=0; i<nt; ++i) {
1035 6064133 : if( checkForTaskForce(i, myval) ) {
1036 2074907 : force_tasks.push_back( i );
1037 : }
1038 : }
1039 : }
1040 90719 : }
1041 :
1042 103595 : void ActionWithVector::getForceTasks( std::vector<unsigned>& force_tasks ) const {
1043 103595 : if( isActive() && checkComponentsForForce() ) {
1044 171577 : for(unsigned k=0; k<values.size(); ++k) {
1045 91329 : updateForceTasksFromValue( getConstPntrToComponent(k), force_tasks );
1046 : }
1047 : }
1048 103595 : if( action_to_do_after ) {
1049 28758 : action_to_do_after->getForceTasks( force_tasks );
1050 : }
1051 103595 : }
1052 :
1053 1565796 : void ActionWithVector::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
1054 : plumed_dbg_assert( myval->storedata && !(myval->getRank()==2 && !myval->hasDerivatives()) );
1055 1565796 : double fforce = myval->getForce(itask);
1056 : unsigned sspos = myval->getPositionInStream();
1057 20912845 : for(unsigned j=0; j<myvals.getNumberActive(sspos); ++j) {
1058 19347049 : unsigned jder=myvals.getActiveIndex(sspos, j);
1059 : plumed_dbg_assert( jder<forces.size() );
1060 19347049 : forces[jder] += fforce*myvals.getDerivative( sspos, jder );
1061 : }
1062 1565796 : }
1063 :
1064 2101182 : void ActionWithVector::gatherForces( const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
1065 2101182 : if( isActive() && checkComponentsForForce() ) {
1066 3587802 : for(unsigned k=0; k<getNumberOfComponents(); ++k) {
1067 1859776 : const Value* myval=getConstPntrToComponent(k);
1068 1859776 : if( myval->getRank()>0 && myval->forcesWereAdded() ) {
1069 1694828 : gatherForcesOnStoredValue( myval, itask, myvals, forces );
1070 : }
1071 : }
1072 : }
1073 2101182 : if( action_to_do_after ) {
1074 680392 : action_to_do_after->gatherForces( itask, myvals, forces );
1075 : }
1076 2101182 : }
1077 :
1078 254548 : void ActionWithVector::apply() {
1079 254548 : if( !checkForForces() ) {
1080 136681 : return;
1081 : }
1082 : // Find the top of the chain and add forces
1083 117867 : unsigned ind=0;
1084 117867 : getFirstActionInChain()->addForcesToInput( getForcesToApply(), ind );
1085 : }
1086 :
1087 186439 : void ActionWithVector::addForcesToInput( const std::vector<double>& forcesToApply, unsigned& ind ) {
1088 186439 : if( ind>=forcesToApply.size() ) {
1089 : return;
1090 : }
1091 136857 : addForcesOnArguments( 0, forcesToApply, ind, getFirstActionInChain()->getLabel() );
1092 136857 : setForcesOnAtoms( forcesToApply, ind );
1093 136857 : if( action_to_do_after ) {
1094 68572 : action_to_do_after->addForcesToInput( forcesToApply, ind );
1095 : }
1096 : }
1097 :
1098 : }
|
