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 "PlumedMain.h"
23 : #include "ActionAtomistic.h"
24 : #include "ActionPilot.h"
25 : #include "ActionForInterface.h"
26 : #include "ActionRegister.h"
27 : #include "ActionSet.h"
28 : #include "ActionWithValue.h"
29 : #include "ActionWithVirtualAtom.h"
30 : #include "ActionToGetData.h"
31 : #include "ActionToPutData.h"
32 : #include "CLToolMain.h"
33 : #include "ExchangePatterns.h"
34 : #include "GREX.h"
35 : #include "DomainDecomposition.h"
36 : #include "config/Config.h"
37 : #include "tools/Citations.h"
38 : #include "tools/Communicator.h"
39 : #include "tools/DLLoader.h"
40 : #include "tools/Exception.h"
41 : #include "tools/IFile.h"
42 : #include "tools/Log.h"
43 : #include "tools/OpenMP.h"
44 : #include "tools/Tools.h"
45 : #include "tools/Stopwatch.h"
46 : #include "tools/TypesafePtr.h"
47 : #include "lepton/Exception.h"
48 : #include "DataPassingTools.h"
49 : #include "small_vector/small_vector.h"
50 : #include <cstdlib>
51 : #include <cstdio>
52 : #include <cstring>
53 : #include <set>
54 : #include <exception>
55 : #include <stdexcept>
56 : #include <ios>
57 : #include <new>
58 : #include <typeinfo>
59 : #include <iostream>
60 : #include <algorithm>
61 : #include <system_error>
62 : #include <future>
63 : #include <memory>
64 : #include <functional>
65 : #include <regex>
66 : #include <any>
67 : #include <optional>
68 : #include <variant>
69 : #include <filesystem>
70 :
71 : namespace PLMD {
72 :
73 : /// Small utility just used in this file to throw arbitrary exceptions
74 61 : [[noreturn]] static void testThrow(const char* what) {
75 61 : auto words=Tools::getWords(what);
76 61 : plumed_assert(words.size()>0);
77 : #define __PLUMED_THROW_NOMSG(type) if(words[0]==#type) throw type()
78 : #define __PLUMED_THROW_MSG(type) if(words[0]==#type) throw type(what)
79 62 : __PLUMED_THROW_MSG(PLMD::ExceptionError);
80 61 : __PLUMED_THROW_MSG(PLMD::ExceptionDebug);
81 60 : __PLUMED_THROW_MSG(PLMD::Exception);
82 59 : __PLUMED_THROW_MSG(PLMD::lepton::Exception);
83 57 : __PLUMED_THROW_NOMSG(std::bad_exception);
84 56 : __PLUMED_THROW_NOMSG(std::bad_array_new_length);
85 55 : __PLUMED_THROW_NOMSG(std::bad_alloc);
86 54 : __PLUMED_THROW_NOMSG(std::bad_function_call);
87 53 : __PLUMED_THROW_NOMSG(std::bad_weak_ptr);
88 52 : __PLUMED_THROW_NOMSG(std::bad_cast);
89 51 : __PLUMED_THROW_NOMSG(std::bad_typeid);
90 50 : __PLUMED_THROW_NOMSG(std::bad_variant_access);
91 49 : __PLUMED_THROW_NOMSG(std::bad_optional_access);
92 48 : __PLUMED_THROW_NOMSG(std::bad_any_cast);
93 47 : __PLUMED_THROW_MSG(std::underflow_error);
94 46 : __PLUMED_THROW_MSG(std::overflow_error);
95 45 : __PLUMED_THROW_MSG(std::range_error);
96 44 : __PLUMED_THROW_MSG(std::runtime_error);
97 43 : __PLUMED_THROW_MSG(std::out_of_range);
98 42 : __PLUMED_THROW_MSG(std::length_error);
99 41 : __PLUMED_THROW_MSG(std::domain_error);
100 40 : __PLUMED_THROW_MSG(std::invalid_argument);
101 39 : __PLUMED_THROW_MSG(std::logic_error);
102 :
103 :
104 :
105 38 : if(words[0]=="std::system_error") {
106 4 : plumed_assert(words.size()>2);
107 : int error_code;
108 4 : Tools::convert(words[2],error_code);
109 4 : if(words[1]=="std::generic_category") throw std::system_error(error_code,std::generic_category(),what);
110 3 : if(words[1]=="std::system_category") throw std::system_error(error_code,std::system_category(),what);
111 2 : if(words[1]=="std::iostream_category") throw std::system_error(error_code,std::iostream_category(),what);
112 1 : if(words[1]=="std::future_category") throw std::system_error(error_code,std::future_category(),what);
113 : }
114 :
115 34 : if(words[0]=="std::filesystem::filesystem_error") {
116 : int error_code;
117 3 : plumed_assert(words.size()>2);
118 3 : Tools::convert(words[2],error_code);
119 : std::error_code x_error_code;
120 3 : if(words[1]=="std::generic_category") x_error_code=::std::error_code(error_code,::std::generic_category());
121 3 : if(words[1]=="std::system_category") x_error_code=::std::error_code(error_code,::std::system_category());
122 3 : if(words[1]=="std::iostream_category") x_error_code=::std::error_code(error_code,::std::iostream_category());
123 3 : if(words[1]=="std::future_category") x_error_code=::std::error_code(error_code,::std::future_category());
124 :
125 4 : if(words.size()<4) throw std::filesystem::filesystem_error(what,x_error_code);
126 4 : if(words.size()<5) throw std::filesystem::filesystem_error(what,std::filesystem::path(words[3]),x_error_code);
127 4 : throw std::filesystem::filesystem_error(what,std::filesystem::path(words[3]),std::filesystem::path(words[4]),x_error_code);
128 : }
129 :
130 : #define __PLUMED_THROW_REGEX(name) if(words[1]=="std::regex_constants::error_" #name) throw std::regex_error(std::regex_constants::error_ ##name)
131 31 : if(words[0]=="std::regex_error") {
132 13 : plumed_assert(words.size()>1);
133 13 : __PLUMED_THROW_REGEX(collate);
134 12 : __PLUMED_THROW_REGEX(ctype);
135 11 : __PLUMED_THROW_REGEX(escape);
136 10 : __PLUMED_THROW_REGEX(backref);
137 9 : __PLUMED_THROW_REGEX(brack);
138 8 : __PLUMED_THROW_REGEX(paren);
139 7 : __PLUMED_THROW_REGEX(brace);
140 6 : __PLUMED_THROW_REGEX(badbrace);
141 5 : __PLUMED_THROW_REGEX(range);
142 4 : __PLUMED_THROW_REGEX(space);
143 3 : __PLUMED_THROW_REGEX(badrepeat);
144 2 : __PLUMED_THROW_REGEX(complexity);
145 1 : __PLUMED_THROW_REGEX(stack);
146 : }
147 :
148 : #define __PLUMED_THROW_FUTURE(name) if(words[1]=="std::future_errc::" #name) throw std::future_error(::std::future_errc::name)
149 18 : if(words[0]=="std::future_error") {
150 4 : plumed_assert(words.size()>1);
151 4 : __PLUMED_THROW_FUTURE(broken_promise);
152 3 : __PLUMED_THROW_FUTURE(future_already_retrieved);
153 2 : __PLUMED_THROW_FUTURE(promise_already_satisfied);
154 1 : __PLUMED_THROW_FUTURE(no_state);
155 : }
156 :
157 14 : if(words[0]=="std::ios_base::failure") {
158 1 : int error_code=0;
159 1 : if(words.size()>2) Tools::convert(words[2],error_code);
160 1 : if(words.size()>1 && words[1]=="std::generic_category") throw std::ios_base::failure(what,std::error_code(error_code,std::generic_category()));
161 1 : if(words.size()>1 && words[1]=="std::system_category") throw std::ios_base::failure(what,std::error_code(error_code,std::system_category()));
162 1 : if(words.size()>1 && words[1]=="std::iostream_category") throw std::ios_base::failure(what,std::error_code(error_code,std::iostream_category()));
163 1 : if(words.size()>1 && words[1]=="std::future_category") throw std::ios_base::failure(what,std::error_code(error_code,std::future_category()));
164 2 : throw std::ios_base::failure(what);
165 : }
166 :
167 13 : if(words[0]=="int") {
168 1 : int value=0;
169 1 : if(words.size()>1) Tools::convert(words[1],value);
170 1 : throw value;
171 : }
172 :
173 12 : if(words[0]=="test_nested1") {
174 : try {
175 12 : throw Exception(std::string("inner ")+what);
176 6 : } catch(...) {
177 : try {
178 18 : std::throw_with_nested(Exception(std::string("middle ")+what));
179 6 : } catch(...) {
180 18 : std::throw_with_nested(Exception(std::string("outer ")+what));
181 6 : }
182 6 : }
183 : }
184 :
185 6 : if(words[0]=="test_nested2") {
186 : try {
187 3 : throw std::bad_alloc();
188 3 : } catch(...) {
189 : try {
190 9 : std::throw_with_nested(Exception(std::string("middle ")+what));
191 3 : } catch(...) {
192 9 : std::throw_with_nested(Exception(std::string("outer ")+what));
193 3 : }
194 3 : }
195 : }
196 :
197 3 : if(words[0]=="test_nested3") {
198 : try {
199 2 : throw "inner";
200 2 : } catch(...) {
201 : try {
202 6 : std::throw_with_nested(Exception(std::string("middle ")+what));
203 2 : } catch(...) {
204 6 : std::throw_with_nested(Exception(std::string("outer ")+what));
205 2 : }
206 2 : }
207 : }
208 :
209 2 : plumed_error() << "unknown exception " << what;
210 61 : }
211 :
212 : namespace {
213 : /// This is an internal tool used to count how many PlumedMain objects have been created
214 : /// and if they were correctly destroyed.
215 : /// When using debug options, it leads to a crash
216 : /// Otherwise, it just prints a message
217 : class CountInstances {
218 : std::atomic<int> counter{};
219 : // private constructor to avoid direct usage
220 5293 : CountInstances() noexcept {}
221 5293 : ~CountInstances() {
222 5293 : if(counter!=0) {
223 0 : std::cerr<<"WARNING: internal inconsistency in allocated PlumedMain instances (" <<counter<< ")\n";
224 0 : std::cerr<<"Might be a consequence of incorrectly paired plumed_create/plumed_finalize in the C interface\n";
225 0 : std::cerr<<"Or it could be due to incorrect calls to std::exit, without properly destroying all PlumedMain objects\n";
226 : #ifndef NDEBUG
227 : std::cerr<<"This is a debug build, so the warning will make PLUMED abort\n";
228 : std::abort();
229 : #endif
230 : }
231 5293 : }
232 1613396 : static CountInstances & instance() {
233 1613396 : static CountInstances counter;
234 1613396 : return counter;
235 : }
236 : public:
237 : /// Only access through these static functions
238 : /// The first call to increase() ensures the instance is constructed
239 : /// This should provide the correct construction and destruction order
240 : /// also in cases where the PlumedMain object is constructed in the
241 : /// constructor of a static object
242 806698 : static void increase() noexcept {
243 806698 : ++instance().counter;
244 806698 : }
245 : /// See increase()
246 806698 : static void decrease() noexcept {
247 806698 : --instance().counter;
248 806698 : }
249 : };
250 :
251 : }
252 :
253 :
254 806698 : PlumedMain::PlumedMain():
255 806698 : datoms_fwd(*this),
256 : // automatically write on log in destructor
257 806698 : stopwatch_fwd(log),
258 806698 : actionSet_fwd(*this),
259 806698 : passtools(DataPassingTools::create(sizeof(double)))
260 : {
261 806698 : passtools->usingNaturalUnits=false;
262 806698 : increaseReferenceCounter();
263 806698 : log.link(comm);
264 806698 : log.setLinePrefix("PLUMED: ");
265 : // this is at last so as to avoid inconsistencies if an exception is thrown
266 806698 : CountInstances::increase(); // noexcept
267 806698 : }
268 :
269 : // destructor needed to delete forward declarated objects
270 1612419 : PlumedMain::~PlumedMain() {
271 806698 : CountInstances::decrease();
272 4839211 : }
273 :
274 : /////////////////////////////////////////////////////////////
275 : // MAIN INTERPRETER
276 :
277 : #define CHECK_INIT(ini,word) plumed_assert(ini)<<"cmd(\"" << word << "\") should be only used after plumed initialization"
278 : #define CHECK_NOTINIT(ini,word) plumed_assert(!(ini))<<"cmd(\"" << word << "\") should be only used before plumed initialization"
279 : #define CHECK_NOTNULL(val,word) plumed_assert(val)<<"NULL pointer received in cmd(\"" << word << "\")"
280 :
281 :
282 1331107 : void PlumedMain::cmd(std::string_view word,const TypesafePtr & val) {
283 :
284 : // Enumerate all possible commands:
285 : enum {
286 : #include "PlumedMainEnum.inc"
287 : };
288 :
289 : // Static object (initialized once) containing the map of commands:
290 : const static Tools::FastStringUnorderedMap<int> word_map = {
291 : #include "PlumedMainMap.inc"
292 1331107 : };
293 :
294 : try {
295 :
296 1331107 : auto ss=stopwatch.startPause();
297 :
298 : gch::small_vector<std::string_view> words;
299 1331107 : Tools::getWordsSimple(words,word);
300 :
301 1331107 : unsigned nw=words.size();
302 1331107 : if(nw==0) {
303 : // do nothing
304 : } else {
305 : int iword=-1;
306 : const auto it=word_map.find(words[0]);
307 1331106 : if(it!=word_map.end()) iword=it->second;
308 :
309 1331099 : switch(iword) {
310 68481 : case cmd_setBox:
311 68481 : CHECK_INIT(initialized,word);
312 68481 : CHECK_NOTNULL(val,word);
313 68481 : setInputValue( "Box", 0, 1, val );
314 68480 : break;
315 72658 : case cmd_setPositions:
316 72658 : CHECK_INIT(initialized,word);
317 72658 : setInputValue("posx", 0, 3, val );
318 72658 : setInputValue("posy", 1, 3, val );
319 72658 : setInputValue("posz", 2, 3, val );
320 72658 : break;
321 72621 : case cmd_setMasses:
322 72621 : CHECK_INIT(initialized,word);
323 72621 : setInputValue("Masses", 0, 1, val );
324 72618 : break;
325 64546 : case cmd_setCharges:
326 64546 : CHECK_INIT(initialized,word);
327 64546 : setInputValue("Charges", 0, 1, val);
328 64546 : break;
329 1 : case cmd_setPositionsX:
330 1 : CHECK_INIT(initialized,word);
331 1 : setInputValue("posx", 0, 1, val);
332 1 : break;
333 1 : case cmd_setPositionsY:
334 1 : CHECK_INIT(initialized,word);
335 1 : setInputValue("posy", 0, 1, val);
336 1 : break;
337 1 : case cmd_setPositionsZ:
338 1 : CHECK_INIT(initialized,word);
339 1 : setInputValue("posz", 0, 1, val);
340 1 : break;
341 64795 : case cmd_setVirial:
342 64795 : CHECK_INIT(initialized,word);
343 64795 : CHECK_NOTNULL(val,word);
344 64795 : setInputForce("Box",val);
345 64777 : break;
346 7842 : case cmd_setEnergy:
347 7842 : CHECK_INIT(initialized,word);
348 7842 : CHECK_NOTNULL(val,word);
349 7842 : if( name_of_energy!="" ) setInputValue( name_of_energy, 0, 1, val );
350 : break;
351 72625 : case cmd_setForces:
352 72625 : CHECK_INIT(initialized,word);
353 72625 : setInputForce("posx",val);
354 72625 : setInputForce("posy",val);
355 72625 : setInputForce("posz",val);
356 72625 : break;
357 1 : case cmd_setForcesX:
358 1 : CHECK_INIT(initialized,word);
359 1 : setInputForce("posx",val);
360 1 : break;
361 1 : case cmd_setForcesY:
362 1 : CHECK_INIT(initialized,word);
363 1 : setInputForce("posy",val);
364 1 : break;
365 1 : case cmd_setForcesZ:
366 1 : CHECK_INIT(initialized,word);
367 1 : setInputForce("posz",val);
368 1 : break;
369 278008 : case cmd_calc:
370 278008 : CHECK_INIT(initialized,word);
371 278008 : calc();
372 : break;
373 99 : case cmd_prepareDependencies:
374 99 : CHECK_INIT(initialized,word);
375 99 : prepareDependencies();
376 : break;
377 84 : case cmd_shareData:
378 84 : CHECK_INIT(initialized,word);
379 84 : shareData();
380 : break;
381 6951 : case cmd_prepareCalc:
382 6951 : CHECK_INIT(initialized,word);
383 6951 : prepareCalc();
384 : break;
385 10 : case cmd_performCalc:
386 10 : CHECK_INIT(initialized,word);
387 10 : performCalc();
388 : break;
389 7011 : case cmd_performCalcNoUpdate:
390 7011 : CHECK_INIT(initialized,word);
391 7011 : performCalcNoUpdate();
392 : break;
393 10 : case cmd_performCalcNoForces:
394 10 : CHECK_INIT(initialized,word);
395 10 : performCalcNoForces();
396 : break;
397 79 : case cmd_update:
398 79 : CHECK_INIT(initialized,word);
399 79 : update();
400 : break;
401 16009 : case cmd_setStep:
402 16009 : CHECK_INIT(initialized,word);
403 16012 : CHECK_NOTNULL(val,word);
404 16006 : step=val.get<int>();
405 16005 : startStep();
406 : break;
407 0 : case cmd_setStepLong:
408 0 : CHECK_INIT(initialized,word);
409 0 : CHECK_NOTNULL(val,word);
410 0 : step=val.get<long int>();
411 0 : startStep();
412 : break;
413 268951 : case cmd_setStepLongLong:
414 268951 : CHECK_INIT(initialized,word);
415 268951 : CHECK_NOTNULL(val,word);
416 268951 : step=val.get<long long int>();
417 268951 : startStep();
418 : break;
419 23456 : case cmd_setValue:
420 : {
421 23456 : CHECK_INIT(initialized,words[0]); plumed_assert(nw==2); setInputValue( std::string(words[1]), 0, 1, val );
422 : }
423 23456 : break;
424 : /* ADDED WITH API=7 */
425 0 : case cmd_setValueForces:
426 : {
427 0 : CHECK_INIT(initialized,words[0]); plumed_assert(nw==2); setInputForce( std::string(words[1]), val );
428 : }
429 0 : break;
430 : // words used less frequently:
431 1140 : case cmd_setAtomsNlocal:
432 1140 : CHECK_INIT(initialized,word);
433 1140 : CHECK_NOTNULL(val,word);
434 4582 : for(const auto & pp : inputs ) {
435 3442 : plumed_assert(pp);
436 3442 : DomainDecomposition* dd=pp->castToDomainDecomposition();
437 3442 : if( dd ) dd->setAtomsNlocal(val.get<int>());
438 : }
439 : break;
440 988 : case cmd_setAtomsGatindex:
441 988 : CHECK_INIT(initialized,word);
442 3974 : for(const auto & pp : inputs ) {
443 2986 : plumed_assert(pp);
444 2986 : DomainDecomposition* dd=pp->castToDomainDecomposition();
445 2986 : if( dd ) dd->setAtomsGatindex(val,false);
446 : }
447 : break;
448 0 : case cmd_setAtomsFGatindex:
449 0 : CHECK_INIT(initialized,word);
450 0 : for(const auto & pp : inputs ) {
451 0 : plumed_assert(pp);
452 0 : DomainDecomposition* dd=pp->castToDomainDecomposition();
453 0 : if( dd ) dd->setAtomsGatindex(val,false);
454 : }
455 : break;
456 152 : case cmd_setAtomsContiguous:
457 152 : CHECK_INIT(initialized,word);
458 152 : CHECK_NOTNULL(val,word);
459 608 : for(const auto & pp : inputs ) {
460 456 : plumed_assert(pp);
461 456 : DomainDecomposition* dd=pp->castToDomainDecomposition();
462 456 : if( dd ) dd->setAtomsContiguous(val.get<int>());
463 : }
464 : break;
465 116 : case cmd_createFullList:
466 116 : CHECK_INIT(initialized,word);
467 116 : CHECK_NOTNULL(val,word);
468 553 : for(const auto & pp : inputs ) {
469 437 : plumed_assert(pp);
470 437 : DomainDecomposition* dd=pp->castToDomainDecomposition();
471 437 : if( dd ) dd->createFullList(val);
472 : }
473 : break;
474 116 : case cmd_getFullList:
475 : {
476 116 : CHECK_INIT(initialized,word);
477 116 : CHECK_NOTNULL(val,word);
478 : unsigned nlists=0;
479 553 : for(const auto & pp : inputs ) {
480 437 : plumed_assert(pp);
481 437 : DomainDecomposition* dd=pp->castToDomainDecomposition();
482 437 : if( dd ) { dd->getFullList(val); nlists++; }
483 : }
484 116 : plumed_assert( nlists==1 );
485 : }
486 : break;
487 116 : case cmd_clearFullList:
488 116 : CHECK_INIT(initialized,word);
489 553 : for(const auto & pp : inputs ) {
490 437 : plumed_assert(pp);
491 437 : DomainDecomposition* dd=pp->castToDomainDecomposition();
492 437 : if( dd ) dd->clearFullList();
493 : }
494 : break;
495 : /* ADDED WITH API==6 */
496 115 : case cmd_getDataRank:
497 : {
498 115 : CHECK_INIT(initialized,words[0]); plumed_assert(nw==2 || nw==3);
499 115 : std::string vtype=""; if( nw==3 ) vtype=" TYPE="+std::string(words[2]);
500 345 : readInputLine( "grab_" + std::string(words[1]) + ": GET ARG=" + std::string(words[1]) + vtype );
501 230 : ActionToGetData* as=actionSet.selectWithLabel<ActionToGetData*>("grab_"+std::string(words[1]));
502 115 : plumed_assert( as ); as->get_rank( val );
503 : }
504 115 : break;
505 : /* ADDED WITH API==6 */
506 51 : case cmd_getDataShape:
507 : {
508 51 : CHECK_INIT(initialized,std::string(words[0]));
509 102 : ActionToGetData* as1=actionSet.selectWithLabel<ActionToGetData*>("grab_"+std::string(words[1]));
510 51 : plumed_assert( as1 ); as1->get_shape( val );
511 : }
512 : break;
513 : /* ADDED WITH API==6 */
514 115 : case cmd_setMemoryForData:
515 : {
516 115 : CHECK_INIT(initialized,words[0]); plumed_assert(nw==2 || nw==3);
517 230 : ActionToGetData* as2=actionSet.selectWithLabel<ActionToGetData*>("grab_"+std::string(words[1]));
518 115 : plumed_assert( as2 ); as2->set_memory( val );
519 : }
520 : break;
521 : /* ADDED WITH API==6 */
522 : case cmd_setErrorHandler:
523 : {
524 0 : if(val) error_handler=*static_cast<const plumed_error_handler*>(val.get<const void*>());
525 0 : else error_handler.handler=NULL;
526 : }
527 : break;
528 0 : case cmd_read:
529 0 : CHECK_INIT(initialized,word);
530 0 : if(val)readInputFile(val.getCString());
531 0 : else readInputFile("plumed.dat");
532 : break;
533 819 : case cmd_readInputLine:
534 819 : CHECK_INIT(initialized,word);
535 819 : CHECK_NOTNULL(val,word);
536 819 : readInputLine(val.getCString());
537 777 : break;
538 2 : case cmd_readInputLines:
539 2 : CHECK_INIT(initialized,word);
540 2 : CHECK_NOTNULL(val,word);
541 2 : readInputLines(val.getCString());
542 2 : break;
543 3 : case cmd_clear:
544 : {
545 3 : CHECK_INIT(initialized,word);
546 : std::vector<int> natoms;
547 7 : for(const auto & pp : inputs ) {
548 4 : plumed_assert(pp);
549 4 : DomainDecomposition* dd=pp->castToDomainDecomposition();
550 4 : if ( dd ) natoms.push_back( dd->getNumberOfAtoms() );
551 : }
552 3 : actionSet.clearDelete(); inputs.clear();
553 4 : for(unsigned i=0; i<natoms.size(); ++i) {
554 1 : std::string str_natoms; Tools::convert( natoms[i], str_natoms );
555 3 : readInputLine( MDEngine + ": DOMAIN_DECOMPOSITION NATOMS=" + str_natoms +
556 2 : " VALUE1=posx UNIT1=length PERIODIC1=NO CONSTANT1=False ROLE1=x" +
557 2 : " VALUE2=posy UNIT2=length PERIODIC2=NO CONSTANT2=False ROLE2=y" +
558 2 : " VALUE3=posz UNIT3=length PERIODIC3=NO CONSTANT3=False ROLE3=z" +
559 2 : " VALUE4=Masses UNIT4=mass PERIODIC4=NO CONSTANT4=True ROLE4=m" +
560 : " VALUE5=Charges UNIT5=charge PERIODIC5=NO CONSTANT5=True ROLE5=q");
561 :
562 : }
563 3 : setUnits( passtools->usingNaturalUnits, passtools->units );
564 : }
565 3 : break;
566 : case cmd_getApiVersion:
567 44 : CHECK_NOTNULL(val,word);
568 44 : val.set(int(10));
569 : break;
570 : // commands which can be used only before initialization:
571 1289 : case cmd_init:
572 1289 : CHECK_NOTINIT(initialized,word);
573 1289 : init();
574 : break;
575 1044 : case cmd_setRealPrecision:
576 1044 : CHECK_NOTINIT(initialized,word);
577 1044 : CHECK_NOTNULL(val,word);
578 2087 : passtools=DataPassingTools::create(val.get<int>());
579 1043 : passtools->usingNaturalUnits=false;
580 1043 : break;
581 956 : case cmd_setMDLengthUnits:
582 956 : CHECK_NOTINIT(initialized,word);
583 956 : CHECK_NOTNULL(val,word);
584 956 : passtools->MDUnits.setLength(passtools->MD2double(val));
585 : break;
586 956 : case cmd_setMDChargeUnits:
587 956 : CHECK_NOTINIT(initialized,word);
588 956 : CHECK_NOTNULL(val,word);
589 956 : passtools->MDUnits.setCharge(passtools->MD2double(val));
590 : break;
591 956 : case cmd_setMDMassUnits:
592 956 : CHECK_NOTINIT(initialized,word);
593 956 : CHECK_NOTNULL(val,word);
594 956 : passtools->MDUnits.setMass(passtools->MD2double(val));
595 : break;
596 45 : case cmd_setMDEnergyUnits:
597 45 : CHECK_NOTINIT(initialized,word);
598 45 : CHECK_NOTNULL(val,word);
599 45 : passtools->MDUnits.setEnergy(passtools->MD2double(val));
600 : break;
601 6 : case cmd_setMDTimeUnits:
602 6 : CHECK_NOTINIT(initialized,word);
603 6 : CHECK_NOTNULL(val,word);
604 6 : passtools->MDUnits.setTime(passtools->MD2double(val));
605 : break;
606 0 : case cmd_setNaturalUnits:
607 : // set the boltzman constant for MD in natural units (kb=1)
608 : // only needed in LJ codes if the MD is passing temperatures to plumed (so, not yet...)
609 : // use as cmd("setNaturalUnits")
610 0 : CHECK_NOTINIT(initialized,word);
611 0 : passtools->usingNaturalUnits=true;
612 0 : break;
613 62 : case cmd_setNoVirial:
614 : {
615 62 : CHECK_NOTINIT(initialized,word);
616 62 : ActionToPutData* ap=actionSet.selectWithLabel<ActionToPutData*>("Box");
617 62 : if( ap ) ap->noforce=true;
618 : else {
619 10 : ActionForInterface* af = actionSet.selectWithLabel<ActionForInterface*>(MDEngine);
620 10 : if( af ) plumed_merror("setNoVirial should be called after number of atoms have been set");
621 : }
622 : }
623 : break;
624 1006 : case cmd_setPlumedDat:
625 1006 : CHECK_NOTINIT(initialized,word);
626 1006 : CHECK_NOTNULL(val,word);
627 1006 : plumedDat=val.getCString();
628 : break;
629 353 : case cmd_setMPIComm:
630 353 : CHECK_NOTINIT(initialized,word);
631 353 : comm.Set_comm(val);
632 355 : for(const auto & pp : inputs ) pp->Set_comm(comm);
633 : break;
634 0 : case cmd_setMPIFComm:
635 0 : CHECK_NOTINIT(initialized,word);
636 0 : comm.Set_fcomm(val);
637 0 : for(const auto & pp : inputs ) pp->Set_comm(comm);
638 : break;
639 0 : case cmd_setMPImultiSimComm:
640 0 : CHECK_NOTINIT(initialized,word);
641 0 : multi_sim_comm.Set_comm(val);
642 : break;
643 1273 : case cmd_setNatoms:
644 : {
645 1273 : CHECK_NOTINIT(initialized,word);
646 1273 : CHECK_NOTNULL(val,word);
647 1273 : int natoms = val.get<int>(); std::string str_natoms; Tools::convert( natoms, str_natoms );
648 1263 : ActionForInterface* dd=actionSet.selectWithLabel<ActionForInterface*>(MDEngine);
649 4845 : if( !dd && natoms>0 ) readInputLine( MDEngine + ": DOMAIN_DECOMPOSITION NATOMS=" + str_natoms + +
650 2388 : " VALUE1=posx UNIT1=length PERIODIC1=NO CONSTANT1=False ROLE1=x" +
651 2388 : " VALUE2=posy UNIT2=length PERIODIC2=NO CONSTANT2=False ROLE2=y" +
652 2388 : " VALUE3=posz UNIT3=length PERIODIC3=NO CONSTANT3=False ROLE3=z" +
653 2388 : " VALUE4=Masses UNIT4=mass PERIODIC4=NO CONSTANT4=True ROLE4=m" +
654 2388 : " VALUE5=Charges UNIT5=charge PERIODIC5=NO CONSTANT5=True ROLE5=q", true );
655 : }
656 1263 : break;
657 1024 : case cmd_setTimestep:
658 : {
659 1024 : CHECK_NOTINIT(initialized,word);
660 1024 : CHECK_NOTNULL(val,word);
661 1024 : ActionToPutData* ts = actionSet.selectWithLabel<ActionToPutData*>("timestep");
662 1024 : if( !ts ) {
663 1022 : readInputLine("timestep: PUT UNIT=time PERIODIC=NO CONSTANT", true);
664 2044 : ts = actionSet.selectWithLabel<ActionToPutData*>("timestep");
665 : }
666 2048 : if( !ts->setValuePointer("timestep", val ) ) plumed_error();
667 : // The following is to avoid extra digits in case the MD code uses floats
668 : // e.g.: float f=0.002 when converted to double becomes 0.002000000094995
669 : // To avoid this, we keep only up to 6 significant digits after first one
670 1024 : if( getRealPrecision()<=4 ) {
671 166 : Value* tstepv = ts->copyOutput(0);
672 0 : double magnitude=std::pow(10,std::floor(std::log10(tstepv->get())));
673 0 : tstepv->set( std::round(tstepv->get()/magnitude*1e6)/1e6*magnitude );
674 : }
675 1024 : ts->updateUnits( passtools.get() );
676 : }
677 : break;
678 : /* ADDED WITH API==2 */
679 61 : case cmd_setKbT:
680 : {
681 61 : CHECK_NOTINIT(initialized,word);
682 61 : CHECK_NOTNULL(val,word);
683 61 : readInputLine("kBT: PUT CONSTANT PERIODIC=NO UNIT=energy", true);
684 61 : ActionToPutData* kb = actionSet.selectWithLabel<ActionToPutData*>("kBT");
685 122 : if( !kb->setValuePointer("kBT", val ) ) plumed_error();
686 61 : kb->updateUnits( passtools.get() );
687 : }
688 : break;
689 : /* ADDED WITH API==3 */
690 8 : case cmd_setRestart:
691 8 : CHECK_NOTINIT(initialized,word);
692 8 : CHECK_NOTNULL(val,word);
693 8 : if(val.get<int>()!=0) restart=true;
694 : break;
695 : /* ADDED WITH API==4 */
696 0 : case cmd_doCheckPoint:
697 0 : CHECK_INIT(initialized,word);
698 0 : CHECK_NOTNULL(val,word);
699 0 : doCheckPoint = false;
700 0 : if(val.get<int>()!=0) doCheckPoint = true;
701 : break;
702 : /* ADDED WITH API==6 */
703 : case cmd_setNumOMPthreads:
704 0 : CHECK_NOTNULL(val,word);
705 : {
706 0 : auto nt=val.get<unsigned>();
707 : if(nt==0) nt=1;
708 0 : OpenMP::setNumThreads(nt);
709 : }
710 : break;
711 : /* ADDED WITH API==10 */
712 : case cmd_setGpuDeviceId:
713 0 : CHECK_NOTNULL(val,word);
714 : {
715 0 : auto id=val.get<int>();
716 0 : if(id>=0) gpuDeviceId=id;
717 : }
718 : break;
719 : /* ADDED WITH API==6 */
720 : /* only used for testing */
721 : case cmd_throw:
722 61 : CHECK_NOTNULL(val,word);
723 61 : testThrow(val.getCString());
724 : /* ADDED WITH API==10 */
725 : case cmd_setNestedExceptions:
726 27 : CHECK_NOTNULL(val,word);
727 27 : if(val.get<int>()!=0) nestedExceptions=true;
728 1 : else nestedExceptions=false;
729 : break;
730 : /* STOP API */
731 1021 : case cmd_setMDEngine:
732 1021 : CHECK_NOTINIT(initialized,word);
733 1021 : CHECK_NOTNULL(val,word);
734 1021 : MDEngine=val.getCString();
735 : break;
736 996 : case cmd_setLog:
737 996 : CHECK_NOTINIT(initialized,word);
738 996 : log.link(val.get<FILE*>());
739 : break;
740 169 : case cmd_setLogFile:
741 169 : CHECK_NOTINIT(initialized,word);
742 169 : CHECK_NOTNULL(val,word);
743 169 : log.open(val.getCString());
744 169 : break;
745 : // other commands that should be used after initialization:
746 266094 : case cmd_setStopFlag:
747 266094 : CHECK_INIT(initialized,word);
748 266094 : CHECK_NOTNULL(val,word);
749 266094 : val.get<int*>(); // just check type and discard pointer
750 266093 : stopFlag=val.copy();
751 266093 : break;
752 0 : case cmd_getExchangesFlag:
753 0 : CHECK_INIT(initialized,word);
754 0 : CHECK_NOTNULL(val,word);
755 0 : exchangePatterns.getFlag(*val.get<int*>()); // note: getFlag changes the value of the reference!
756 : break;
757 0 : case cmd_setExchangesSeed:
758 0 : CHECK_INIT(initialized,word);
759 0 : CHECK_NOTNULL(val,word);
760 0 : exchangePatterns.setSeed(val.get<int>());
761 : break;
762 0 : case cmd_setNumberOfReplicas:
763 0 : CHECK_INIT(initialized,word);
764 0 : CHECK_NOTNULL(val,word);
765 0 : exchangePatterns.setNofR(val.get<int>());
766 : break;
767 0 : case cmd_getExchangesList:
768 0 : CHECK_INIT(initialized,word);
769 0 : CHECK_NOTNULL(val,word);
770 0 : exchangePatterns.getList(val);
771 : break;
772 946 : case cmd_runFinalJobs:
773 946 : CHECK_INIT(initialized,word);
774 946 : runJobsAtEndOfCalculation();
775 : break;
776 96 : case cmd_isEnergyNeeded:
777 : {
778 96 : CHECK_INIT(initialized,word);
779 96 : CHECK_NOTNULL(val,word);
780 96 : if( name_of_energy =="" ) {
781 96 : val.set(int(0));
782 : } else {
783 0 : ActionToPutData* ap=actionSet.selectWithLabel<ActionToPutData*>(name_of_energy);
784 0 : if(ap->isActive()) val.set(int(1));
785 0 : else val.set(int(0));
786 : }
787 : }
788 : break;
789 7051 : case cmd_getBias:
790 7051 : CHECK_INIT(initialized,word);
791 7051 : CHECK_NOTNULL(val,word);
792 7051 : plumedQuantityToMD( "energy", getBias(), val );
793 7051 : break;
794 : case cmd_checkAction:
795 2 : CHECK_NOTNULL(val,word);
796 2 : plumed_assert(nw==2);
797 3 : val.set(int(actionRegister().check(dlloader.getHandles(), std::string(words[1])) ? 1:0));
798 2 : break;
799 : case cmd_setExtraCV:
800 : {
801 30 : CHECK_NOTNULL(val,word);
802 30 : plumed_assert(nw==2);
803 90 : if( valueExists(std::string(words[1])) ) setInputValue( std::string(words[1]), 0, 1, val );
804 : }
805 : break;
806 : case cmd_setExtraCVForce:
807 : {
808 30 : CHECK_NOTNULL(val,word); plumed_assert(nw==2);
809 120 : if( valueExists(std::string(words[1])) ) setInputForce( std::string(words[1]), val );
810 : }
811 : break;
812 : /* ADDED WITH API==10 */
813 : case cmd_isExtraCVNeeded:
814 10 : CHECK_NOTNULL(val,word);
815 10 : plumed_assert(nw==2); val.set(int(0));
816 56 : for(const auto & p : inputs) {
817 50 : if( p->getLabel()==words[1] && p->isActive() ) { val.set(int(1)); break; }
818 : }
819 : break;
820 1077 : case cmd_GREX:
821 1281 : if(!grex) grex=Tools::make_unique<GREX>(*this);
822 1077 : plumed_massert(grex,"error allocating grex");
823 : {
824 1077 : std::string kk=std::string(words[1]);
825 1419 : for(unsigned i=2; i<words.size(); i++) kk+=" "+std::string(words[i]);
826 1077 : grex->cmd(kk.c_str(),val);
827 : }
828 1077 : break;
829 16374 : case cmd_CLTool:
830 16374 : CHECK_NOTINIT(initialized,word);
831 21737 : if(!cltool) cltool=Tools::make_unique<CLToolMain>();
832 : {
833 16374 : std::string kk(words[1]);
834 16374 : for(unsigned i=2; i<words.size(); i++) kk+=" "+std::string(words[i]);
835 16374 : cltool->cmd(kk.c_str(),val);
836 : }
837 16374 : break;
838 : break;
839 : /* ADDED WITH API==7 */
840 : case cmd_convert:
841 : {
842 : double v;
843 57 : plumed_assert(words.size()==2);
844 280 : if(Tools::convertNoexcept(std::string(words[1]),v)) passtools->double2MD(v,val);
845 : }
846 57 : break;
847 7 : default:
848 14 : plumed_error() << "cannot interpret cmd(\"" << word << "\"). check plumed developers manual to see the available commands.";
849 : break;
850 : }
851 : }
852 :
853 1331273 : } catch (...) {
854 166 : if(log.isOpen()) {
855 : try {
856 92 : log<<"\n################################################################################\n";
857 92 : log<<Tools::concatenateExceptionMessages();
858 92 : log<<"\n################################################################################\n";
859 92 : log.flush();
860 0 : } catch(...) {
861 : // ignore errors here.
862 : // in any case, we are rethrowing this below
863 0 : }
864 : }
865 166 : throw;
866 166 : }
867 1330941 : }
868 :
869 : ////////////////////////////////////////////////////////////////////////
870 :
871 1289 : void PlumedMain::init() {
872 : // check that initialization just happens once
873 1289 : initialized=true;
874 1289 : if(!log.isOpen()) log.link(stdout);
875 1289 : log<<"PLUMED is starting\n";
876 3867 : log<<"Version: "<<config::getVersionLong()<<" (git: "<<config::getVersionGit()<<") "
877 5156 : <<"compiled on " <<config::getCompilationDate() << " at " << config::getCompilationTime() << "\n";
878 1289 : log<<"Please cite these papers when using PLUMED ";
879 2578 : log<<cite("The PLUMED consortium, Nat. Methods 16, 670 (2019)");
880 2578 : log<<cite("Tribello, Bonomi, Branduardi, Camilloni, and Bussi, Comput. Phys. Commun. 185, 604 (2014)");
881 1289 : log<<"\n";
882 1289 : log<<"For further information see the PLUMED web page at http://www.plumed.org\n";
883 1289 : log<<"Root: "<<config::getPlumedRoot()<<"\n";
884 1289 : log<<"LibraryPath: "<<config::getLibraryPath()<<"\n";
885 2578 : log<<"For installed feature, see "<<config::getPlumedRoot() + "/src/config/config.txt\n";
886 1289 : log.printf("Molecular dynamics engine: %s\n",MDEngine.c_str());
887 1289 : log.printf("Precision of reals: %d\n",passtools->getRealPrecision());
888 2338 : log.printf("Running over %d %s\n",comm.Get_size(),(comm.Get_size()>1?"nodes":"node"));
889 1289 : log<<"Number of threads: "<<OpenMP::getNumThreads()<<"\n";
890 1289 : log<<"Cache line size: "<<OpenMP::getCachelineSize()<<"\n";
891 4778 : for(const auto & pp : inputs ) {
892 3489 : plumed_assert(pp);
893 3489 : DomainDecomposition* dd=pp->castToDomainDecomposition();
894 3489 : if ( dd ) log.printf("Number of atoms: %d\n",dd->getNumberOfAtoms());
895 : }
896 1289 : if(grex) log.printf("GROMACS-like replica exchange is on\n");
897 1289 : log.printf("File suffix: %s\n",getSuffix().c_str());
898 1289 : if(plumedDat.length()>0) {
899 1006 : readInputFile(plumedDat);
900 : plumedDat="";
901 : }
902 1289 : setUnits( passtools->usingNaturalUnits, passtools->units );
903 1289 : ActionToPutData* ts = actionSet.selectWithLabel<ActionToPutData*>("timestep");
904 1289 : if(ts) log.printf("Timestep: %f\n",(ts->copyOutput(0))->get());
905 1289 : ActionToPutData* kb = actionSet.selectWithLabel<ActionToPutData*>("kBT");
906 1289 : if(kb)
907 61 : log.printf("KbT: %f\n",(kb->copyOutput(0))->get());
908 : else {
909 1228 : log.printf("KbT has not been set by the MD engine\n");
910 1228 : log.printf("It should be set by hand where needed\n");
911 : }
912 1289 : log<<"Relevant bibliography:\n";
913 1289 : log<<citations;
914 1289 : log<<"Please read and cite where appropriate!\n";
915 1289 : log<<"Finished setup\n";
916 1289 : }
917 :
918 52185 : void PlumedMain::setupInterfaceActions() {
919 103099 : inputs.clear(); std::vector<ActionForInterface*> ap=actionSet.select<ActionForInterface*>();
920 423003 : for(unsigned i=0; i<ap.size(); ++i) {
921 370818 : if( ap[i]->getName()=="ENERGY" || ap[i]->getDependencies().size()==0 ) inputs.push_back( ap[i] );
922 : }
923 52185 : }
924 :
925 1029 : void PlumedMain::readInputFile(const std::string & str) {
926 1029 : plumed_assert(initialized);
927 1029 : log<<"FILE: "<<str<<"\n";
928 1029 : IFile ifile;
929 1029 : ifile.link(*this);
930 1029 : ifile.open(str);
931 1029 : ifile.allowNoEOL();
932 1029 : readInputFile(ifile);
933 1029 : log<<"END FILE: "<<str<<"\n";
934 1029 : log.flush();
935 :
936 1029 : }
937 :
938 1031 : void PlumedMain::readInputFile(IFile & ifile) {
939 : std::vector<std::string> words;
940 19501 : while(Tools::getParsedLine(ifile,words) && !endPlumed) readInputWords(words,false);
941 1031 : endPlumed=false;
942 1031 : pilots=actionSet.select<ActionPilot*>();
943 1031 : setupInterfaceActions();
944 1031 : }
945 :
946 32861 : void PlumedMain::readInputLine(const std::string & str, const bool& before_init) {
947 32861 : if( !before_init ) plumed_assert(initialized);
948 32861 : if(str.empty()) return;
949 32861 : std::vector<std::string> words=Tools::getWords(str);
950 32861 : citations.clear();
951 32861 : readInputWords(words,before_init);
952 32806 : if(!citations.empty()) {
953 487 : log<<"Relevant bibliography:\n";
954 487 : log<<citations;
955 487 : log<<"Please read and cite where appropriate!\n";
956 : }
957 32861 : }
958 :
959 2 : void PlumedMain::readInputLines(const std::string & str) {
960 2 : plumed_assert(initialized);
961 2 : if(str.empty()) return;
962 :
963 2 : log<<"FILE: (temporary)\n";
964 :
965 : // Open a temporary file
966 2 : auto fp=std::tmpfile();
967 2 : plumed_assert(fp);
968 :
969 : // make sure file is closed (and thus deleted) also if an exception occurs
970 2 : auto deleter=[](auto fp) { std::fclose(fp); };
971 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
972 :
973 2 : auto ret=std::fputs(str.c_str(),fp);
974 2 : plumed_assert(ret!=EOF);
975 :
976 2 : std::rewind(fp);
977 :
978 2 : IFile ifile;
979 2 : ifile.link(*this);
980 2 : ifile.link(fp);
981 2 : ifile.allowNoEOL();
982 :
983 2 : readInputFile(ifile);
984 2 : log<<"END FILE: (temporary)\n";
985 2 : }
986 :
987 51381 : void PlumedMain::readInputWords(const std::vector<std::string> & words, const bool& before_init) {
988 51381 : if( !before_init ) plumed_assert(initialized);
989 51381 : if(words.empty())return;
990 51209 : else if(words[0]=="_SET_SUFFIX") {
991 3 : plumed_assert(words.size()==2);
992 : setSuffix(words[1]);
993 : } else {
994 51206 : std::vector<std::string> interpreted(words);
995 51206 : Tools::interpretLabel(interpreted);
996 102357 : auto action=actionRegister().create(dlloader.getHandles(),ActionOptions(*this,interpreted));
997 51151 : if(!action) {
998 : std::string msg;
999 : msg ="ERROR\nI cannot understand line:";
1000 0 : for(unsigned i=0; i<interpreted.size(); ++i) msg+=" "+interpreted[i];
1001 : msg+="\nMaybe a missing space or a typo?";
1002 0 : log << msg;
1003 0 : log.flush();
1004 0 : plumed_merror(msg);
1005 : }
1006 51151 : action->checkRead();
1007 51151 : actionSet.emplace_back(std::move(action));
1008 51206 : };
1009 :
1010 51154 : pilots=actionSet.select<ActionPilot*>();
1011 51154 : setupInterfaceActions();
1012 : }
1013 :
1014 : ////////////////////////////////////////////////////////////////////////
1015 :
1016 0 : void PlumedMain::exit(int c) {
1017 0 : comm.Abort(c);
1018 0 : }
1019 :
1020 51204 : Log& PlumedMain::getLog() {
1021 51204 : return log;
1022 : }
1023 :
1024 278008 : void PlumedMain::calc() {
1025 278008 : prepareCalc();
1026 277992 : performCalc();
1027 277990 : }
1028 :
1029 284959 : void PlumedMain::prepareCalc() {
1030 284959 : prepareDependencies();
1031 284959 : shareData();
1032 284943 : }
1033 :
1034 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1035 : // here we have the main steps in "calc()"
1036 : // they can be called individually, but the standard thing is to
1037 : // traverse them in this order:
1038 285286 : void PlumedMain::prepareDependencies() {
1039 :
1040 : // Stopwatch is stopped when sw goes out of scope
1041 285286 : auto sw=stopwatch.startStop("1 Prepare dependencies");
1042 :
1043 : // activate all the actions which are on step
1044 : // activation is recursive and enables also the dependencies
1045 : // before doing that, the prepare() method is called to see if there is some
1046 : // new/changed dependency (up to now, only useful for dependences on virtual atoms,
1047 : // which can be dynamically changed).
1048 :
1049 : // First switch off all actions
1050 4648601 : for(const auto & p : actionSet) {
1051 4363315 : p->deactivate();
1052 : }
1053 :
1054 : // for optimization, an "active" flag remains false if no action at all is active
1055 285286 : active=false;
1056 1720265 : for(unsigned i=0; i<pilots.size(); ++i) {
1057 1434979 : if(pilots[i]->onStep()) {
1058 1341728 : pilots[i]->activate();
1059 1341728 : active=true;
1060 : }
1061 : };
1062 :
1063 : // This stops the driver calculation if there is not a read action
1064 285286 : if( !active && !inputsAreActive() ) stopFlag.set(int(1));
1065 :
1066 : // also, if one of them is the total energy, tell to atoms that energy should be collected
1067 4648601 : for(const auto & p : actionSet) {
1068 4363315 : if(p->isActive()) {
1069 2622876 : if(p->checkNeedsGradients()) p->setOption("GRADIENTS");
1070 : }
1071 : }
1072 :
1073 285286 : }
1074 :
1075 1763 : bool PlumedMain::inputsAreActive() const {
1076 5178 : for(const auto & ip : inputs) {
1077 5122 : if( ip->onStep() ) return true;
1078 : }
1079 : return false;
1080 : }
1081 :
1082 114 : void PlumedMain::shareAll() {
1083 456 : for(const auto & ip : inputs) ip->shareAll();
1084 114 : }
1085 :
1086 285043 : void PlumedMain::shareData() {
1087 : // atom positions are shared (but only if there is something to do)
1088 285043 : if(!active)return;
1089 : // Stopwatch is stopped when sw goes out of scope
1090 283340 : auto sw=stopwatch.startStop("2 Sharing data");
1091 743782 : for(const auto & ip : inputs) ip->share();
1092 283340 : }
1093 :
1094 7011 : void PlumedMain::performCalcNoUpdate() {
1095 7011 : waitData();
1096 7011 : justCalculate();
1097 7011 : backwardPropagate();
1098 7011 : resetInputs();
1099 7011 : }
1100 :
1101 10 : void PlumedMain::performCalcNoForces() {
1102 10 : waitData();
1103 10 : justCalculate();
1104 10 : }
1105 :
1106 278002 : void PlumedMain::performCalc() {
1107 278002 : waitData();
1108 278002 : justCalculate();
1109 278002 : backwardPropagate();
1110 278000 : update();
1111 278000 : resetInputs();
1112 278000 : }
1113 :
1114 285137 : void PlumedMain::waitData() {
1115 285137 : if(!active)return;
1116 : // Stopwatch is stopped when sw goes out of scope
1117 283434 : auto sw=stopwatch.startStop("3 Waiting for data");
1118 744194 : for(const auto & ip : inputs) {
1119 460760 : if( ip->isActive() && ip->hasBeenSet() ) ip->wait();
1120 300385 : else if( ip->isActive() ) ip->warning("input requested but this quantity has not been set");
1121 : }
1122 283434 : }
1123 :
1124 285137 : void PlumedMain::justCalculate() {
1125 285137 : if(!active)return;
1126 : // Stopwatch is stopped when sw goes out of scope
1127 283434 : auto sw=stopwatch.startStop("4 Calculating (forward loop)");
1128 283434 : bias=0.0;
1129 283434 : work=0.0;
1130 :
1131 : // Check the input actions to determine if we need to calculate constants that
1132 : // depend on masses and charges
1133 : bool firststep=false;
1134 744194 : for(const auto & ip : inputs) {
1135 460760 : if( ip->firststep ) firststep=true;
1136 : }
1137 286741 : if( firststep ) { for(const auto & ip : inputs) ip->firststep=false; }
1138 :
1139 : int iaction=0;
1140 : // calculate the active actions in order (assuming *backward* dependence)
1141 4615751 : for(const auto & pp : actionSet) {
1142 : Action* p(pp.get());
1143 4332317 : plumed_assert(p);
1144 : try {
1145 4332317 : if(p->isActive()) {
1146 : // Stopwatch is stopped when sw goes out of scope.
1147 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
1148 2620795 : Stopwatch::Handler sw;
1149 2620795 : if(detailedTimers) {
1150 2481 : auto actionNumberLabel=std::to_string(iaction);
1151 2481 : const unsigned m=actionSet.size();
1152 7443 : unsigned k=0; unsigned n=1; while(n<m) { n*=10; k++; }
1153 2481 : auto spaces=std::string(k-actionNumberLabel.length(),' ');
1154 2481 : sw=stopwatch.startStop("4A " + spaces + actionNumberLabel+" "+p->getLabel());
1155 : }
1156 2620795 : ActionWithValue*av=p->castToActionWithValue();
1157 2620795 : ActionAtomistic*aa=p->castToActionAtomistic();
1158 : {
1159 2620795 : if(av) av->clearInputForces();
1160 2620795 : if(av) av->clearDerivatives();
1161 2620795 : if( av && av->calculateOnUpdate() ) continue ;
1162 : }
1163 : {
1164 2616780 : if(aa) if(aa->isActive()) aa->retrieveAtoms();
1165 : }
1166 2616780 : if(p->checkNumericalDerivatives()) p->calculateNumericalDerivatives();
1167 2616530 : else p->calculate();
1168 : // This retrieves components called bias
1169 2616780 : if(av) {
1170 2339469 : bias+=av->getOutputQuantity("bias");
1171 2339469 : work+=av->getOutputQuantity("work");
1172 2339469 : av->setGradientsIfNeeded();
1173 : }
1174 : // This makes all values that depend on the (fixed) masses and charges constant
1175 2616780 : if( firststep ) p->setupConstantValues( true );
1176 2616780 : ActionWithVirtualAtom*avv=p->castToActionWithVirtualAtom();
1177 2616780 : if(avv)avv->setGradientsIfNeeded();
1178 2620795 : }
1179 0 : } catch(...) {
1180 0 : plumed_error_nested() << "An error happened while calculating " << p->getLabel();
1181 0 : }
1182 4328302 : iaction++;
1183 : }
1184 283434 : }
1185 :
1186 0 : void PlumedMain::justApply() {
1187 0 : backwardPropagate();
1188 0 : update();
1189 0 : resetInputs();
1190 0 : }
1191 :
1192 285013 : void PlumedMain::backwardPropagate() {
1193 285013 : if(!active)return;
1194 : int iaction=0;
1195 : // Stopwatch is stopped when sw goes out of scope
1196 283310 : auto sw=stopwatch.startStop("5 Applying (backward loop)");
1197 : // apply them in reverse order
1198 4613285 : for(auto pp=actionSet.rbegin(); pp!=actionSet.rend(); ++pp) {
1199 : const auto & p(pp->get());
1200 4329977 : if(p->isActive()) {
1201 :
1202 : // Stopwatch is stopped when sw goes out of scope.
1203 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
1204 2619119 : Stopwatch::Handler sw;
1205 2619119 : if(detailedTimers) {
1206 2481 : auto actionNumberLabel=std::to_string(iaction);
1207 2481 : const unsigned m=actionSet.size();
1208 7443 : unsigned k=0; unsigned n=1; while(n<m) { n*=10; k++; }
1209 2481 : auto spaces=std::string(k-actionNumberLabel.length(),' ');
1210 2481 : sw=stopwatch.startStop("5A " + spaces + actionNumberLabel+" "+p->getLabel());
1211 : }
1212 :
1213 2619119 : p->apply();
1214 2619119 : }
1215 4329975 : iaction++;
1216 : }
1217 :
1218 : // Stopwatch is stopped when sw goes out of scope.
1219 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
1220 283308 : Stopwatch::Handler sw1;
1221 283308 : if(detailedTimers) sw1=stopwatch.startStop("5B Update forces");
1222 283310 : }
1223 :
1224 278079 : void PlumedMain::update() {
1225 278079 : if(!active)return;
1226 :
1227 : // Stopwatch is stopped when sw goes out of scope
1228 276376 : auto sw=stopwatch.startStop("6 Update");
1229 :
1230 : // update step (for statistics, etc)
1231 276376 : updateFlags.push(true);
1232 4389066 : for(const auto & p : actionSet) {
1233 4112690 : p->beforeUpdate();
1234 6636800 : if(p->isActive() && p->checkUpdate() && updateFlagsTop()) {
1235 2524090 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(p.get());
1236 2524090 : if( av && av->calculateOnUpdate() ) { p->prepare(); p->calculate(); }
1237 2520075 : else p->update();
1238 : }
1239 : }
1240 552756 : while(!updateFlags.empty()) updateFlags.pop();
1241 : if(!updateFlags.empty()) plumed_merror("non matching changes in the update flags");
1242 : // Check that no action has told the calculation to stop
1243 276376 : if(stopNow) {
1244 65 : if(stopFlag) stopFlag.set(int(1));
1245 0 : else plumed_merror("your md code cannot handle plumed stop events - add a call to plumed.comm(stopFlag,stopCondition)");
1246 : }
1247 :
1248 : // flush by default every 10000 steps
1249 : // hopefully will not affect performance
1250 : // also if receive checkpointing signal
1251 276376 : if(step%10000==0||doCheckPoint) {
1252 1204 : fflush();
1253 1204 : log.flush();
1254 53332 : for(const auto & p : actionSet) p->fflush();
1255 : }
1256 276376 : }
1257 :
1258 42 : void PlumedMain::load(const std::string& fileName) {
1259 42 : if(DLLoader::installed()) {
1260 42 : std::string libName=fileName;
1261 42 : size_t n=libName.find_last_of(".");
1262 43 : std::string extension="";
1263 42 : std::string base=libName;
1264 42 : if(n!=std::string::npos && n<libName.length()-1)
1265 84 : extension=libName.substr(n+1);
1266 42 : if(n!=std::string::npos && n<libName.length())
1267 84 : base=libName.substr(0,n);
1268 :
1269 42 : if(extension=="cpp") {
1270 46 : libName="./"+base+"."+config::getVersionLong()+"."+config::getSoExt();
1271 : // full path command, including environment setup
1272 : // this will work even if plumed is not in the execution path or if it has been
1273 : // installed with a name different from "plumed"
1274 46 : std::string cmd=config::getEnvCommand()+" \""+config::getPlumedRoot()+"\"/scripts/mklib.sh -n -o "+libName+" "+fileName;
1275 :
1276 23 : if(std::getenv("PLUMED_LOAD_ACTION_DEBUG")) log<<"Executing: "<<cmd;
1277 23 : else log<<"Compiling: "<<fileName<<" to "<<libName;
1278 :
1279 23 : if(comm.Get_size()>0) log<<" (only on master node)";
1280 23 : log<<"\n";
1281 :
1282 : // On MPI process (intracomm), we use Get_rank to make sure a single process does the compilation
1283 : // Processes from multiple replicas might simultaneously do the compilation.
1284 23 : if(comm.Get_rank()==0) {
1285 23 : static Tools::CriticalSectionWithKey<std::string> section;
1286 : // This is only locking commands that are running with identical arguments.
1287 : // It is not necessary for correctness (a second mklib would just result in a no op since
1288 : // the library is already there, even if running simultaneously).
1289 : // It however decreases the system load if many threads are used.
1290 : auto s=section.startStop(cmd);
1291 23 : int ret=std::system(cmd.c_str());
1292 24 : if(ret!=0) plumed_error() <<"An error happened while executing command "<<cmd<<"\n";
1293 23 : }
1294 22 : comm.Barrier();
1295 : } else {
1296 38 : libName=base+"."+config::getSoExt();
1297 : }
1298 :
1299 : // If we have multiple threads (each holding a Plumed object), each of them
1300 : // will load the library, but each of them will only see actions registered
1301 : // from the owned library
1302 41 : auto *p=dlloader.load(libName);
1303 42 : log<<"Loading shared library "<<libName.c_str()<<" at "<<p<<"\n";
1304 41 : log<<"Here is the list of new actions\n";
1305 41 : log<<"\n";
1306 82 : for(const auto & a : actionRegister().getKeysWithDLHandle(p)) log<<a<<"\n";
1307 41 : log<<"\n";
1308 : } else {
1309 0 : plumed_error()<<"While loading library "<< fileName << " loading was not enabled, please check if dlopen was found at configure time";
1310 : }
1311 41 : }
1312 :
1313 285011 : void PlumedMain::resetInputs() {
1314 750530 : for(const auto & ip : inputs) {
1315 465519 : if( ip->isActive() && ip->hasBeenSet() ) ip->reset();
1316 : }
1317 285011 : }
1318 :
1319 7684 : double PlumedMain::getBias() const {
1320 7684 : return bias;
1321 : }
1322 :
1323 450 : double PlumedMain::getWork() const {
1324 450 : return work;
1325 : }
1326 :
1327 82 : FILE* PlumedMain::fopen(const char *path, const char *mode) {
1328 82 : std::string mmode(mode);
1329 83 : std::string ppath(path);
1330 82 : std::string suffix(getSuffix());
1331 82 : std::string ppathsuf=ppath+suffix;
1332 82 : FILE*fp=std::fopen(const_cast<char*>(ppathsuf.c_str()),const_cast<char*>(mmode.c_str()));
1333 82 : if(!fp) fp=std::fopen(const_cast<char*>(ppath.c_str()),const_cast<char*>(mmode.c_str()));
1334 84 : plumed_massert(fp,"file " + ppath + " cannot be found");
1335 81 : return fp;
1336 : }
1337 :
1338 99 : int PlumedMain::fclose(FILE*fp) {
1339 99 : return std::fclose(fp);
1340 : }
1341 :
1342 4694 : std::string PlumedMain::cite(const std::string&item) {
1343 4694 : return citations.cite(item);
1344 : }
1345 :
1346 1765 : void PlumedMain::fflush() {
1347 5282 : for(const auto & p : files) {
1348 3517 : p->flush();
1349 : }
1350 1765 : }
1351 :
1352 5086 : void PlumedMain::insertFile(FileBase&f) {
1353 5086 : files.insert(&f);
1354 5086 : }
1355 :
1356 5352 : void PlumedMain::eraseFile(FileBase&f) {
1357 5352 : files.erase(&f);
1358 5352 : }
1359 :
1360 65 : void PlumedMain::stop() {
1361 65 : stopNow=true;
1362 65 : }
1363 :
1364 946 : void PlumedMain::runJobsAtEndOfCalculation() {
1365 47667 : for(const auto & p : actionSet) {
1366 46721 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(p.get());
1367 46721 : if( av && av->calculateOnUpdate() ) p->activate();
1368 : }
1369 47667 : for(const auto & p : actionSet) {
1370 46721 : ActionPilot* ap=dynamic_cast<ActionPilot*>(p.get());
1371 46721 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(p.get());
1372 46721 : if( av && av->calculateOnUpdate() ) { p->calculate(); }
1373 46307 : else if( ap && !av && ap->getStride()==0 ) { p->update(); }
1374 46237 : else p->runFinalJobs();
1375 : }
1376 946 : }
1377 :
1378 8807369 : unsigned PlumedMain::increaseReferenceCounter() noexcept {
1379 8807369 : return ++referenceCounter;
1380 : }
1381 :
1382 8806290 : unsigned PlumedMain::decreaseReferenceCounter() noexcept {
1383 8806290 : return --referenceCounter;
1384 : }
1385 :
1386 42 : unsigned PlumedMain::useCountReferenceCounter() const noexcept {
1387 42 : return referenceCounter;
1388 : }
1389 :
1390 60 : bool PlumedMain::valueExists( const std::string& name ) const {
1391 240 : for(const auto & p : inputs) {
1392 240 : if( p->getLabel()==name ) return true;
1393 : }
1394 : return false;
1395 : }
1396 :
1397 451100 : void PlumedMain::setInputValue( const std::string& name, const unsigned& start, const unsigned& stride, const TypesafePtr & val ) {
1398 : bool found=false;
1399 1221600 : for(const auto & pp : inputs) {
1400 1221600 : if( pp->setValuePointer( name, val ) ) { pp->setStart(name, start); pp->setStride(name, stride); found=true; break; }
1401 : }
1402 0 : plumed_massert( found, "found no action to set named " + name );
1403 451096 : }
1404 :
1405 282703 : void PlumedMain::setInputForce( const std::string& name, const TypesafePtr & val ) {
1406 : bool found=false;
1407 759059 : for(const auto & pp : inputs) {
1408 759059 : if( pp->setForcePointer( name, val ) ) { found=true; break; }
1409 : }
1410 282685 : plumed_massert( found, "found no action to set named " + name );
1411 282685 : }
1412 :
1413 1326 : void PlumedMain::setUnits( const bool& natural, const Units& u ) {
1414 1326 : passtools->usingNaturalUnits = natural; passtools->units=u;
1415 1326 : std::vector<ActionToPutData*> idata = actionSet.select<ActionToPutData*>();
1416 9873 : for(const auto & ip : idata) ip->updateUnits( passtools.get() );
1417 51012 : for(const auto & p : actionSet ) p->resetStoredTimestep();
1418 1326 : }
1419 :
1420 284956 : void PlumedMain::startStep() {
1421 750267 : for(const auto & ip : inputs) ip->resetForStepStart();
1422 284956 : }
1423 :
1424 114 : void PlumedMain::writeBinary(std::ostream&o)const {
1425 456 : for(const auto & ip : inputs) ip->writeBinary(o);
1426 114 : }
1427 :
1428 114 : void PlumedMain::readBinary(std::istream&i) {
1429 456 : for(const auto & ip : inputs) ip->readBinary(i);
1430 114 : }
1431 :
1432 40 : void PlumedMain::setEnergyValue( const std::string& name ) {
1433 40 : name_of_energy = name;
1434 40 : }
1435 :
1436 9457 : int PlumedMain::getRealPrecision() const {
1437 9457 : return passtools->getRealPrecision();
1438 : }
1439 :
1440 5026 : bool PlumedMain::usingNaturalUnits() const {
1441 5026 : return passtools->usingNaturalUnits;
1442 : }
1443 :
1444 14038994 : const Units& PlumedMain::getUnits() {
1445 14038994 : return passtools->units;
1446 : }
1447 :
1448 4 : PlumedMain::DeprecatedAtoms& PlumedMain::getAtoms() {
1449 4 : return datoms;
1450 : }
1451 :
1452 7222 : void PlumedMain::plumedQuantityToMD( const std::string& unit, const double& eng, const TypesafePtr & m) const {
1453 7222 : passtools->double2MD( eng/passtools->getUnitConversion(unit),m );
1454 7222 : }
1455 :
1456 0 : double PlumedMain::MDQuantityToPLUMED( const std::string& unit, const TypesafePtr & m) const {
1457 0 : double x=passtools->MD2double(m);
1458 0 : return x*passtools->getUnitConversion(unit);
1459 : }
1460 :
1461 1 : double PlumedMain::DeprecatedAtoms::getKBoltzmann() const {
1462 1 : if( plumed.usingNaturalUnits() ) return 1.0;
1463 1 : return kBoltzmann/plumed.getUnits().getEnergy();
1464 : }
1465 :
1466 1 : double PlumedMain::DeprecatedAtoms::getKbT() const {
1467 1 : ActionForInterface* kb=plumed.getActionSet().selectWithLabel<ActionForInterface*>("kBT");
1468 1 : if( kb ) return (kb->copyOutput(0))->get();
1469 : return 0.0;
1470 : }
1471 :
1472 1 : int PlumedMain::DeprecatedAtoms::getNatoms() const {
1473 1 : std::vector<ActionToPutData*> inputs=plumed.getActionSet().select<ActionToPutData*>();
1474 1 : for(const auto & pp : inputs ) {
1475 2 : if( pp->getRole()=="x" ) return (pp->copyOutput(0))->getShape()[0];
1476 : }
1477 : return 0;
1478 : }
1479 :
1480 1 : bool PlumedMain::DeprecatedAtoms::usingNaturalUnits() const {
1481 1 : return plumed.usingNaturalUnits();
1482 : }
1483 :
1484 0 : void PlumedMain::DeprecatedAtoms::setCollectEnergy(bool b) const {
1485 0 : plumed.readInputLine( plumed.MDEngine + "_energy: ENERGY" );
1486 0 : plumed.setEnergyValue( plumed.MDEngine + "_energy" );
1487 0 : }
1488 :
1489 0 : double PlumedMain::DeprecatedAtoms::getEnergy() const {
1490 0 : ActionToPutData* av = plumed.getActionSet().selectWithLabel<ActionToPutData*>( plumed.MDEngine + "_energy" );
1491 0 : return (av->copyOutput(0))->get();
1492 : }
1493 :
1494 5 : void PlumedMain::activateParseOnlyMode() {
1495 5 : doParseOnly=true;
1496 5 : }
1497 :
1498 2166 : bool PlumedMain::parseOnlyMode() const {
1499 2166 : return doParseOnly;
1500 : }
1501 :
1502 96 : void PlumedMain::getKeywordsForAction( const std::string& action, Keywords& keys ) const {
1503 96 : actionRegister().getKeywords( dlloader.getHandles(), action, keys );
1504 96 : }
1505 :
1506 : #ifdef __PLUMED_HAS_PYTHON
1507 : // This is here to stop cppcheck throwing an error
1508 : #endif
1509 :
1510 : #ifdef __PLUMED_HAS_DLADDR
1511 : // This is here to stop cppcheck throwing an error
1512 : #endif
1513 :
1514 : }
1515 :
1516 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
