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 5380 : void ActionWithVector::registerKeywords( Keywords& keys ) {
31 5380 : Action::registerKeywords( keys );
32 5380 : ActionAtomistic::registerKeywords( keys );
33 5380 : ActionWithValue::registerKeywords( keys ); keys.remove("NUMERICAL_DERIVATIVES");
34 5380 : ActionWithArguments::registerKeywords( keys );
35 10760 : keys.addFlag("SERIAL",false,"do the calculation in serial. Do not parallelize");
36 5380 : }
37 :
38 4916 : ActionWithVector::ActionWithVector(const ActionOptions&ao):
39 : Action(ao),
40 : ActionAtomistic(ao),
41 : ActionWithValue(ao),
42 : ActionWithArguments(ao),
43 4916 : serial(false),
44 4916 : action_to_do_before(NULL),
45 4916 : action_to_do_after(NULL),
46 4916 : never_reduce_tasks(false),
47 4916 : reduce_tasks(false),
48 4916 : atomsWereRetrieved(false),
49 4916 : done_in_chain(false)
50 : {
51 12272 : if( keywords.exists("SERIAL") ) parseFlag("SERIAL",serial);
52 4916 : }
53 :
54 4916 : ActionWithVector::~ActionWithVector() {
55 4916 : if( action_to_do_before ) action_to_do_before->action_to_do_after=NULL;
56 4916 : }
57 :
58 273524 : void ActionWithVector::lockRequests() {
59 : ActionAtomistic::lockRequests();
60 : ActionWithArguments::lockRequests();
61 273524 : }
62 :
63 273524 : void ActionWithVector::unlockRequests() {
64 : ActionAtomistic::unlockRequests();
65 : ActionWithArguments::unlockRequests();
66 273524 : }
67 :
68 0 : void ActionWithVector::calculateNumericalDerivatives(ActionWithValue* av) {
69 0 : plumed_merror("cannot calculate numerical derivative for action " + getName() + " with label " + getLabel() );
70 : }
71 :
72 457632 : void ActionWithVector::clearDerivatives( const bool& force ) {
73 457632 : if( !force && actionInChain() ) return;
74 338953 : ActionWithValue::clearDerivatives();
75 338953 : if( action_to_do_after ) action_to_do_after->clearDerivatives( true );
76 : }
77 :
78 457632 : void ActionWithVector::clearInputForces( const bool& force ) {
79 457632 : if( !force && actionInChain() ) return;
80 338953 : ActionWithValue::clearInputForces();
81 338953 : if( action_to_do_after ) action_to_do_after->clearInputForces( true );
82 : }
83 :
84 46 : const ActionWithVector* ActionWithVector::getFirstActionInChain() const {
85 46 : if( !actionInChain() ) return this;
86 46 : return action_to_do_before->getFirstActionInChain();
87 : }
88 :
89 209535206 : ActionWithVector* ActionWithVector::getFirstActionInChain() {
90 209535206 : if( !actionInChain() ) return this;
91 198646985 : return action_to_do_before->getFirstActionInChain();
92 : }
93 :
94 456395 : void ActionWithVector::retrieveAtoms( const bool& force ) {
95 456395 : if( !force && actionInChain() || atomsWereRetrieved ) return;
96 337377 : ActionAtomistic::retrieveAtoms(); atomsWereRetrieved = !actionInChain();
97 337377 : if( action_to_do_after ) action_to_do_after->retrieveAtoms( true );
98 : }
99 :
100 46 : bool ActionWithVector::hasStoredArguments() const {
101 46 : std::string headstr=getFirstActionInChain()->getLabel();
102 100 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
103 68 : if( !getPntrToArgument(i)->ignoreStoredValue(headstr) ) return true;
104 : }
105 : return false;
106 : }
107 :
108 344 : bool ActionWithVector::argumentDependsOn( const std::string& headstr, ActionWithVector* faction, Value* thearg ) {
109 826 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
110 488 : if( this!=faction && thearg==getPntrToArgument(i) ) return true;
111 485 : ActionWithVector* av = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
112 485 : if( av && (av->getFirstActionInChain())->getLabel()==headstr ) {
113 180 : if( av->argumentDependsOn( headstr, faction, thearg ) ) return true;;
114 : }
115 : }
116 : return false;
117 : }
118 :
119 3401 : unsigned ActionWithVector::buildArgumentStore( const unsigned& argstart ) {
120 : // Don't use chains for grids
121 10627 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
122 7608 : if( getPntrToArgument(i)->isConstant() ) continue;
123 6707 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction());
124 6707 : if( !av || getPntrToArgument(i)->getRank()>0 && getPntrToArgument(i)->hasDerivatives() ) { done_in_chain=false; break; }
125 : }
126 3401 : if( done_in_chain ) {
127 : std::vector<std::string> alabels; std::vector<ActionWithVector*> f_actions;
128 8623 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
129 6200 : bool found=false; std::string mylab = (getPntrToArgument(i)->getPntrToAction())->getLabel();
130 12803 : for(unsigned j=0; j<alabels.size(); ++j) {
131 8749 : if( alabels[j]==mylab ) { found=true; break; }
132 : }
133 6200 : if( !found ) alabels.push_back( mylab );
134 :
135 : // If this is calculated in setup we never need to add to chain
136 6200 : if( getPntrToArgument(i)->isConstant() ) continue;
137 : // Find the chain we need to add this to from the arguments
138 5783 : ActionWithVector* av=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction()); plumed_assert( av );
139 5783 : found=false; ActionWithVector* myact = av->getFirstActionInChain();
140 5783 : if( getPntrToArgument(i)->getRank()==1 && getPntrToArgument(i)->storedata ) {
141 5540 : for(unsigned j=argstart; j<getNumberOfArguments(); ++j) {
142 5252 : ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getPntrToArgument(j)->getPntrToAction() );
143 5252 : if( !aarg || i==j ) continue;
144 14182 : for(unsigned k=0; k<aarg->getNumberOfArguments(); ++k) {
145 9263 : if( aarg->getPntrToArgument(k)==getPntrToArgument(i) ) { done_in_chain=false; return reallyBuildArgumentStore( argstart ); }
146 : }
147 : }
148 : }
149 6129 : for(unsigned j=0; j<f_actions.size(); ++j) {
150 3320 : if( f_actions[j]==myact ) { found=true; break; }
151 : }
152 5763 : if( !found ) {
153 2809 : if( f_actions.size()>0 ) {
154 336 : if( f_actions[0]->checkForDependency(myact) ) getPntrToArgument(i)->buildDataStore();
155 336 : if( myact->checkForDependency(f_actions[0]) ) error("cannot deal with arguments in this order. Try swapping argument order");
156 : }
157 2809 : if( !getPntrToArgument(i)->storedata && getPntrToArgument(i)->getRank()>0 ) f_actions.push_back( myact );
158 : }
159 : }
160 : // Now make sure that everything we need is in the chain
161 2423 : if( f_actions.size()>0 ) {
162 : // Check everything for later f_actions is done before f_actions[0]
163 2448 : for(unsigned i=1; i<f_actions.size(); ++i) {
164 122 : ActionWithArguments* aarg = dynamic_cast<ActionWithArguments*>( f_actions[i] );
165 122 : if( !aarg || aarg->getNumberOfArguments()==0 ) continue;
166 179 : for(unsigned j=0; j<aarg->getNumberOfArguments(); ++j) {
167 134 : if( (aarg->getPntrToArgument(j))->isConstant() ) continue ;
168 123 : bool found=false; std::string dep_argname = (aarg->getPntrToArgument(j))->getPntrToAction()->getLabel();
169 17501 : for(const auto & pp : plumed.getActionSet()) {
170 : Action* p(pp.get());
171 : // Check if this is the dependency
172 17501 : if( p->getLabel()==dep_argname ) { found=true; break; }
173 : // Check if this is the first of the arguments that will appear in this chain
174 17389 : else if( p->getLabel()==f_actions[0]->getLabel() ) break;
175 : }
176 123 : if( !found ) { done_in_chain=false; break; }
177 : }
178 : // Stop trying to add things in the chain if we cannot
179 56 : if( !done_in_chain ) return reallyBuildArgumentStore( argstart );
180 : }
181 2326 : std::vector<std::string> empty(1); empty[0] = f_actions[0]->getLabel();
182 2437 : for(unsigned i=1; i<f_actions.size(); ++i) f_actions[0]->addActionToChain( empty, f_actions[i] );
183 2326 : }
184 : // Now add this argument to the chain
185 : bool added=false;
186 2461 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) {
187 : // Add this function to jobs to do in recursive loop in previous action
188 2461 : if( getPntrToArgument(i)->getRank()>0 && !getPntrToArgument(i)->isConstant() ) {
189 2461 : ActionWithVector* av=dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
190 2461 : if( av && av->addActionToChain( alabels, this ) ) { added=true; break; }
191 : }
192 : }
193 2412 : plumed_massert(added, "could not add action " + getLabel() + " to chain of any of its arguments");
194 : // And get the number of derivatives
195 2412 : ActionWithVector* head=getFirstActionInChain();
196 2412 : unsigned nder=0; arg_deriv_starts.resize( getNumberOfArguments() );
197 8565 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
198 : // Ensures that we ignore the grid (first argument) if we are evaluating a function on a grid
199 6153 : if( i==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) continue ;
200 6152 : ActionWithVector* iaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(i)->getPntrToAction());
201 6152 : if( actionInChain() && !getPntrToArgument(i)->ignoreStoredValue(head->getLabel()) ) {
202 675 : arg_deriv_starts[i] = 0; head->getNumberOfStreamedDerivatives( arg_deriv_starts[i], getPntrToArgument(i) );
203 5477 : } else if( iaction && iaction->isInSubChain(nder) ) {
204 2255 : arg_deriv_starts[i] = nder;
205 : // Add the total number of derivatives that we have by this point in the chain to nder
206 2255 : if( iaction ) { nder=0; head->getNumberOfStreamedDerivatives( nder, getPntrToArgument(i) ); }
207 : } else {
208 : // Check if we have already found this action
209 : int k=-1;
210 3222 : if( iaction ) {
211 3222 : ActionWithVector* ider_action=iaction->getActionWithDerivatives( iaction );
212 3443 : for(unsigned j=0; j<i; ++j) {
213 1628 : if( j==0 && getName().find("EVALUATE_FUNCTION_FROM_GRID")!=std::string::npos ) continue ;
214 1627 : ActionWithVector* jaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(j)->getPntrToAction());
215 1627 : if( jaction->getActionWithDerivatives(jaction)==ider_action || jaction->checkForDependency(ider_action) ) { k=j; break; }
216 : }
217 3222 : if( k>=0 ) { arg_deriv_starts[i] = arg_deriv_starts[k]; continue; }
218 : }
219 :
220 1815 : if( i>0 ) {
221 : // This is a fudge so that inputs like this work:
222 : // c: CONTACT_MATRIX ATOMS=1-100
223 : // d: MATRIX_PRODUCT ARG=mat1,mat2
224 : // e: CUSTOM ARG=c,d
225 : // f: MATRIX_PRODUCT ARG=mat3,mat4
226 : // g: CUSTOM ARG=c,f
227 : // See symfunc rt-nbonds-q6 for an example
228 : // In this example when we set arg_deriv_starts[1] for f in g nder=number of derivatives of c
229 : // mder is equal to the number of derivatives by the time you get to f minus the number of derivatives for c
230 136 : unsigned mder=0;
231 136 : ActionWithVector* jaction=dynamic_cast<ActionWithVector*>(getPntrToArgument(i-1)->getPntrToAction());
232 136 : if( jaction->action_to_do_after && !(jaction->action_to_do_after)->getNumberOfStoredValues( getPntrToArgument(i-1), mder, i, getArguments() ) ) mder=0;
233 136 : if( mder>0 ) nder = nder + mder;
234 : }
235 :
236 1815 : arg_deriv_starts[i] = nder;
237 : // Add the total number of derivatives that we have by this point in the chain to nder
238 1815 : if( iaction ) {
239 1815 : nder=0;
240 1815 : if( (getPntrToArgument(i)->getPntrToAction())->getName().find("DIFFERENCE")!=std::string::npos ) {
241 0 : ActionWithArguments* aarg=dynamic_cast<ActionWithArguments*>( getPntrToArgument(i)->getPntrToAction() );
242 0 : plumed_assert( aarg && aarg->getNumberOfArguments()==2 );
243 0 : head->getNumberOfStreamedDerivatives( nder, aarg->getPntrToArgument(0) );
244 0 : nder += (aarg->getPntrToArgument(1))->getNumberOfValues();
245 1815 : } else head->getNumberOfStreamedDerivatives( nder, getPntrToArgument(i) );
246 : }
247 : }
248 : }
249 2412 : nder=0; head->getNumberOfStreamedDerivatives( nder, NULL );
250 2412 : return nder;
251 2443 : }
252 958 : return reallyBuildArgumentStore( argstart );
253 : }
254 :
255 989 : unsigned ActionWithVector::reallyBuildArgumentStore( const unsigned& argstart ) {
256 2703 : for(unsigned i=argstart; i<getNumberOfArguments(); ++i) { if( getPntrToArgument(i)->getRank()>0 ) getPntrToArgument(i)->buildDataStore(); }
257 989 : unsigned nder=0; arg_deriv_starts.resize( getNumberOfArguments() );
258 2703 : for(unsigned i=0; i<getNumberOfArguments(); ++i) { arg_deriv_starts[i] = nder; nder += getPntrToArgument(i)->getNumberOfValues(); }
259 989 : return nder;
260 : }
261 :
262 17691 : ActionWithVector* ActionWithVector::getActionWithDerivatives( ActionWithVector* depaction ) {
263 17691 : if( depaction==this || depaction->checkForDependency(this) ) {
264 11635 : if( getNumberOfAtoms()>0 ) return this;
265 6786 : std::string c=getFirstActionInChain()->getLabel();
266 42648 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
267 36322 : if( !getPntrToArgument(i)->ignoreStoredValue(c) && !getPntrToArgument(i)->isConstant() ) return this;
268 : }
269 : }
270 12842 : plumed_assert( action_to_do_before );
271 12842 : return action_to_do_before->getActionWithDerivatives(depaction);
272 : }
273 :
274 11762 : bool ActionWithVector::addActionToChain( const std::vector<std::string>& alabels, ActionWithVector* act ) {
275 11762 : if( action_to_do_after ) { bool state=action_to_do_after->addActionToChain( alabels, act ); return state; }
276 : // Check action is not already in chain
277 2572 : std::vector<std::string> mylabels; getFirstActionInChain()->getAllActionLabelsInChain( mylabels );
278 21211 : for(unsigned i=0; i<mylabels.size(); ++i) {
279 18639 : if( act->getLabel()==mylabels[i] ) return true;
280 : }
281 :
282 : // Check that everything that is required has been calculated
283 6756 : for(unsigned i=0; i<alabels.size(); ++i) {
284 : bool found=false;
285 20767 : for(unsigned j=0; j<mylabels.size(); ++j) {
286 20037 : if( alabels[i]==mylabels[j] ) { found=true; break; }
287 : }
288 4233 : if( !found ) {
289 730 : ActionWithValue* av=plumed.getActionSet().selectWithLabel<ActionWithValue*>( alabels[i] );
290 730 : plumed_massert( av, "could not cast " + alabels[i] ); bool storingall=true;
291 1522 : for(int j=0; j<av->getNumberOfComponents(); ++j) {
292 792 : if( !(av->getPntrToComponent(j))->storedata ) storingall=false;
293 : }
294 730 : if( !storingall ) return false;
295 : }
296 : }
297 : // This checks that there is nothing that will cause problems in the chain
298 2523 : mylabels.resize(0); getFirstActionInChain()->getAllActionLabelsInChain( mylabels );
299 20709 : for(unsigned i=0; i<mylabels.size(); ++i) {
300 18186 : ActionWithVector* av1=plumed.getActionSet().selectWithLabel<ActionWithVector*>( mylabels[i] );
301 176648 : for(unsigned j=0; j<i; ++j) {
302 158462 : ActionWithVector* av2=plumed.getActionSet().selectWithLabel<ActionWithVector*>( mylabels[j] );
303 158462 : if( !av1->canBeAfterInChain( av2 ) ) error("must calculate " + mylabels[j] + " before " + mylabels[i] );
304 : }
305 : }
306 2523 : action_to_do_after=act; act->action_to_do_before=this; updateTaskListReductionStatus();
307 2523 : ActionWithVector* head = getFirstActionInChain();
308 2523 : head->broadcastThatTasksAreReduced( head ); head->finishChainBuild( act );
309 : return true;
310 2572 : }
311 :
312 2674 : void ActionWithVector::updateTaskListReductionStatus() {
313 2674 : ActionWithVector* head = getFirstActionInChain();
314 2674 : std::vector<ActionWithVector*> task_reducing_actions; head->canReduceTasks( task_reducing_actions );
315 2674 : if( task_reducing_actions.size()>0 ) head->reduce_tasks=true;
316 2674 : }
317 :
318 3843518 : void ActionWithVector::broadcastThatTasksAreReduced( ActionWithVector* aselect ) {
319 3843518 : std::string c=getFirstActionInChain()->getLabel();
320 18085391 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
321 14241873 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
322 71210 : ActionWithVector* av = dynamic_cast<ActionWithVector*>( getPntrToArgument(i)->getPntrToAction() );
323 71210 : if( av ) {
324 : bool found=false;
325 12111 : ActionWithVector* av_head = av->getFirstActionInChain();
326 17731 : for(unsigned i=0; i<av_head->task_control_list.size(); ++i) {
327 17268 : if( aselect==av_head->task_control_list[i] ) { found=true; break; }
328 : }
329 12111 : if( !found ) av_head->task_control_list.insert( av_head->task_control_list.begin(), aselect );
330 :
331 12111 : av_head->reduce_tasks=true; av_head->updateTaskReductionFlag( av_head->reduce_tasks );
332 : }
333 : }
334 : }
335 3843518 : if( action_to_do_after ) action_to_do_after->broadcastThatTasksAreReduced( aselect );
336 3843518 : }
337 :
338 229464 : void ActionWithVector::updateTaskReductionFlag( bool& head_reduce_tasks ) {
339 229464 : if( actionInChain() ) {
340 217353 : plumed_assert( task_control_list.size()==0 );
341 : } else {
342 30767 : for(unsigned i=0; i<task_control_list.size(); ++i) {
343 18656 : if( !(task_control_list[i]->getFirstActionInChain())->reduce_tasks ) head_reduce_tasks=false;
344 : }
345 : }
346 229464 : broadcastThatTasksAreReduced( getFirstActionInChain() );
347 229464 : if( action_to_do_after ) action_to_do_after->updateTaskReductionFlag( head_reduce_tasks );
348 229464 : }
349 :
350 20916 : void ActionWithVector::canReduceTasks( std::vector<ActionWithVector*>& task_reducing_actions ) {
351 20916 : areAllTasksRequired( task_reducing_actions );
352 20916 : if( action_to_do_after ) action_to_do_after->canReduceTasks( task_reducing_actions );
353 20916 : }
354 :
355 3544 : void ActionWithVector::finishChainBuild( ActionWithVector* act ) {
356 3544 : if( action_to_do_after ) action_to_do_after->finishChainBuild( act );
357 3544 : }
358 :
359 36905 : void ActionWithVector::getAllActionLabelsInChain( std::vector<std::string>& mylabels ) const {
360 : bool found = false ;
361 359315 : for(unsigned i=0; i<mylabels.size(); ++i) {
362 322410 : if( getLabel()==mylabels[i] ) { found=true; }
363 : }
364 36905 : if( !found ) mylabels.push_back( getLabel() );
365 36905 : if( action_to_do_after ) action_to_do_after->getAllActionLabelsInChain( mylabels );
366 36905 : }
367 :
368 6793058 : void ActionWithVector::taskIsActive( const unsigned& current, int& flag ) const {
369 6793058 : if( isActive() ) flag = checkTaskStatus( current, flag );
370 6793058 : if( flag<=0 && action_to_do_after ) action_to_do_after->taskIsActive( current, flag );
371 6793058 : }
372 :
373 1309 : void ActionWithVector::getAdditionalTasksRequired( ActionWithVector* action, std::vector<unsigned>& atasks ) {
374 1341 : for(unsigned i=0; i<task_control_list.size(); ++i ) task_control_list[i]->getAdditionalTasksRequired( action, atasks );
375 1309 : }
376 :
377 274530 : void ActionWithVector::prepare() {
378 274530 : active_tasks.resize(0); atomsWereRetrieved=false;
379 274530 : }
380 :
381 144849 : std::vector<unsigned>& ActionWithVector::getListOfActiveTasks( ActionWithVector* action ) {
382 144849 : if( active_tasks.size()>0 ) return active_tasks;
383 132859 : unsigned ntasks=0; getNumberOfTasks( ntasks );
384 :
385 132859 : unsigned stride=comm.Get_size();
386 132859 : unsigned rank=comm.Get_rank();
387 132859 : if(serial) { stride=1; rank=0; }
388 :
389 : // Get number of threads for OpenMP
390 132859 : unsigned nt=OpenMP::getNumThreads();
391 132859 : if( nt*stride*10>ntasks ) nt=ntasks/stride/10;
392 : if( nt==0 ) nt=1;
393 :
394 132859 : if( !never_reduce_tasks && reduce_tasks ) {
395 1493 : if( task_control_list.size()>0 ) {
396 : // Get the list of tasks that are active in the action that uses the output of this action
397 238 : for(unsigned i=0; i<task_control_list.size(); ++i) {
398 130 : task_control_list[i]->retrieveAtoms();
399 130 : active_tasks = task_control_list[i]->getListOfActiveTasks( action );
400 : }
401 : // Now work out else we need from here to calculate the later action
402 108 : getAdditionalTasksRequired( action, active_tasks );
403 : } else {
404 1385 : std::vector<int> taskFlags( ntasks, -1 );
405 :
406 1385 : #pragma omp parallel num_threads(nt)
407 : {
408 : #pragma omp for nowait
409 : for(unsigned i=rank; i<ntasks; i+=stride ) {
410 : taskIsActive( i, taskFlags[i] );
411 : }
412 : }
413 2492457 : for(unsigned i=0; i<ntasks; ++i) taskFlags[i] = std::abs( taskFlags[i] );
414 1385 : if( !serial ) comm.Sum( taskFlags );
415 :
416 : unsigned nt=0;
417 2492457 : for(unsigned i=0; i<ntasks; ++i) {
418 2491072 : if( taskFlags[i]>=stride ) nt++;
419 : }
420 1385 : active_tasks.resize(nt); nt=0;
421 2492457 : for(unsigned i=0; i<ntasks; ++i) {
422 2491072 : if( taskFlags[i]>=stride ) { active_tasks[nt]=i; nt++; }
423 : }
424 1385 : getAdditionalTasksRequired( this, active_tasks );
425 : }
426 : } else {
427 131366 : active_tasks.resize( ntasks );
428 5116683 : for(unsigned i=0; i<ntasks; ++i) active_tasks[i]=i;
429 : }
430 132859 : return active_tasks;
431 : }
432 :
433 133004 : void ActionWithVector::runAllTasks() {
434 : // Skip this if this is done elsewhere
435 133004 : if( action_to_do_before ) return;
436 :
437 132969 : unsigned stride=comm.Get_size();
438 132969 : unsigned rank=comm.Get_rank();
439 132969 : if(serial) { stride=1; rank=0; }
440 :
441 : // Get the list of active tasks
442 132969 : std::vector<unsigned> & partialTaskList( getListOfActiveTasks( this ) );
443 132969 : unsigned nactive_tasks=partialTaskList.size();
444 :
445 : // Get number of threads for OpenMP
446 132969 : unsigned nt=OpenMP::getNumThreads();
447 132969 : if( nt*stride*10>nactive_tasks ) nt=nactive_tasks/stride/10;
448 132969 : if( nt==0 ) nt=1;
449 :
450 : // Now do all preparations required to run all the tasks
451 : // prepareForTaskLoop();
452 :
453 : // Get the total number of streamed quantities that we need
454 132969 : unsigned nquants=0, nmatrices=0, maxcol=0, nbooks=0;
455 132969 : getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol, nbooks );
456 : // Get size for buffer
457 132969 : unsigned bufsize=0; getSizeOfBuffer( nactive_tasks, bufsize );
458 132969 : if( buffer.size()!=bufsize ) buffer.resize( bufsize );
459 : // Clear buffer
460 132969 : buffer.assign( buffer.size(), 0.0 );
461 :
462 : // Recover the number of derivatives we require
463 132969 : unsigned nderivatives = 0; bool gridsInStream=checkForGrids(nderivatives);
464 132969 : if( !doNotCalculateDerivatives() && !gridsInStream ) getNumberOfStreamedDerivatives( nderivatives, NULL );
465 :
466 132969 : #pragma omp parallel num_threads(nt)
467 : {
468 : std::vector<double> omp_buffer;
469 : if( nt>1 ) omp_buffer.resize( bufsize, 0.0 );
470 : MultiValue myvals( nquants, nderivatives, nmatrices, maxcol, nbooks );
471 : myvals.clearAll();
472 :
473 : #pragma omp for nowait
474 : for(unsigned i=rank; i<nactive_tasks; i+=stride) {
475 : // Calculate the stuff in the loop for this action
476 : runTask( partialTaskList[i], myvals );
477 :
478 : // Now transfer the data to the actions that accumulate values from the calculated quantities
479 : if( nt>1 ) gatherAccumulators( partialTaskList[i], myvals, omp_buffer );
480 : else gatherAccumulators( partialTaskList[i], myvals, buffer );
481 :
482 : // Clear the value
483 : myvals.clearAll();
484 : }
485 : #pragma omp critical
486 : gatherThreads( nt, bufsize, omp_buffer, buffer, myvals );
487 : }
488 :
489 : // MPI Gather everything
490 132969 : if( !serial && buffer.size()>0 ) gatherProcesses( buffer );
491 132969 : finishComputations( buffer );
492 : }
493 :
494 209939 : void ActionWithVector::gatherThreads( const unsigned& nt, const unsigned& bufsize, const std::vector<double>& omp_buffer, std::vector<double>& buffer, MultiValue& myvals ) {
495 100317013 : if( nt>1 ) for(unsigned i=0; i<bufsize; ++i) buffer[i]+=omp_buffer[i];
496 209939 : }
497 :
498 132962 : void ActionWithVector::gatherProcesses( std::vector<double>& buffer ) {
499 132962 : comm.Sum( buffer );
500 132962 : }
501 :
502 317683 : bool ActionWithVector::checkForGrids( unsigned& nder ) const {
503 650581 : for(int i=0; i<getNumberOfComponents(); ++i) {
504 341631 : if( getConstPntrToComponent(i)->getRank()>0 && getConstPntrToComponent(i)->hasDerivatives() ) {
505 8733 : nder=getConstPntrToComponent(i)->getNumberOfGridDerivatives(); return true;
506 : }
507 : }
508 308950 : if( action_to_do_after ) return action_to_do_after->checkForGrids(nder);
509 : return false;
510 : }
511 :
512 314964 : void ActionWithVector::getNumberOfTasks( unsigned& ntasks ) {
513 314964 : if( ntasks==0 ) {
514 130206 : if( getNumberOfArguments()==1 && getNumberOfComponents()==1 && getPntrToComponent(0)->getRank()==0 ) {
515 2285 : if( !getPntrToArgument(0)->hasDerivatives() && getPntrToArgument(0)->getRank()==2 ) ntasks = getPntrToArgument(0)->getShape()[0];
516 2268 : else ntasks = getPntrToArgument(0)->getNumberOfValues();
517 : } else {
518 127921 : plumed_assert( getNumberOfComponents()>0 && getPntrToComponent(0)->getRank()>0 );
519 127921 : if( getPntrToComponent(0)->hasDerivatives() ) ntasks = getPntrToComponent(0)->getNumberOfValues();
520 121860 : else ntasks = getPntrToComponent(0)->getShape()[0];
521 : }
522 : }
523 653872 : for(int i=0; i<getNumberOfComponents(); ++i) {
524 338908 : if( getPntrToComponent(i)->getRank()==0 ) {
525 65937 : if( getNumberOfArguments()!=1 ) error("mismatched numbers of tasks in streamed quantities");
526 65937 : if( getPntrToArgument(0)->hasDerivatives() && ntasks!=getPntrToArgument(0)->getNumberOfValues() ) error("mismatched numbers of tasks in streamed quantities");
527 65937 : else if ( !getPntrToArgument(0)->hasDerivatives() && ntasks!=getPntrToArgument(0)->getShape()[0] ) error("mismatched numbers of tasks in streamed quantities");
528 272971 : } else if( getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getNumberOfValues() ) error("mismatched numbers of tasks in streamed quantities");
529 272971 : else if( !getPntrToComponent(i)->hasDerivatives() && ntasks!=getPntrToComponent(i)->getShape()[0] ) error("mismatched numbers of tasks in streamed quantities");
530 : }
531 314964 : if( action_to_do_after ) action_to_do_after->getNumberOfTasks( ntasks );
532 314964 : }
533 :
534 421278 : void ActionWithVector::setupStreamedComponents( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
535 1075219 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
536 653941 : if( !getPntrToArgument(i)->ignoreStoredValue(headstr) ) { getPntrToArgument(i)->streampos=nquants; nquants++; }
537 : }
538 881595 : for(int i=0; i<getNumberOfComponents(); ++i) { getPntrToComponent(i)->streampos=nquants; nquants++; }
539 421278 : }
540 :
541 421278 : void ActionWithVector::getNumberOfStreamedQuantities( const std::string& headstr, unsigned& nquants, unsigned& nmat, unsigned& maxcol, unsigned& nbookeeping ) {
542 421278 : setupStreamedComponents( headstr, nquants, nmat, maxcol, nbookeeping );
543 421278 : if( action_to_do_after ) action_to_do_after->getNumberOfStreamedQuantities( headstr, nquants, nmat, maxcol, nbookeeping );
544 421278 : }
545 :
546 317683 : void ActionWithVector::getSizeOfBuffer( const unsigned& nactive_tasks, unsigned& bufsize ) {
547 659314 : for(int i=0; i<getNumberOfComponents(); ++i) { getPntrToComponent(i)->bufstart=bufsize; bufsize += getPntrToComponent(i)->data.size(); }
548 317683 : if( action_to_do_after ) action_to_do_after->getSizeOfBuffer( nactive_tasks, bufsize );
549 317683 : }
550 :
551 424411 : void ActionWithVector::getNumberOfStreamedDerivatives( unsigned& nderivatives, Value* stopat ) {
552 424411 : std::string c=getFirstActionInChain()->getLabel();
553 1144906 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
554 721170 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
555 385865 : if( getPntrToArgument(i)==stopat ) return;
556 385190 : nderivatives += getPntrToArgument(i)->getNumberOfValues();
557 : }
558 : }
559 423736 : if( getNumberOfAtoms()>0 ) nderivatives += 3*getNumberOfAtoms() + 9;
560 : // Don't do the whole chain if we have been told to stop early
561 423736 : if( stopat && stopat->getPntrToAction()==this ) return;
562 :
563 419562 : if( action_to_do_after ) action_to_do_after->getNumberOfStreamedDerivatives( nderivatives, stopat );
564 : }
565 :
566 136 : bool ActionWithVector::getNumberOfStoredValues( Value* startat, unsigned& nvals, const unsigned& astart, const std::vector<Value*>& stopat ) {
567 181 : for(unsigned j=astart; j<stopat.size(); ++j) {
568 151 : if( stopat[j] && (stopat[j]->getPntrToAction()==this || (stopat[j]->getPntrToAction())->checkForDependency(this)) ) return true;
569 : }
570 :
571 30 : std::string c=getFirstActionInChain()->getLabel();
572 77 : for(unsigned i=0; i<getNumberOfArguments(); ++i) {
573 47 : if( !getPntrToArgument(i)->ignoreStoredValue(c) ) {
574 4 : for(unsigned j=astart; j<stopat.size(); ++j) {
575 2 : if( getPntrToArgument(i)==stopat[j] ) return true;
576 : }
577 2 : nvals += getPntrToArgument(i)->getNumberOfValues();
578 : }
579 : }
580 30 : if( startat->getPntrToAction()!=this && getNumberOfAtoms()>0 ) return false;
581 :
582 30 : if( action_to_do_after ) return action_to_do_after->getNumberOfStoredValues( startat, nvals, astart, stopat );
583 : return false;
584 : }
585 :
586 11340462 : void ActionWithVector::runTask( const unsigned& current, MultiValue& myvals ) const {
587 11340462 : if( isActive() ) {
588 10224413 : myvals.setTaskIndex(current); myvals.vector_call=true; performTask( current, myvals );
589 : }
590 11340462 : if( action_to_do_after ) action_to_do_after->runTask( current, myvals );
591 11340462 : }
592 :
593 9239280 : void ActionWithVector::gatherAccumulators( const unsigned& taskCode, const MultiValue& myvals, std::vector<double>& buffer ) const {
594 9239280 : if( isActive() ) {
595 18562303 : for(int i=0; i<getNumberOfComponents(); ++i) {
596 10284069 : unsigned bufstart = getConstPntrToComponent(i)->bufstart;
597 : // This looks after storing of scalars that are summed from vectors/matrices
598 10284069 : if( getConstPntrToComponent(i)->getRank()==0 ) {
599 : plumed_dbg_massert( bufstart<buffer.size(), "problem in " + getLabel() );
600 1146054 : unsigned sind = getConstPntrToComponent(i)->streampos; buffer[bufstart] += myvals.get(sind);
601 1146054 : if( getConstPntrToComponent(i)->hasDerivatives() ) {
602 26904718 : for(unsigned k=0; k<myvals.getNumberActive(sind); ++k) {
603 25758664 : unsigned kindex = myvals.getActiveIndex(sind,k);
604 : plumed_dbg_massert( bufstart+1+kindex<buffer.size(), "problem in " + getLabel() );
605 25758664 : buffer[bufstart + 1 + kindex] += myvals.getDerivative(sind,kindex);
606 : }
607 : }
608 : // This looks after storing of vectors
609 9138015 : } else if( getConstPntrToComponent(i)->storedata ) gatherStoredValue( i, taskCode, myvals, bufstart, buffer );
610 : }
611 : }
612 9239280 : if( action_to_do_after ) action_to_do_after->gatherAccumulators( taskCode, myvals, buffer );
613 9239280 : }
614 :
615 3571815 : void ActionWithVector::gatherStoredValue( const unsigned& valindex, const unsigned& taskCode, const MultiValue& myvals,
616 : const unsigned& bufstart, std::vector<double>& buffer ) const {
617 : plumed_dbg_assert( getConstPntrToComponent(valindex)->getRank()==1 && !getConstPntrToComponent(valindex)->hasDeriv );
618 3571815 : unsigned vindex = getConstPntrToComponent(valindex)->bufstart + taskCode; plumed_dbg_massert( vindex<buffer.size(), "failing in " + getLabel() );
619 3571815 : buffer[vindex] += myvals.get(getConstPntrToComponent(valindex)->streampos);
620 3571815 : }
621 :
622 317683 : void ActionWithVector::finishComputations( const std::vector<double>& buf ) {
623 317683 : if( isActive() ) {
624 527282 : for(int i=0; i<getNumberOfComponents(); ++i) {
625 : // This gathers vectors and grids at the end of the calculation
626 275615 : unsigned bufstart = getPntrToComponent(i)->bufstart;
627 275615 : getPntrToComponent(i)->data.assign( getPntrToComponent(i)->data.size(), 0 );
628 275615 : if( (getPntrToComponent(i)->getRank()>0 && getPntrToComponent(i)->hasDerivatives()) || getPntrToComponent(i)->storedata ) {
629 130918 : unsigned sz_v = getPntrToComponent(i)->data.size();
630 46583551 : for(unsigned j=0; j<sz_v; ++j) {
631 : plumed_dbg_assert( bufstart+j<buf.size() );
632 46452633 : getPntrToComponent(i)->add( j, buf[bufstart+j] );
633 : }
634 : // Make sure single values are set
635 144697 : } else if( getPntrToComponent(i)->getRank()==0 ) getPntrToComponent(i)->set( buf[bufstart] );
636 : // This gathers derivatives of scalars
637 275615 : if( !doNotCalculateDerivatives() && getPntrToComponent(i)->hasDeriv && getPntrToComponent(i)->getRank()==0 ) {
638 16229167 : for(unsigned j=0; j<getPntrToComponent(i)->getNumberOfDerivatives(); ++j) getPntrToComponent(i)->setDerivative( j, buf[bufstart+1+j] );
639 : }
640 : }
641 : }
642 317683 : if( action_to_do_after ) action_to_do_after->finishComputations( buf );
643 317683 : }
644 :
645 292884 : bool ActionWithVector::checkChainForNonScalarForces() const {
646 523180 : for(unsigned i=0; i<getNumberOfComponents(); ++i) {
647 305133 : if( getConstPntrToComponent(i)->getRank()>0 && getConstPntrToComponent(i)->forcesWereAdded() ) return true;
648 : }
649 218047 : if( action_to_do_after ) return action_to_do_after->checkChainForNonScalarForces();
650 : return false;
651 : }
652 :
653 254422 : bool ActionWithVector::checkForForces() {
654 254422 : if( getPntrToComponent(0)->getRank()==0 ) return ActionWithValue::checkForForces();
655 202650 : else if( actionInChain() ) return false;
656 :
657 : // Check if there are any forces
658 125264 : if( !checkChainForNonScalarForces() ) return false;
659 :
660 : // Setup MPI parallel loop
661 74837 : unsigned stride=comm.Get_size();
662 74837 : unsigned rank=comm.Get_rank();
663 74837 : if(serial) { stride=1; rank=0; }
664 :
665 : // Get the number of tasks
666 74837 : std::vector<unsigned> force_tasks; getForceTasks( force_tasks );
667 74837 : Tools::removeDuplicates(force_tasks); unsigned nf_tasks=force_tasks.size();
668 :
669 : // Get number of threads for OpenMP
670 74837 : unsigned nt=OpenMP::getNumThreads();
671 74837 : if( nt*stride*10>nf_tasks ) nt=nf_tasks/stride/10;
672 : if( nt==0 ) nt=1;
673 :
674 : // Now determine how big the multivalue needs to be
675 74837 : unsigned nquants=0, nmatrices=0, maxcol=0, nbooks=0;
676 74837 : getNumberOfStreamedQuantities( getLabel(), nquants, nmatrices, maxcol, nbooks );
677 : // Recover the number of derivatives we require (this should be equal to the number of forces)
678 74837 : unsigned nderiv=0; getNumberOfStreamedDerivatives( nderiv, NULL );
679 74837 : if( forcesForApply.size()!=nderiv ) forcesForApply.resize( nderiv );
680 : // Clear force buffer
681 74837 : forcesForApply.assign( forcesForApply.size(), 0.0 );
682 :
683 74837 : #pragma omp parallel num_threads(nt)
684 : {
685 : std::vector<double> omp_forces;
686 : if( nt>1 ) omp_forces.resize( forcesForApply.size(), 0.0 );
687 : MultiValue myvals( nquants, nderiv, nmatrices, maxcol, nbooks );
688 : myvals.clearAll();
689 :
690 : #pragma omp for nowait
691 : for(unsigned i=rank; i<nf_tasks; i+=stride) {
692 : runTask( force_tasks[i], myvals );
693 :
694 : // Now get the forces
695 : if( nt>1 ) gatherForces( force_tasks[i], myvals, omp_forces );
696 : else gatherForces( force_tasks[i], myvals, forcesForApply );
697 :
698 : myvals.clearAll();
699 : }
700 : #pragma omp critical
701 : if(nt>1) for(unsigned i=0; i<forcesForApply.size(); ++i) forcesForApply[i]+=omp_forces[i];
702 : }
703 : // MPI Gather on forces
704 74837 : if( !serial ) comm.Sum( forcesForApply );
705 : return true;
706 : }
707 :
708 2039456 : bool ActionWithVector::checkComponentsForForce() const {
709 2777915 : for(unsigned i=0; i<values.size(); ++i) {
710 2546733 : if( getConstPntrToComponent(i)->forcesWereAdded() ) return true;
711 : }
712 : return false;
713 : }
714 :
715 5561312 : bool ActionWithVector::checkForTaskForce( const unsigned& itask, const Value* myval ) const {
716 : plumed_dbg_assert( !(myval->getRank()==2 && !myval->hasDerivatives()) );
717 5561312 : return fabs(myval->getForce(itask))>epsilon;
718 : }
719 :
720 90719 : void ActionWithVector::updateForceTasksFromValue( const Value* myval, std::vector<unsigned>& force_tasks ) const {
721 90719 : if( myval->getRank()>0 && myval->forcesWereAdded() ) {
722 85478 : unsigned nt = myval->getNumberOfValues();
723 85478 : if( !myval->hasDerivatives() ) nt = myval->getShape()[0];
724 6149611 : for(unsigned i=0; i<nt; ++i) {
725 6064133 : if( checkForTaskForce(i, myval) ) force_tasks.push_back( i );
726 : }
727 : }
728 90719 : }
729 :
730 103595 : void ActionWithVector::getForceTasks( std::vector<unsigned>& force_tasks ) const {
731 103595 : if( isActive() && checkComponentsForForce() ) {
732 171577 : for(unsigned k=0; k<values.size(); ++k) {
733 91329 : updateForceTasksFromValue( getConstPntrToComponent(k), force_tasks );
734 : }
735 : }
736 103595 : if( action_to_do_after ) action_to_do_after->getForceTasks( force_tasks );
737 103595 : }
738 :
739 1565796 : void ActionWithVector::gatherForcesOnStoredValue( const Value* myval, const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
740 : plumed_dbg_assert( myval->storedata && !(myval->getRank()==2 && !myval->hasDerivatives()) );
741 1565796 : double fforce = myval->getForce(itask);
742 : unsigned sspos = myval->getPositionInStream();
743 20912845 : for(unsigned j=0; j<myvals.getNumberActive(sspos); ++j) {
744 19347049 : unsigned jder=myvals.getActiveIndex(sspos, j); plumed_dbg_assert( jder<forces.size() );
745 19347049 : forces[jder] += fforce*myvals.getDerivative( sspos, jder );
746 : }
747 1565796 : }
748 :
749 2101182 : void ActionWithVector::gatherForces( const unsigned& itask, const MultiValue& myvals, std::vector<double>& forces ) const {
750 2101182 : if( isActive() && checkComponentsForForce() ) {
751 3587802 : for(unsigned k=0; k<getNumberOfComponents(); ++k) {
752 1859776 : const Value* myval=getConstPntrToComponent(k);
753 1859776 : if( myval->getRank()>0 && myval->forcesWereAdded() ) gatherForcesOnStoredValue( myval, itask, myvals, forces );
754 : }
755 : }
756 2101182 : if( action_to_do_after ) action_to_do_after->gatherForces( itask, myvals, forces );
757 2101182 : }
758 :
759 254422 : void ActionWithVector::apply() {
760 254422 : if( !checkForForces() ) return;
761 : // Find the top of the chain and add forces
762 117831 : unsigned ind=0; getFirstActionInChain()->addForcesToInput( getForcesToApply(), ind );
763 : }
764 :
765 186367 : void ActionWithVector::addForcesToInput( const std::vector<double>& forcesToApply, unsigned& ind ) {
766 186367 : if( ind>=forcesToApply.size() ) return;
767 136821 : addForcesOnArguments( 0, forcesToApply, ind, getFirstActionInChain()->getLabel() ); setForcesOnAtoms( forcesToApply, ind );
768 136821 : if( action_to_do_after ) action_to_do_after->addForcesToInput( forcesToApply, ind );
769 : }
770 :
771 : }
|
