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 <cmath>
51 : #include <cstdlib>
52 : #include <cstdio>
53 : #include <cstring>
54 : #include <set>
55 : #include <exception>
56 : #include <iomanip>
57 : #include <sstream>
58 : #include <stdexcept>
59 : #include <ios>
60 : #include <new>
61 : #include <typeinfo>
62 : #include <iostream>
63 : #include <algorithm>
64 : #include <system_error>
65 : #include <future>
66 : #include <memory>
67 : #include <functional>
68 : #include <regex>
69 : #include <any>
70 : #include <optional>
71 : #include <variant>
72 : #include <filesystem>
73 :
74 : namespace PLMD {
75 :
76 : namespace {
77 :
78 7371152 : const std::set<std::string>& nonFiniteActionTraceLabels() {
79 2608 : static const std::set<std::string> labels=[]() {
80 : std::set<std::string> result;
81 1304 : const char* environment=std::getenv("PLUMED_NONFINITE_ACTION_TRACE");
82 1304 : if(!environment) {
83 : return result;
84 : }
85 5 : std::string selection(environment);
86 : std::replace(selection.begin(),selection.end(),',',' ');
87 5 : std::istringstream words(selection);
88 : std::string label;
89 12 : while(words >> label) {
90 7 : result.insert(label);
91 : }
92 : return result;
93 5 : }
94 7371152 : ();
95 7371152 : return labels;
96 : }
97 :
98 2402 : bool traceActionSelected(const std::string& label) {
99 2402 : const std::set<std::string>& labels=nonFiniteActionTraceLabels();
100 2402 : return labels.count("*")>0 || labels.count(label)>0;
101 : }
102 :
103 187 : void traceActionNonFiniteLocal(const ActionSet& actions, const Action& trigger,
104 : const long long step, const char* phase,
105 : const bool checkValues, const bool checkForces) {
106 187 : if(nonFiniteActionTraceLabels().empty()) {
107 : return;
108 : }
109 2583 : for(const auto& actionPointer : actions) {
110 : Action* action=actionPointer.get();
111 : // Forward data are only current after their owning action has run.
112 : // Backward atom/box interface forces must also be checked: finite scalar
113 : // forces and derivatives can overflow when multiplied or accumulated.
114 2402 : const bool interfaceForce=checkForces &&
115 713 : dynamic_cast<ActionForInterface*>(action);
116 2402 : if(!traceActionSelected(action->getLabel()) && !interfaceForce) {
117 1672 : continue;
118 : }
119 730 : if(checkValues && action!=&trigger) {
120 192 : continue;
121 : }
122 538 : ActionWithValue* actionWithValue=action->castToActionWithValue();
123 538 : if(!actionWithValue) {
124 0 : continue;
125 : }
126 1037 : for(unsigned component=0;
127 1037 : component<actionWithValue->getNumberOfComponents(); ++component) {
128 505 : Value* value=actionWithValue->copyOutput(component);
129 505 : const std::size_t numberOfValues=value->getNumberOfValues();
130 505 : if(checkValues && value->valueHasBeenSet()) {
131 54 : for(std::size_t element=0; element<numberOfValues; ++element) {
132 27 : const double current=value->get(element);
133 27 : if(std::isfinite(current)) {
134 : continue;
135 : }
136 0 : std::ostringstream message;
137 : message << std::setprecision(17)
138 0 : << "PLUMED action non-finite trace: step=" << step
139 : << " phase=" << phase
140 : << " after_action=" << trigger.getLabel()
141 : << " after_action_type=" << trigger.getName()
142 : << " target=" << value->getName()
143 0 : << " field=value element=" << element
144 0 : << " value=" << current;
145 0 : plumed_merror(message.str());
146 0 : }
147 : // Scalar actions store a derivative array. Grid derivatives have a
148 : // different layout and are deliberately not treated as scalar data.
149 27 : if(value->hasDerivatives() && value->getRank()==0) {
150 : for(unsigned derivative=0;
151 144 : derivative<value->getNumberOfDerivatives(); ++derivative) {
152 : const double current=value->getDerivative(derivative);
153 124 : if(std::isfinite(current)) {
154 : continue;
155 : }
156 3 : std::ostringstream message;
157 : message << std::setprecision(17)
158 3 : << "PLUMED action non-finite trace: step=" << step
159 : << " phase=" << phase
160 : << " after_action=" << trigger.getLabel()
161 : << " after_action_type=" << trigger.getName()
162 : << " target=" << value->getName()
163 12 : << " field=derivative element=" << derivative
164 3 : << " value=" << current;
165 6 : plumed_merror(message.str());
166 3 : }
167 : }
168 : }
169 502 : if(checkForces && value->forcesWereAdded()) {
170 896 : for(std::size_t element=0; element<numberOfValues; ++element) {
171 : const double current=value->getForce(element);
172 671 : if(std::isfinite(current)) {
173 : continue;
174 : }
175 3 : std::ostringstream message;
176 : message << std::setprecision(17)
177 3 : << "PLUMED action non-finite trace: step=" << step
178 : << " phase=" << phase
179 : << " after_action=" << trigger.getLabel()
180 : << " after_action_type=" << trigger.getName()
181 : << " target=" << value->getName()
182 12 : << " field=force element=" << element
183 3 : << " value=" << current;
184 6 : plumed_merror(message.str());
185 3 : }
186 : }
187 : }
188 : }
189 : }
190 :
191 :
192 7368563 : void traceActionNonFinite(const ActionSet& actions, const Action& trigger,
193 : const long long step, const char* phase,
194 : const bool checkValues, const bool checkForces,
195 : Communicator& comm) {
196 7368563 : if(nonFiniteActionTraceLabels().empty()) {
197 7368376 : return;
198 : }
199 : std::string message;
200 : try {
201 187 : traceActionNonFiniteLocal(actions,trigger,step,phase,checkValues,checkForces);
202 6 : } catch(const std::exception& error) {
203 6 : message=error.what();
204 6 : }
205 : // A derivative or force may be non-finite on only one spatial rank.
206 : // Report the first failing rank to every participant before throwing.
207 187 : int reporter=message.empty() ? comm.Get_size() : comm.Get_rank();
208 187 : comm.Min(reporter);
209 187 : if(reporter<comm.Get_size()) {
210 9 : unsigned length=message.size();
211 9 : comm.Bcast(length,reporter);
212 9 : message.resize(length);
213 9 : comm.Bcast(message,reporter);
214 36 : plumed_merror("non-finite action trace on MPI rank="+std::to_string(reporter)+": "+message);
215 : }
216 : }
217 :
218 : }
219 :
220 : /// Small utility just used in this file to throw arbitrary exceptions
221 61 : [[noreturn]] static void testThrow(const char* what) {
222 61 : auto words=Tools::getWords(what);
223 61 : plumed_assert(words.size()>0);
224 : #define __PLUMED_THROW_NOMSG(type) if(words[0]==#type) throw type()
225 : #define __PLUMED_THROW_MSG(type) if(words[0]==#type) throw type(what)
226 62 : __PLUMED_THROW_MSG(PLMD::ExceptionError);
227 61 : __PLUMED_THROW_MSG(PLMD::ExceptionDebug);
228 60 : __PLUMED_THROW_MSG(PLMD::Exception);
229 59 : __PLUMED_THROW_MSG(PLMD::lepton::Exception);
230 57 : __PLUMED_THROW_NOMSG(std::bad_exception);
231 56 : __PLUMED_THROW_NOMSG(std::bad_array_new_length);
232 55 : __PLUMED_THROW_NOMSG(std::bad_alloc);
233 54 : __PLUMED_THROW_NOMSG(std::bad_function_call);
234 53 : __PLUMED_THROW_NOMSG(std::bad_weak_ptr);
235 52 : __PLUMED_THROW_NOMSG(std::bad_cast);
236 51 : __PLUMED_THROW_NOMSG(std::bad_typeid);
237 50 : __PLUMED_THROW_NOMSG(std::bad_variant_access);
238 49 : __PLUMED_THROW_NOMSG(std::bad_optional_access);
239 48 : __PLUMED_THROW_NOMSG(std::bad_any_cast);
240 47 : __PLUMED_THROW_MSG(std::underflow_error);
241 46 : __PLUMED_THROW_MSG(std::overflow_error);
242 45 : __PLUMED_THROW_MSG(std::range_error);
243 44 : __PLUMED_THROW_MSG(std::runtime_error);
244 43 : __PLUMED_THROW_MSG(std::out_of_range);
245 42 : __PLUMED_THROW_MSG(std::length_error);
246 41 : __PLUMED_THROW_MSG(std::domain_error);
247 40 : __PLUMED_THROW_MSG(std::invalid_argument);
248 39 : __PLUMED_THROW_MSG(std::logic_error);
249 :
250 :
251 :
252 38 : if(words[0]=="std::system_error") {
253 4 : plumed_assert(words.size()>2);
254 : int error_code;
255 4 : Tools::convert(words[2],error_code);
256 4 : if(words[1]=="std::generic_category") {
257 1 : throw std::system_error(error_code,std::generic_category(),what);
258 : }
259 3 : if(words[1]=="std::system_category") {
260 1 : throw std::system_error(error_code,std::system_category(),what);
261 : }
262 2 : if(words[1]=="std::iostream_category") {
263 1 : throw std::system_error(error_code,std::iostream_category(),what);
264 : }
265 1 : if(words[1]=="std::future_category") {
266 1 : throw std::system_error(error_code,std::future_category(),what);
267 : }
268 : }
269 :
270 34 : if(words[0]=="std::filesystem::filesystem_error") {
271 : int error_code;
272 3 : plumed_assert(words.size()>2);
273 3 : Tools::convert(words[2],error_code);
274 : std::error_code x_error_code;
275 3 : if(words[1]=="std::generic_category") {
276 3 : x_error_code=::std::error_code(error_code,::std::generic_category());
277 : }
278 3 : if(words[1]=="std::system_category") {
279 0 : x_error_code=::std::error_code(error_code,::std::system_category());
280 : }
281 3 : if(words[1]=="std::iostream_category") {
282 0 : x_error_code=::std::error_code(error_code,::std::iostream_category());
283 : }
284 3 : if(words[1]=="std::future_category") {
285 0 : x_error_code=::std::error_code(error_code,::std::future_category());
286 : }
287 :
288 3 : if(words.size()<4) {
289 2 : throw std::filesystem::filesystem_error(what,x_error_code);
290 : }
291 2 : if(words.size()<5) {
292 3 : throw std::filesystem::filesystem_error(what,std::filesystem::path(words[3]),x_error_code);
293 : }
294 4 : throw std::filesystem::filesystem_error(what,std::filesystem::path(words[3]),std::filesystem::path(words[4]),x_error_code);
295 : }
296 :
297 : #define __PLUMED_THROW_REGEX(name) if(words[1]=="std::regex_constants::error_" #name) throw std::regex_error(std::regex_constants::error_ ##name)
298 31 : if(words[0]=="std::regex_error") {
299 13 : plumed_assert(words.size()>1);
300 13 : __PLUMED_THROW_REGEX(collate);
301 12 : __PLUMED_THROW_REGEX(ctype);
302 11 : __PLUMED_THROW_REGEX(escape);
303 10 : __PLUMED_THROW_REGEX(backref);
304 9 : __PLUMED_THROW_REGEX(brack);
305 8 : __PLUMED_THROW_REGEX(paren);
306 7 : __PLUMED_THROW_REGEX(brace);
307 6 : __PLUMED_THROW_REGEX(badbrace);
308 5 : __PLUMED_THROW_REGEX(range);
309 4 : __PLUMED_THROW_REGEX(space);
310 3 : __PLUMED_THROW_REGEX(badrepeat);
311 2 : __PLUMED_THROW_REGEX(complexity);
312 1 : __PLUMED_THROW_REGEX(stack);
313 : }
314 :
315 : #define __PLUMED_THROW_FUTURE(name) if(words[1]=="std::future_errc::" #name) throw std::future_error(::std::future_errc::name)
316 18 : if(words[0]=="std::future_error") {
317 4 : plumed_assert(words.size()>1);
318 4 : __PLUMED_THROW_FUTURE(broken_promise);
319 3 : __PLUMED_THROW_FUTURE(future_already_retrieved);
320 2 : __PLUMED_THROW_FUTURE(promise_already_satisfied);
321 1 : __PLUMED_THROW_FUTURE(no_state);
322 : }
323 :
324 14 : if(words[0]=="std::ios_base::failure") {
325 1 : int error_code=0;
326 1 : if(words.size()>2) {
327 0 : Tools::convert(words[2],error_code);
328 : }
329 1 : if(words.size()>1 && words[1]=="std::generic_category") {
330 0 : throw std::ios_base::failure(what,std::error_code(error_code,std::generic_category()));
331 : }
332 1 : if(words.size()>1 && words[1]=="std::system_category") {
333 0 : throw std::ios_base::failure(what,std::error_code(error_code,std::system_category()));
334 : }
335 1 : if(words.size()>1 && words[1]=="std::iostream_category") {
336 0 : throw std::ios_base::failure(what,std::error_code(error_code,std::iostream_category()));
337 : }
338 1 : if(words.size()>1 && words[1]=="std::future_category") {
339 0 : throw std::ios_base::failure(what,std::error_code(error_code,std::future_category()));
340 : }
341 1 : throw std::ios_base::failure(what);
342 : }
343 :
344 13 : if(words[0]=="int") {
345 1 : int value=0;
346 1 : if(words.size()>1) {
347 0 : Tools::convert(words[1],value);
348 : }
349 1 : throw value;
350 : }
351 :
352 12 : if(words[0]=="test_nested1") {
353 : try {
354 18 : throw Exception(std::string("inner ")+what);
355 6 : } catch(...) {
356 : try {
357 12 : std::throw_with_nested(Exception(std::string("middle ")+what));
358 6 : } catch(...) {
359 12 : std::throw_with_nested(Exception(std::string("outer ")+what));
360 6 : }
361 6 : }
362 : }
363 :
364 6 : if(words[0]=="test_nested2") {
365 : try {
366 3 : throw std::bad_alloc();
367 3 : } catch(...) {
368 : try {
369 6 : std::throw_with_nested(Exception(std::string("middle ")+what));
370 3 : } catch(...) {
371 6 : std::throw_with_nested(Exception(std::string("outer ")+what));
372 3 : }
373 3 : }
374 : }
375 :
376 3 : if(words[0]=="test_nested3") {
377 : try {
378 2 : throw "inner";
379 2 : } catch(...) {
380 : try {
381 4 : std::throw_with_nested(Exception(std::string("middle ")+what));
382 2 : } catch(...) {
383 4 : std::throw_with_nested(Exception(std::string("outer ")+what));
384 2 : }
385 2 : }
386 : }
387 :
388 2 : plumed_error() << "unknown exception " << what;
389 61 : }
390 :
391 : namespace {
392 : /// This is an internal tool used to count how many PlumedMain objects have been created
393 : /// and if they were correctly destroyed.
394 : /// When using debug options, it leads to a crash
395 : /// Otherwise, it just prints a message
396 : class CountInstances {
397 : std::atomic<int> counter{};
398 : // private constructor to avoid direct usage
399 7264 : CountInstances() noexcept {}
400 7264 : ~CountInstances() {
401 7264 : if(counter!=0) {
402 0 : std::cerr<<"WARNING: internal inconsistency in allocated PlumedMain instances (" <<counter<< ")\n";
403 0 : std::cerr<<"Might be a consequence of incorrectly paired plumed_create/plumed_finalize in the C interface\n";
404 0 : std::cerr<<"Or it could be due to incorrect calls to std::exit, without properly destroying all PlumedMain objects\n";
405 : #ifndef NDEBUG
406 : std::cerr<<"This is a debug build, so the warning will make PLUMED abort\n";
407 : std::abort();
408 : #endif
409 : }
410 7264 : }
411 1619248 : static CountInstances & instance() {
412 1619248 : static CountInstances counter;
413 1619248 : return counter;
414 : }
415 : public:
416 : /// Only access through these static functions
417 : /// The first call to increase() ensures the instance is constructed
418 : /// This should provide the correct construction and destruction order
419 : /// also in cases where the PlumedMain object is constructed in the
420 : /// constructor of a static object
421 809624 : static void increase() noexcept {
422 809624 : ++instance().counter;
423 809624 : }
424 : /// See increase()
425 809624 : static void decrease() noexcept {
426 809624 : --instance().counter;
427 809624 : }
428 : };
429 :
430 : }
431 :
432 :
433 809624 : PlumedMain::PlumedMain():
434 809624 : datoms_fwd(*this),
435 : // automatically write on log in destructor
436 809624 : stopwatch_fwd(log),
437 809624 : actionSet_fwd(*this),
438 809624 : passtools(DataPassingTools::create(sizeof(double))) {
439 809624 : passtools->usingNaturalUnits=false;
440 809624 : increaseReferenceCounter();
441 809624 : log.link(comm);
442 809624 : log.setLinePrefix("PLUMED: ");
443 : // this is at last so as to avoid inconsistencies if an exception is thrown
444 809624 : CountInstances::increase(); // noexcept
445 809624 : }
446 :
447 : // destructor needed to delete forward declarated objects
448 1617864 : PlumedMain::~PlumedMain() {
449 809624 : CountInstances::decrease();
450 4046736 : }
451 :
452 : /////////////////////////////////////////////////////////////
453 : // MAIN INTERPRETER
454 :
455 : #define CHECK_INIT(ini,word) plumed_assert(ini)<<"cmd(\"" << word << "\") should be only used after plumed initialization"
456 : #define CHECK_NOTINIT(ini,word) plumed_assert(!(ini))<<"cmd(\"" << word << "\") should be only used before plumed initialization"
457 : #define CHECK_NOTNULL(val,word) plumed_assert(val)<<"NULL pointer received in cmd(\"" << word << "\")"
458 :
459 :
460 2188539 : void PlumedMain::cmd(std::string_view word,const TypesafePtr & val) {
461 :
462 : // Enumerate all possible commands:
463 : enum {
464 : #include "PlumedMainEnum.inc"
465 : };
466 :
467 : // Static object (initialized once) containing the map of commands:
468 : const static Tools::FastStringUnorderedMap<int> word_map = {
469 : #include "PlumedMainMap.inc"
470 2188539 : };
471 :
472 : try {
473 :
474 2188539 : auto ss=stopwatch.startPause();
475 :
476 : gch::small_vector<std::string_view> words;
477 2188539 : Tools::getWordsSimple(words,word);
478 :
479 2188539 : unsigned nw=words.size();
480 2188539 : if(nw==0) {
481 : // do nothing
482 : } else {
483 : int iword=-1;
484 : const auto it=word_map.find(words[0]);
485 2188539 : if(it!=word_map.end()) {
486 2188531 : iword=it->second;
487 : }
488 :
489 2188531 : switch(iword) {
490 161988 : case cmd_setBox:
491 161988 : CHECK_INIT(initialized,word);
492 161988 : CHECK_NOTNULL(val,word);
493 161988 : setInputValue( "Box", 0, 1, val );
494 161987 : break;
495 166467 : case cmd_setPositions:
496 166467 : CHECK_INIT(initialized,word);
497 166467 : setInputValue("posx", 0, 3, val );
498 166467 : setInputValue("posy", 1, 3, val );
499 166467 : setInputValue("posz", 2, 3, val );
500 166467 : break;
501 166429 : case cmd_setMasses:
502 166429 : CHECK_INIT(initialized,word);
503 166429 : setInputValue("Masses", 0, 1, val );
504 166426 : break;
505 88879 : case cmd_setCharges:
506 88879 : CHECK_INIT(initialized,word);
507 88879 : setInputValue("Charges", 0, 1, val);
508 88879 : break;
509 1 : case cmd_setPositionsX:
510 1 : CHECK_INIT(initialized,word);
511 1 : setInputValue("posx", 0, 1, val);
512 1 : break;
513 1 : case cmd_setPositionsY:
514 1 : CHECK_INIT(initialized,word);
515 1 : setInputValue("posy", 0, 1, val);
516 1 : break;
517 1 : case cmd_setPositionsZ:
518 1 : CHECK_INIT(initialized,word);
519 1 : setInputValue("posz", 0, 1, val);
520 1 : break;
521 158152 : case cmd_setVirial:
522 158152 : CHECK_INIT(initialized,word);
523 158152 : CHECK_NOTNULL(val,word);
524 158152 : setInputForce("Box",val);
525 158134 : break;
526 8294 : case cmd_setEnergy:
527 8294 : CHECK_INIT(initialized,word);
528 8294 : CHECK_NOTNULL(val,word);
529 8294 : if( name_of_energy!="" ) {
530 4089 : setInputValue( name_of_energy, 0, 1, val );
531 : }
532 : break;
533 166431 : case cmd_setForces:
534 166431 : CHECK_INIT(initialized,word);
535 166431 : setInputForce("posx",val);
536 166431 : setInputForce("posy",val);
537 166431 : setInputForce("posz",val);
538 166431 : break;
539 1 : case cmd_setForcesX:
540 1 : CHECK_INIT(initialized,word);
541 1 : setInputForce("posx",val);
542 1 : break;
543 1 : case cmd_setForcesY:
544 1 : CHECK_INIT(initialized,word);
545 1 : setInputForce("posy",val);
546 1 : break;
547 1 : case cmd_setForcesZ:
548 1 : CHECK_INIT(initialized,word);
549 1 : setInputForce("posz",val);
550 1 : break;
551 379230 : case cmd_calc:
552 379230 : CHECK_INIT(initialized,word);
553 379230 : calc();
554 : break;
555 123 : case cmd_prepareDependencies:
556 123 : CHECK_INIT(initialized,word);
557 123 : prepareDependencies();
558 : break;
559 108 : case cmd_shareData:
560 108 : CHECK_INIT(initialized,word);
561 108 : shareData();
562 : break;
563 14331 : case cmd_prepareCalc:
564 14331 : CHECK_INIT(initialized,word);
565 14331 : prepareCalc();
566 : break;
567 2517 : case cmd_performCalc:
568 2517 : CHECK_INIT(initialized,word);
569 2517 : performCalc();
570 : break;
571 11907 : case cmd_performCalcNoUpdate:
572 11907 : CHECK_INIT(initialized,word);
573 11907 : performCalcNoUpdate();
574 : break;
575 10 : case cmd_performCalcNoForces:
576 10 : CHECK_INIT(initialized,word);
577 10 : performCalcNoForces();
578 : break;
579 503 : case cmd_update:
580 503 : CHECK_INIT(initialized,word);
581 503 : update();
582 : break;
583 100451 : case cmd_setStep:
584 100451 : CHECK_INIT(initialized,word);
585 100454 : CHECK_NOTNULL(val,word);
586 100448 : step=val.get<int>();
587 100447 : startStep();
588 : break;
589 0 : case cmd_setStepLong:
590 0 : CHECK_INIT(initialized,word);
591 0 : CHECK_NOTNULL(val,word);
592 0 : step=val.get<long int>();
593 0 : startStep();
594 : break;
595 293068 : case cmd_setStepLongLong:
596 293068 : CHECK_INIT(initialized,word);
597 293068 : CHECK_NOTNULL(val,word);
598 293068 : step=val.get<long long int>();
599 293068 : startStep();
600 : break;
601 52953 : case cmd_setValue: {
602 52953 : CHECK_INIT(initialized,words[0]);
603 52953 : plumed_assert(nw==2);
604 105906 : setInputValue( std::string(words[1]), 0, 1, val );
605 : }
606 52953 : break;
607 : /* ADDED WITH API=7 */
608 0 : case cmd_setValueForces: {
609 0 : CHECK_INIT(initialized,words[0]);
610 0 : plumed_assert(nw==2);
611 0 : setInputForce( std::string(words[1]), val );
612 : }
613 0 : break;
614 : // words used less frequently:
615 1448 : case cmd_setAtomsNlocal:
616 1448 : CHECK_INIT(initialized,word);
617 1448 : CHECK_NOTNULL(val,word);
618 5809 : for(const auto & pp : inputs ) {
619 4361 : plumed_assert(pp);
620 4361 : DomainDecomposition* dd=pp->castToDomainDecomposition();
621 4361 : if( dd ) {
622 1448 : dd->setAtomsNlocal(val.get<int>());
623 : }
624 : }
625 : break;
626 1294 : case cmd_setAtomsGatindex:
627 1294 : CHECK_INIT(initialized,word);
628 5195 : for(const auto & pp : inputs ) {
629 3901 : plumed_assert(pp);
630 3901 : DomainDecomposition* dd=pp->castToDomainDecomposition();
631 3901 : if( dd ) {
632 1294 : dd->setAtomsGatindex(val,false);
633 : }
634 : }
635 : break;
636 2 : case cmd_setAtomsFGatindex:
637 2 : CHECK_INIT(initialized,word);
638 6 : for(const auto & pp : inputs ) {
639 4 : plumed_assert(pp);
640 4 : DomainDecomposition* dd=pp->castToDomainDecomposition();
641 4 : if( dd ) {
642 2 : dd->setAtomsGatindex(val,true);
643 : }
644 : }
645 : break;
646 152 : case cmd_setAtomsContiguous:
647 152 : CHECK_INIT(initialized,word);
648 152 : CHECK_NOTNULL(val,word);
649 608 : for(const auto & pp : inputs ) {
650 456 : plumed_assert(pp);
651 456 : DomainDecomposition* dd=pp->castToDomainDecomposition();
652 456 : if( dd ) {
653 152 : dd->setAtomsContiguous(val.get<int>());
654 : }
655 : }
656 : break;
657 140 : case cmd_createFullList:
658 140 : CHECK_INIT(initialized,word);
659 140 : CHECK_NOTNULL(val,word);
660 649 : for(const auto & pp : inputs ) {
661 509 : plumed_assert(pp);
662 509 : DomainDecomposition* dd=pp->castToDomainDecomposition();
663 509 : if( dd ) {
664 140 : dd->createFullList(val);
665 : }
666 : }
667 : break;
668 140 : case cmd_getFullList: {
669 140 : CHECK_INIT(initialized,word);
670 140 : CHECK_NOTNULL(val,word);
671 : unsigned nlists=0;
672 649 : for(const auto & pp : inputs ) {
673 509 : plumed_assert(pp);
674 509 : DomainDecomposition* dd=pp->castToDomainDecomposition();
675 509 : if( dd ) {
676 140 : dd->getFullList(val);
677 140 : nlists++;
678 : }
679 : }
680 140 : plumed_assert( nlists==1 );
681 : }
682 : break;
683 116 : case cmd_clearFullList:
684 116 : CHECK_INIT(initialized,word);
685 553 : for(const auto & pp : inputs ) {
686 437 : plumed_assert(pp);
687 437 : DomainDecomposition* dd=pp->castToDomainDecomposition();
688 437 : if( dd ) {
689 116 : dd->clearFullList();
690 : }
691 : }
692 : break;
693 : /* ADDED WITH API==6 */
694 124 : case cmd_getDataRank: {
695 124 : CHECK_INIT(initialized,words[0]);
696 124 : plumed_assert(nw==2 || nw==3);
697 124 : std::string vtype="";
698 124 : if( nw==3 ) {
699 0 : vtype=" TYPE="+std::string(words[2]);
700 : }
701 496 : readInputLine( "grab_" + std::string(words[1]) + ": GET ARG=" + std::string(words[1]) + vtype );
702 248 : ActionToGetData* as=actionSet.selectWithLabel<ActionToGetData*>("grab_"+std::string(words[1]));
703 124 : plumed_assert( as );
704 124 : as->get_rank( val );
705 : }
706 124 : break;
707 : /* ADDED WITH API==6 */
708 60 : case cmd_getDataShape: {
709 60 : CHECK_INIT(initialized,std::string(words[0]));
710 120 : ActionToGetData* as1=actionSet.selectWithLabel<ActionToGetData*>("grab_"+std::string(words[1]));
711 60 : plumed_assert( as1 );
712 60 : as1->get_shape( val );
713 : }
714 : break;
715 : /* ADDED WITH API==6 */
716 124 : case cmd_setMemoryForData: {
717 124 : CHECK_INIT(initialized,words[0]);
718 124 : plumed_assert(nw==2 || nw==3);
719 248 : ActionToGetData* as2=actionSet.selectWithLabel<ActionToGetData*>("grab_"+std::string(words[1]));
720 124 : plumed_assert( as2 );
721 124 : as2->set_memory( val );
722 : }
723 : break;
724 : /* ADDED WITH API==6 */
725 : case cmd_setErrorHandler: {
726 0 : if(val) {
727 0 : error_handler=*static_cast<const plumed_error_handler*>(val.get<const void*>());
728 : } else {
729 0 : error_handler.handler=NULL;
730 : }
731 : }
732 : break;
733 0 : case cmd_read:
734 0 : CHECK_INIT(initialized,word);
735 0 : if(val) {
736 0 : readInputFile(val.getCString());
737 : } else {
738 0 : readInputFile("plumed.dat");
739 : }
740 : break;
741 3794 : case cmd_readInputLine:
742 3794 : CHECK_INIT(initialized,word);
743 3794 : CHECK_NOTNULL(val,word);
744 3794 : readInputLine(val.getCString());
745 3656 : break;
746 26 : case cmd_readInputLines:
747 26 : CHECK_INIT(initialized,word);
748 26 : CHECK_NOTNULL(val,word);
749 26 : readInputLines(val.getCString());
750 26 : break;
751 3 : case cmd_clear: {
752 3 : CHECK_INIT(initialized,word);
753 3 : std::vector<int> natoms;
754 7 : for(const auto & pp : inputs ) {
755 4 : plumed_assert(pp);
756 4 : DomainDecomposition* dd=pp->castToDomainDecomposition();
757 4 : if ( dd ) {
758 1 : natoms.push_back( dd->getNumberOfAtoms() );
759 : }
760 : }
761 3 : actionSet.clearDelete();
762 3 : inputs.clear();
763 4 : for(unsigned i=0; i<natoms.size(); ++i) {
764 : std::string str_natoms;
765 1 : Tools::convert( natoms[i], str_natoms );
766 2 : readInputLine( MDEngine + ": DOMAIN_DECOMPOSITION NATOMS=" + str_natoms +
767 1 : " VALUE1=posx UNIT1=length PERIODIC1=NO CONSTANT1=False ROLE1=x" +
768 1 : " VALUE2=posy UNIT2=length PERIODIC2=NO CONSTANT2=False ROLE2=y" +
769 1 : " VALUE3=posz UNIT3=length PERIODIC3=NO CONSTANT3=False ROLE3=z" +
770 1 : " VALUE4=Masses UNIT4=mass PERIODIC4=NO CONSTANT4=True ROLE4=m" +
771 : " VALUE5=Charges UNIT5=charge PERIODIC5=NO CONSTANT5=True ROLE5=q");
772 :
773 : }
774 3 : setUnits( passtools->usingNaturalUnits, passtools->units );
775 3 : }
776 3 : break;
777 : case cmd_getApiVersion:
778 50 : CHECK_NOTNULL(val,word);
779 50 : val.set(int(10));
780 : break;
781 : // commands which can be used only before initialization:
782 2285 : case cmd_init:
783 2285 : CHECK_NOTINIT(initialized,word);
784 2285 : init();
785 : break;
786 1518 : case cmd_setRealPrecision:
787 1518 : CHECK_NOTINIT(initialized,word);
788 1518 : CHECK_NOTNULL(val,word);
789 3035 : passtools=DataPassingTools::create(val.get<int>());
790 1517 : passtools->usingNaturalUnits=false;
791 1517 : break;
792 1347 : case cmd_setMDLengthUnits:
793 1347 : CHECK_NOTINIT(initialized,word);
794 1347 : CHECK_NOTNULL(val,word);
795 1347 : passtools->MDUnits.setLength(passtools->MD2double(val));
796 : break;
797 1347 : case cmd_setMDChargeUnits:
798 1347 : CHECK_NOTINIT(initialized,word);
799 1347 : CHECK_NOTNULL(val,word);
800 1347 : passtools->MDUnits.setCharge(passtools->MD2double(val));
801 : break;
802 1347 : case cmd_setMDMassUnits:
803 1347 : CHECK_NOTINIT(initialized,word);
804 1347 : CHECK_NOTNULL(val,word);
805 1347 : passtools->MDUnits.setMass(passtools->MD2double(val));
806 : break;
807 45 : case cmd_setMDEnergyUnits:
808 45 : CHECK_NOTINIT(initialized,word);
809 45 : CHECK_NOTNULL(val,word);
810 45 : passtools->MDUnits.setEnergy(passtools->MD2double(val));
811 : break;
812 6 : case cmd_setMDTimeUnits:
813 6 : CHECK_NOTINIT(initialized,word);
814 6 : CHECK_NOTNULL(val,word);
815 6 : passtools->MDUnits.setTime(passtools->MD2double(val));
816 : break;
817 0 : case cmd_setNaturalUnits:
818 : // set the boltzman constant for MD in natural units (kb=1)
819 : // only needed in LJ codes if the MD is passing temperatures to plumed (so, not yet...)
820 : // use as cmd("setNaturalUnits")
821 0 : CHECK_NOTINIT(initialized,word);
822 0 : passtools->usingNaturalUnits=true;
823 0 : break;
824 75 : case cmd_setNoVirial: {
825 75 : CHECK_NOTINIT(initialized,word);
826 75 : ActionToPutData* ap=actionSet.selectWithLabel<ActionToPutData*>("Box");
827 75 : if( ap ) {
828 56 : ap->noforce=true;
829 : } else {
830 19 : ActionForInterface* af = actionSet.selectWithLabel<ActionForInterface*>(MDEngine);
831 19 : if( af ) {
832 0 : plumed_merror("setNoVirial should be called after number of atoms have been set");
833 : }
834 : }
835 : }
836 : break;
837 1401 : case cmd_setPlumedDat:
838 1401 : CHECK_NOTINIT(initialized,word);
839 1401 : CHECK_NOTNULL(val,word);
840 1401 : plumedDat=val.getCString();
841 : break;
842 1015 : case cmd_setMPIComm:
843 1015 : CHECK_NOTINIT(initialized,word);
844 1015 : comm.Set_comm(val);
845 1017 : for(const auto & pp : inputs ) {
846 2 : pp->Set_comm(comm);
847 : }
848 : break;
849 0 : case cmd_setMPIFComm:
850 0 : CHECK_NOTINIT(initialized,word);
851 0 : comm.Set_fcomm(val);
852 0 : for(const auto & pp : inputs ) {
853 0 : pp->Set_comm(comm);
854 : }
855 : break;
856 0 : case cmd_setMPImultiSimComm:
857 0 : CHECK_NOTINIT(initialized,word);
858 0 : multi_sim_comm.Set_comm(val);
859 : break;
860 2267 : case cmd_setNatoms: {
861 2267 : CHECK_NOTINIT(initialized,word);
862 2267 : CHECK_NOTNULL(val,word);
863 2267 : int natoms = val.get<int>();
864 : std::string str_natoms;
865 2257 : Tools::convert( natoms, str_natoms );
866 2257 : ActionForInterface* dd=actionSet.selectWithLabel<ActionForInterface*>(MDEngine);
867 2257 : if( !dd && natoms>0 )
868 4350 : readInputLine( MDEngine + ": DOMAIN_DECOMPOSITION NATOMS=" + str_natoms + +
869 2175 : " VALUE1=posx UNIT1=length PERIODIC1=NO CONSTANT1=False ROLE1=x" +
870 2175 : " VALUE2=posy UNIT2=length PERIODIC2=NO CONSTANT2=False ROLE2=y" +
871 2175 : " VALUE3=posz UNIT3=length PERIODIC3=NO CONSTANT3=False ROLE3=z" +
872 2175 : " VALUE4=Masses UNIT4=mass PERIODIC4=NO CONSTANT4=True ROLE4=m" +
873 4350 : " VALUE5=Charges UNIT5=charge PERIODIC5=NO CONSTANT5=True ROLE5=q", true );
874 : }
875 2257 : break;
876 1983 : case cmd_setTimestep: {
877 1983 : CHECK_NOTINIT(initialized,word);
878 1983 : CHECK_NOTNULL(val,word);
879 1983 : ActionToPutData* ts = actionSet.selectWithLabel<ActionToPutData*>("timestep");
880 1983 : if( !ts ) {
881 1981 : readInputLine("timestep: PUT UNIT=time PERIODIC=NO CONSTANT ROUND_TO_BASE_TEN", true);
882 3962 : ts = actionSet.selectWithLabel<ActionToPutData*>("timestep");
883 : }
884 3966 : if( !ts->setValuePointer("timestep", val ) ) {
885 0 : plumed_error();
886 : }
887 1983 : ts->updateUnits( passtools.get() );
888 : }
889 : break;
890 : /* ADDED WITH API==2 */
891 330 : case cmd_setKbT: {
892 330 : CHECK_NOTINIT(initialized,word);
893 330 : CHECK_NOTNULL(val,word);
894 330 : readInputLine("kBT: PUT CONSTANT PERIODIC=NO UNIT=energy", true);
895 330 : ActionToPutData* kb = actionSet.selectWithLabel<ActionToPutData*>("kBT");
896 660 : if( !kb->setValuePointer("kBT", val ) ) {
897 0 : plumed_error();
898 : }
899 330 : kb->updateUnits( passtools.get() );
900 : }
901 : break;
902 : /* ADDED WITH API==3 */
903 112 : case cmd_setRestart:
904 112 : CHECK_NOTINIT(initialized,word);
905 112 : CHECK_NOTNULL(val,word);
906 112 : if(val.get<int>()!=0) {
907 54 : restart=true;
908 : }
909 : break;
910 : /* ADDED WITH API==4 */
911 0 : case cmd_doCheckPoint:
912 0 : CHECK_INIT(initialized,word);
913 0 : CHECK_NOTNULL(val,word);
914 0 : doCheckPoint = false;
915 0 : if(val.get<int>()!=0) {
916 0 : doCheckPoint = true;
917 : }
918 : break;
919 : /* ADDED WITH API==6 */
920 : case cmd_setNumOMPthreads:
921 3 : CHECK_NOTNULL(val,word);
922 : {
923 3 : auto nt=val.get<unsigned>();
924 : if(nt==0) {
925 : nt=1;
926 : }
927 3 : OpenMP::setNumThreads(nt);
928 : }
929 : break;
930 : /* ADDED WITH API==10 */
931 : case cmd_setGpuDeviceId:
932 0 : CHECK_NOTNULL(val,word);
933 : {
934 0 : auto id=val.get<int>();
935 0 : if(id>=0) {
936 0 : gpuDeviceId=id;
937 : }
938 : }
939 : break;
940 : /* ADDED WITH API==6 */
941 : /* only used for testing */
942 : case cmd_throw:
943 61 : CHECK_NOTNULL(val,word);
944 61 : testThrow(val.getCString());
945 : /* ADDED WITH API==10 */
946 : case cmd_setNestedExceptions:
947 30 : CHECK_NOTNULL(val,word);
948 30 : if(val.get<int>()!=0) {
949 29 : nestedExceptions=true;
950 : } else {
951 1 : nestedExceptions=false;
952 : }
953 : break;
954 : /* STOP API */
955 1434 : case cmd_setMDEngine:
956 1434 : CHECK_NOTINIT(initialized,word);
957 1434 : CHECK_NOTNULL(val,word);
958 1434 : MDEngine=val.getCString();
959 : break;
960 1387 : case cmd_setLog:
961 1387 : CHECK_NOTINIT(initialized,word);
962 1387 : log.link(val.get<FILE*>());
963 : break;
964 227 : case cmd_setLogFile:
965 227 : CHECK_NOTINIT(initialized,word);
966 227 : CHECK_NOTNULL(val,word);
967 227 : log.open(val.getCString());
968 227 : break;
969 : // other commands that should be used after initialization:
970 286469 : case cmd_setStopFlag:
971 286469 : CHECK_INIT(initialized,word);
972 286469 : CHECK_NOTNULL(val,word);
973 286469 : val.get<int*>(); // just check type and discard pointer
974 286468 : stopFlag=val.copy();
975 286468 : break;
976 0 : case cmd_getExchangesFlag:
977 0 : CHECK_INIT(initialized,word);
978 0 : CHECK_NOTNULL(val,word);
979 0 : exchangePatterns.getFlag(*val.get<int*>()); // note: getFlag changes the value of the reference!
980 : break;
981 0 : case cmd_setExchangesSeed:
982 0 : CHECK_INIT(initialized,word);
983 0 : CHECK_NOTNULL(val,word);
984 0 : exchangePatterns.setSeed(val.get<int>());
985 : break;
986 0 : case cmd_setNumberOfReplicas:
987 0 : CHECK_INIT(initialized,word);
988 0 : CHECK_NOTNULL(val,word);
989 0 : exchangePatterns.setNofR(val.get<int>());
990 : break;
991 0 : case cmd_getExchangesList:
992 0 : CHECK_INIT(initialized,word);
993 0 : CHECK_NOTNULL(val,word);
994 0 : exchangePatterns.getList(val);
995 : break;
996 1270 : case cmd_runFinalJobs:
997 1270 : CHECK_INIT(initialized,word);
998 1270 : runJobsAtEndOfCalculation();
999 : break;
1000 96 : case cmd_isEnergyNeeded: {
1001 96 : CHECK_INIT(initialized,word);
1002 96 : CHECK_NOTNULL(val,word);
1003 96 : if( name_of_energy =="" ) {
1004 96 : val.set(int(0));
1005 : } else {
1006 0 : ActionToPutData* ap=actionSet.selectWithLabel<ActionToPutData*>(name_of_energy);
1007 0 : if(ap->isActive()) {
1008 0 : val.set(int(1));
1009 : } else {
1010 0 : val.set(int(0));
1011 : }
1012 : }
1013 : }
1014 : break;
1015 78057 : case cmd_getBias:
1016 78057 : CHECK_INIT(initialized,word);
1017 78057 : CHECK_NOTNULL(val,word);
1018 78057 : plumedQuantityToMD( "energy", getBias(), val );
1019 78057 : break;
1020 : case cmd_checkAction:
1021 2 : CHECK_NOTNULL(val,word);
1022 2 : plumed_assert(nw==2);
1023 3 : val.set(int(checkAction(std::string(words[1])) ? 1:0));
1024 2 : break;
1025 : case cmd_setExtraCV: {
1026 30 : CHECK_NOTNULL(val,word);
1027 30 : plumed_assert(nw==2);
1028 60 : if( valueExists(std::string(words[1])) ) {
1029 90 : setInputValue( std::string(words[1]), 0, 1, val );
1030 : }
1031 : }
1032 : break;
1033 : case cmd_setExtraCVForce: {
1034 30 : CHECK_NOTNULL(val,word);
1035 30 : plumed_assert(nw==2);
1036 60 : if( valueExists(std::string(words[1])) ) {
1037 60 : setInputForce( std::string(words[1]), val );
1038 : }
1039 : }
1040 : break;
1041 : /* ADDED WITH API==10 */
1042 : case cmd_isExtraCVNeeded:
1043 10 : CHECK_NOTNULL(val,word);
1044 10 : plumed_assert(nw==2);
1045 10 : val.set(int(0));
1046 56 : for(const auto & p : inputs) {
1047 50 : if( p->getLabel()==words[1] && p->isActive() ) {
1048 4 : val.set(int(1));
1049 : break;
1050 : }
1051 : }
1052 : break;
1053 2672 : case cmd_GREX:
1054 2672 : if(!grex) {
1055 802 : grex=Tools::make_unique<GREX>(*this);
1056 : }
1057 2672 : plumed_massert(grex,"error allocating grex");
1058 : {
1059 3060 : std::string kk=std::string(words[1]);
1060 2843 : for(unsigned i=2; i<words.size(); i++) {
1061 342 : kk+=" "+std::string(words[i]);
1062 : }
1063 2672 : grex->cmd(kk.c_str(),val);
1064 : }
1065 2672 : break;
1066 22298 : case cmd_CLTool:
1067 22298 : CHECK_NOTINIT(initialized,word);
1068 22298 : if(!cltool) {
1069 7283 : cltool=Tools::make_unique<CLToolMain>();
1070 : }
1071 : {
1072 22686 : std::string kk(words[1]);
1073 22298 : for(unsigned i=2; i<words.size(); i++) {
1074 0 : kk+=" "+std::string(words[i]);
1075 : }
1076 22298 : cltool->cmd(kk.c_str(),val);
1077 : }
1078 22298 : break;
1079 : break;
1080 : /* ADDED WITH API==7 */
1081 : case cmd_convert: {
1082 : double v;
1083 57 : plumed_assert(words.size()==2);
1084 502 : if(Tools::convertNoexcept(std::string(words[1]),v)) {
1085 47 : passtools->double2MD(v,val);
1086 : }
1087 : }
1088 57 : break;
1089 8 : default:
1090 16 : plumed_error() << "cannot interpret cmd(\"" << word << "\"). check plumed developers manual to see the available commands.";
1091 : break;
1092 : }
1093 : }
1094 :
1095 2188927 : } catch (...) {
1096 388 : if(log.isOpen()) {
1097 : try {
1098 314 : log<<"\n################################################################################\n";
1099 314 : log<<Tools::concatenateExceptionMessages();
1100 314 : log<<"\n################################################################################\n";
1101 314 : log.flush();
1102 0 : } catch(...) {
1103 : // ignore errors here.
1104 : // in any case, we are rethrowing this below
1105 0 : }
1106 : }
1107 388 : throw;
1108 388 : }
1109 2188151 : }
1110 :
1111 : ////////////////////////////////////////////////////////////////////////
1112 :
1113 3087 : void PlumedMain::updateReplicaIdentity(Communicator& intra) {
1114 : // The engine normally hands the inter-replica communicator to a single rank of
1115 : // each replica (see GREX setMPIIntercomm); on the other ranks multi_sim_comm keeps
1116 : // its MPI_COMM_SELF default and reports size 1, rank 0. Taking the maximum over
1117 : // the intra-replica communicator gives every rank the same, correct answer,
1118 : // whichever rank or ranks the engine chose. This is a collective over intra.
1119 : // This can run twice (GREX init and init), so clear the stand-in mark first.
1120 3087 : multi_sim_comm.setStandIn(false);
1121 3087 : const int localSize=multi_sim_comm.Get_size();
1122 3087 : nReplicas=localSize;
1123 3087 : replicaIndex=multi_sim_comm.Get_rank();
1124 3087 : intra.Max(nReplicas);
1125 3087 : intra.Max(replicaIndex);
1126 : // With more than one replica, a rank whose multi_sim_comm has size 1 was not
1127 : // given the inter-replica communicator. Mark it, so that using it there
1128 : // raises an error in debug builds instead of silently acting on one replica.
1129 5674 : multi_sim_comm.setStandIn(nReplicas>1 && localSize==1);
1130 3087 : }
1131 :
1132 2285 : void PlumedMain::init() {
1133 : // check that initialization just happens once
1134 2285 : initialized=true;
1135 2285 : if(!log.isOpen()) {
1136 675 : log.link(stdout);
1137 : }
1138 2285 : log<<"PLUMED is starting\n";
1139 6855 : log<<"Version: "<<config::getVersionLong()<<" (git: "<<config::getVersionGit()<<") "
1140 9140 : <<"compiled on " <<config::getCompilationDate() << " at " << config::getCompilationTime() << "\n";
1141 2285 : log<<"Please cite these papers when using PLUMED ";
1142 2285 : log<<cite("The PLUMED consortium, Nat. Methods 16, 670 (2019)");
1143 2285 : log<<cite("Tribello, Bonomi, Branduardi, Camilloni, and Bussi, Comput. Phys. Commun. 185, 604 (2014)");
1144 2285 : log<<"\n";
1145 2285 : log<<"For further information see the PLUMED web page at http://www.plumed.org\n";
1146 2285 : log<<"Root: "<<config::getPlumedRoot()<<"\n";
1147 2285 : log<<"LibraryPath: "<<config::getLibraryPath()<<"\n";
1148 4570 : log<<"For installed feature, see "<<config::getPlumedRoot() + "/src/config/config.txt\n";
1149 2285 : log.printf("Molecular dynamics engine: %s\n",MDEngine.c_str());
1150 2285 : log.printf("Precision of reals: %d\n",passtools->getRealPrecision());
1151 3828 : log.printf("Running over %d %s\n",comm.Get_size(),(comm.Get_size()>1?"nodes":"node"));
1152 2285 : log<<"Number of threads: "<<OpenMP::getNumThreads()<<"\n";
1153 2285 : log<<"Cache line size: "<<OpenMP::getCachelineSize()<<"\n";
1154 8980 : for(const auto & pp : inputs ) {
1155 6695 : plumed_assert(pp);
1156 6695 : DomainDecomposition* dd=pp->castToDomainDecomposition();
1157 6695 : if ( dd ) {
1158 2175 : log.printf("Number of atoms: %d\n",dd->getNumberOfAtoms());
1159 : }
1160 : }
1161 2285 : if(grex) {
1162 802 : log.printf("GROMACS-like replica exchange is on\n");
1163 : }
1164 2285 : log.printf("File suffix: %s\n",getSuffix().c_str());
1165 : // Before any action is created, so that their constructors can rely on it.
1166 2285 : updateReplicaIdentity(comm);
1167 2285 : if(plumedDat.length()>0) {
1168 1397 : readInputFile(plumedDat);
1169 : plumedDat="";
1170 : }
1171 2225 : setUnits( passtools->usingNaturalUnits, passtools->units );
1172 2225 : ActionToPutData* ts = actionSet.selectWithLabel<ActionToPutData*>("timestep");
1173 2225 : if(ts) {
1174 1917 : log.printf("Timestep: %f\n",(ts->copyOutput(0))->get());
1175 : }
1176 2225 : ActionToPutData* kb = actionSet.selectWithLabel<ActionToPutData*>("kBT");
1177 2225 : if(kb) {
1178 330 : log.printf("KbT: %f\n",(kb->copyOutput(0))->get());
1179 : } else {
1180 1895 : log.printf("KbT has not been set by the MD engine\n");
1181 1895 : log.printf("It should be set by hand where needed\n");
1182 : }
1183 2225 : log<<"Relevant bibliography:\n";
1184 2225 : log<<citations;
1185 2225 : log<<"Please read and cite where appropriate!\n";
1186 2225 : log<<"Finished setup\n";
1187 2225 : }
1188 :
1189 19026 : void PlumedMain::setupInterfaceActions() {
1190 19026 : inputs.clear();
1191 19026 : std::vector<ActionForInterface*> ap=actionSet.select<ActionForInterface*>();
1192 109534 : for(unsigned i=0; i<ap.size(); ++i) {
1193 90508 : if( ap[i]->getName()=="ENERGY" || ap[i]->getDependencies().size()==0 ) {
1194 35418 : inputs.push_back( ap[i] );
1195 : }
1196 : }
1197 19026 : }
1198 :
1199 1423 : void PlumedMain::readInputFile(const std::string & str) {
1200 1423 : plumed_assert(initialized);
1201 1423 : log<<"FILE: "<<str<<"\n";
1202 1423 : IFile ifile;
1203 1423 : ifile.link(*this);
1204 1423 : ifile.open(str);
1205 1422 : ifile.allowNoEOL();
1206 1422 : readInputFile(ifile);
1207 1362 : log<<"END FILE: "<<str<<"\n";
1208 1362 : log.flush();
1209 :
1210 1423 : }
1211 :
1212 1448 : void PlumedMain::readInputFile(IFile & ifile) {
1213 1448 : std::vector<std::string> words;
1214 22405 : while(Tools::getParsedLine(ifile,words) && !endPlumed) {
1215 21017 : readInputWords(words,false);
1216 : }
1217 1388 : endPlumed=false;
1218 1388 : pilots=actionSet.select<ActionPilot*>();
1219 1388 : setupInterfaceActions();
1220 1448 : }
1221 :
1222 21567 : void PlumedMain::readInputLine(const std::string & str, const bool& before_init) {
1223 21567 : if( !before_init ) {
1224 3992 : plumed_assert(initialized);
1225 : }
1226 21567 : if(str.empty()) {
1227 0 : return;
1228 : }
1229 21567 : std::vector<std::string> words=Tools::getWords(str);
1230 21567 : if( before_init ) {
1231 17575 : plumed_assert( citations.empty() );
1232 : }
1233 21567 : citations.clear();
1234 21567 : readInputWords(words,before_init);
1235 21431 : if(!citations.empty()) {
1236 107 : log<<"Relevant bibliography:\n";
1237 107 : log<<citations;
1238 107 : log<<"Please read and cite where appropriate!\n";
1239 : }
1240 21567 : }
1241 :
1242 26 : void PlumedMain::readInputLines(const std::string & str) {
1243 26 : plumed_assert(initialized);
1244 26 : if(str.empty()) {
1245 0 : return;
1246 : }
1247 :
1248 26 : log<<"FILE: (temporary)\n";
1249 :
1250 : // Open a temporary file
1251 26 : auto fp=std::tmpfile();
1252 26 : plumed_assert(fp);
1253 :
1254 : // make sure file is closed (and thus deleted) also if an exception occurs
1255 : auto deleter=[](auto FP) {
1256 26 : std::fclose(FP);
1257 : };
1258 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
1259 :
1260 26 : auto ret=std::fputs(str.c_str(),fp);
1261 26 : plumed_assert(ret!=EOF);
1262 :
1263 26 : std::rewind(fp);
1264 :
1265 26 : IFile ifile;
1266 26 : ifile.link(*this);
1267 26 : ifile.link(fp);
1268 26 : ifile.allowNoEOL();
1269 :
1270 26 : readInputFile(ifile);
1271 26 : log<<"END FILE: (temporary)\n";
1272 26 : }
1273 :
1274 66212 : void PlumedMain::readInputWords(const std::vector<std::string> & words, const bool& before_init) {
1275 66212 : if( !before_init ) {
1276 48637 : plumed_assert(initialized);
1277 : }
1278 66212 : if(words.empty()) {
1279 : return;
1280 66038 : } else if(words[0]=="_SET_SUFFIX") {
1281 3 : plumed_assert(words.size()==2);
1282 : setSuffix(words[1]);
1283 : } else {
1284 66035 : std::vector<std::string> interpreted(words);
1285 66035 : Tools::interpretLabel(interpreted);
1286 131846 : auto action=actionRegister().create(dlloader.getHandles(),ActionOptions(*this,interpreted));
1287 65811 : if(!action) {
1288 : std::string msg;
1289 : msg ="ERROR\nI cannot understand line:";
1290 0 : for(unsigned i=0; i<interpreted.size(); ++i) {
1291 0 : msg+=" "+interpreted[i];
1292 : }
1293 : msg+="\nMaybe a missing space or a typo?";
1294 0 : log << msg;
1295 0 : log.flush();
1296 0 : plumed_merror(msg);
1297 : }
1298 65811 : action->checkRead();
1299 65809 : actionSet.emplace_back(std::move(action));
1300 : // update the lists incrementally: rebuilding them by scanning the whole action set
1301 : // after every line makes reading the input quadratic in the number of actions
1302 65809 : Action* last=actionSet.back().get();
1303 65809 : if( auto* pilot=dynamic_cast<ActionPilot*>(last) ) {
1304 5037 : pilots.push_back(pilot);
1305 : }
1306 65809 : if( dynamic_cast<ActionForInterface*>(last) ) {
1307 17638 : setupInterfaceActions();
1308 : }
1309 66037 : };
1310 : }
1311 :
1312 : ////////////////////////////////////////////////////////////////////////
1313 :
1314 0 : void PlumedMain::exit(int c) {
1315 0 : comm.Abort(c);
1316 0 : }
1317 :
1318 66033 : Log& PlumedMain::getLog() {
1319 66033 : return log;
1320 : }
1321 :
1322 379230 : void PlumedMain::calc() {
1323 379230 : prepareCalc();
1324 379211 : performCalc();
1325 379148 : }
1326 :
1327 393561 : void PlumedMain::prepareCalc() {
1328 393561 : prepareDependencies();
1329 393561 : shareData();
1330 393542 : }
1331 :
1332 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
1333 : // here we have the main steps in "calc()"
1334 : // they can be called individually, but the standard thing is to
1335 : // traverse them in this order:
1336 393936 : void PlumedMain::prepareDependencies() {
1337 :
1338 : // Stopwatch is stopped when sw goes out of scope
1339 393936 : auto sw=stopwatch.startStop("1 Prepare dependencies");
1340 :
1341 : // activate all the actions which are on step
1342 : // activation is recursive and enables also the dependencies
1343 : // before doing that, the prepare() method is called to see if there is some
1344 : // new/changed dependency (up to now, only useful for dependences on virtual atoms,
1345 : // which can be dynamically changed).
1346 :
1347 : // First switch off all actions
1348 6172150 : for(const auto & p : actionSet) {
1349 5778214 : p->deactivate();
1350 : }
1351 :
1352 : // for optimization, an "active" flag remains false if no action at all is active
1353 393936 : active=false;
1354 2024775 : for(unsigned i=0; i<pilots.size(); ++i) {
1355 1630839 : if(pilots[i]->onStep()) {
1356 1525031 : pilots[i]->activate();
1357 1525031 : active=true;
1358 : }
1359 : };
1360 :
1361 : // This stops the driver calculation if there is not a read action
1362 393936 : if( !active && !inputsAreActive() ) {
1363 2 : stopFlag.set(int(1));
1364 : }
1365 :
1366 : // also, if one of them is the total energy, tell to atoms that energy should be collected
1367 6172150 : for(const auto & p : actionSet) {
1368 5778214 : if(p->isActive()) {
1369 3689728 : if(p->checkNeedsGradients()) {
1370 2030 : p->setOption("GRADIENTS");
1371 : }
1372 : }
1373 : }
1374 :
1375 393936 : }
1376 :
1377 1915 : bool PlumedMain::inputsAreActive() const {
1378 5498 : for(const auto & ip : inputs) {
1379 5441 : if( ip->onStep() ) {
1380 : return true;
1381 : }
1382 : }
1383 : return false;
1384 : }
1385 :
1386 126 : void PlumedMain::shareAll() {
1387 504 : for(const auto & ip : inputs) {
1388 378 : ip->shareAll();
1389 : }
1390 126 : }
1391 :
1392 393669 : void PlumedMain::shareData() {
1393 : // atom positions are shared (but only if there is something to do)
1394 393669 : if(!active) {
1395 1854 : return;
1396 : }
1397 : // Stopwatch is stopped when sw goes out of scope
1398 391815 : auto sw=stopwatch.startStop("2 Sharing data");
1399 1243175 : for(const auto & ip : inputs) {
1400 851379 : ip->share();
1401 : }
1402 391815 : }
1403 :
1404 11907 : void PlumedMain::performCalcNoUpdate() {
1405 11907 : waitData();
1406 11907 : justCalculate();
1407 11907 : backwardPropagate();
1408 11907 : resetInputs();
1409 11907 : }
1410 :
1411 10 : void PlumedMain::performCalcNoForces() {
1412 10 : waitData();
1413 10 : justCalculate();
1414 10 : }
1415 :
1416 381728 : void PlumedMain::performCalc() {
1417 381728 : waitData();
1418 381728 : justCalculate();
1419 381676 : backwardPropagate();
1420 381665 : update();
1421 381665 : resetInputs();
1422 381665 : }
1423 :
1424 393771 : void PlumedMain::waitData() {
1425 393771 : if(!active) {
1426 1854 : return;
1427 : }
1428 : // Stopwatch is stopped when sw goes out of scope
1429 391917 : auto sw=stopwatch.startStop("3 Waiting for data");
1430 1243626 : for(const auto & ip : inputs) {
1431 : // shareData also starts the first-step domain exchange for inactive
1432 : // inputs, to initialize fixed masses and charges. Complete that exchange
1433 : // before clearing firststep, or later atomic actions receive stale data.
1434 851709 : if( (ip->isActive() || ip->firststep) && ip->hasBeenSet() ) {
1435 383197 : ip->wait();
1436 468512 : } else if( ip->isActive() ) {
1437 8688 : ip->warning("input requested but this quantity has not been set");
1438 : }
1439 : }
1440 391917 : }
1441 :
1442 393771 : void PlumedMain::justCalculate() {
1443 393771 : if(!active) {
1444 1854 : return;
1445 : }
1446 : // Stopwatch is stopped when sw goes out of scope
1447 391917 : auto sw=stopwatch.startStop("4 Calculating (forward loop)");
1448 391917 : bias=0.0;
1449 391917 : work=0.0;
1450 :
1451 : // Check the input actions to determine if we need to calculate constants that
1452 : // depend on masses and charges
1453 : bool firststep=false;
1454 1243626 : for(const auto & ip : inputs) {
1455 851709 : if( ip->firststep ) {
1456 : firststep=true;
1457 : }
1458 : }
1459 391917 : if( firststep ) {
1460 7999 : for(const auto & ip : inputs) {
1461 6011 : ip->firststep=false;
1462 : }
1463 : }
1464 :
1465 : // the counter is incremented at the top of the loop so that the labels of the
1466 : // detailed timers stay aligned also after the actions skipped with continue
1467 : int iaction=-1;
1468 : // calculate the active actions in order (assuming *backward* dependence)
1469 6136928 : for(const auto & pp : actionSet) {
1470 5745063 : iaction++;
1471 : Action* p(pp.get());
1472 5745063 : plumed_assert(p);
1473 : try {
1474 5745063 : if(p->isActive()) {
1475 : // Stopwatch is stopped when sw goes out of scope.
1476 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
1477 3687496 : Stopwatch::Handler swh;
1478 3687496 : if(detailedTimers) {
1479 4074 : auto actionNumberLabel=std::to_string(iaction);
1480 4074 : const unsigned m=actionSet.size();
1481 : unsigned k=0;
1482 : unsigned n=1;
1483 12222 : while(n<m) {
1484 8148 : n*=10;
1485 8148 : k++;
1486 : }
1487 4074 : auto spaces=std::string(k-actionNumberLabel.length(),' ');
1488 12222 : swh=stopwatch.startStop("4A " + spaces + actionNumberLabel+" "+p->getLabel());
1489 : }
1490 3687496 : ActionWithValue*av=p->castToActionWithValue();
1491 3687496 : ActionAtomistic*aa=p->castToActionAtomistic();
1492 : {
1493 3687496 : if(av) {
1494 3354040 : av->clearInputForces();
1495 : }
1496 : if(av) {
1497 3354040 : av->clearDerivatives();
1498 : }
1499 3354040 : if( av && av->calculateOnUpdate() ) {
1500 : continue ;
1501 : }
1502 : }
1503 : {
1504 3683492 : if(aa)
1505 606519 : if(aa->isActive()) {
1506 606519 : aa->retrieveAtoms();
1507 : }
1508 : }
1509 3683492 : if(p->checkNumericalDerivatives()) {
1510 271 : p->calculateNumericalDerivatives();
1511 : } else {
1512 3683221 : p->calculate();
1513 : }
1514 3683446 : traceActionNonFinite(actionSet,*p,step,"forward",true,false,comm);
1515 :
1516 : // This retrieves components called bias
1517 3683440 : if(av) {
1518 3349984 : bias+=av->getOutputQuantity("bias");
1519 3349984 : work+=av->getOutputQuantity("work");
1520 3349984 : av->setGradientsIfNeeded();
1521 : }
1522 : // This makes all values that depend on the (fixed) masses and charges constant
1523 3683440 : if( firststep ) {
1524 38403 : p->setupConstantValues( true );
1525 : }
1526 3683440 : ActionWithVirtualAtom*avv=p->castToActionWithVirtualAtom();
1527 3683440 : if(avv) {
1528 39262 : avv->setGradientsIfNeeded();
1529 : }
1530 3687496 : }
1531 52 : } catch(...) {
1532 104 : plumed_error_nested() << "An error happened while calculating " << p->getLabel();
1533 52 : }
1534 : }
1535 391917 : }
1536 :
1537 0 : void PlumedMain::justApply() {
1538 0 : backwardPropagate();
1539 0 : update();
1540 0 : resetInputs();
1541 0 : }
1542 :
1543 393583 : void PlumedMain::backwardPropagate() {
1544 393583 : if(!active) {
1545 1854 : return;
1546 : }
1547 : int iaction=0;
1548 : // Stopwatch is stopped when sw goes out of scope
1549 391729 : auto sw=stopwatch.startStop("5 Applying (backward loop)");
1550 : // apply them in reverse order
1551 6133534 : for(auto pp=actionSet.rbegin(); pp!=actionSet.rend(); ++pp) {
1552 : const auto & p(pp->get());
1553 5741816 : if(p->isActive()) {
1554 :
1555 : // Stopwatch is stopped when sw goes out of scope.
1556 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
1557 3685125 : Stopwatch::Handler swh;
1558 3685125 : if(detailedTimers) {
1559 4074 : auto actionNumberLabel=std::to_string(iaction);
1560 4074 : const unsigned m=actionSet.size();
1561 : unsigned k=0;
1562 : unsigned n=1;
1563 12222 : while(n<m) {
1564 8148 : n*=10;
1565 8148 : k++;
1566 : }
1567 4074 : auto spaces=std::string(k-actionNumberLabel.length(),' ');
1568 12222 : swh=stopwatch.startStop("5A " + spaces + actionNumberLabel+" "+p->getLabel());
1569 : }
1570 :
1571 3685125 : p->apply();
1572 3685117 : traceActionNonFinite(actionSet,*p,step,"backward",false,true,comm);
1573 3685125 : }
1574 5741805 : iaction++;
1575 : }
1576 :
1577 : // Stopwatch is stopped when sw goes out of scope.
1578 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
1579 391718 : Stopwatch::Handler sw1;
1580 391718 : if(detailedTimers) {
1581 194 : sw1=stopwatch.startStop("5B Update forces");
1582 : }
1583 391729 : }
1584 :
1585 382168 : void PlumedMain::update() {
1586 382168 : if(!active) {
1587 1854 : return;
1588 : }
1589 :
1590 : // Stopwatch is stopped when sw goes out of scope
1591 380314 : auto sw=stopwatch.startStop("6 Update");
1592 :
1593 : // update step (for statistics, etc)
1594 380314 : updateFlags.push(true);
1595 5828347 : for(const auto & p : actionSet) {
1596 5448033 : p->beforeUpdate();
1597 5448033 : if(p->isActive() && p->checkUpdate() && updateFlagsTop()) {
1598 3525821 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(p.get());
1599 3525821 : if( av && av->calculateOnUpdate() ) {
1600 4004 : p->prepare();
1601 4004 : p->calculate();
1602 : } else {
1603 3521817 : p->update();
1604 : }
1605 : }
1606 : }
1607 : // an UPDATE_IF block that is not closed with END extends to the end of the input,
1608 : // so more than one flag can be left here. Extra END lines are rejected when reading the input.
1609 380314 : if(updateFlags.empty()) {
1610 0 : plumed_merror("non matching changes in the update flags");
1611 : }
1612 760632 : while(!updateFlags.empty()) {
1613 : updateFlags.pop();
1614 : }
1615 : // Check that no action has told the calculation to stop
1616 380314 : if(stopNow) {
1617 66 : if(stopFlag) {
1618 66 : stopFlag.set(int(1));
1619 : } else {
1620 0 : plumed_merror("your md code cannot handle plumed stop events - add a call to plumed.comm(stopFlag,stopCondition)");
1621 : }
1622 : }
1623 :
1624 : // flush by default every 10000 steps
1625 : // hopefully will not affect performance
1626 : // also if receive checkpointing signal
1627 380314 : if(step%10000==0||doCheckPoint) {
1628 3595 : fflush();
1629 3595 : log.flush();
1630 97884 : for(const auto & p : actionSet) {
1631 94289 : p->fflush();
1632 : }
1633 : }
1634 : // invalidate the checkpoint status
1635 : // this fix issues like #1428
1636 380314 : doCheckPoint = false;
1637 380314 : }
1638 :
1639 50 : void PlumedMain::load(const std::string& fileName) {
1640 50 : if(DLLoader::installed()) {
1641 50 : std::string libName=fileName;
1642 50 : size_t n=libName.find_last_of(".");
1643 51 : std::string extension="";
1644 50 : std::string base=libName;
1645 50 : if(n!=std::string::npos && n<libName.length()-1) {
1646 100 : extension=libName.substr(n+1);
1647 : }
1648 50 : if(n!=std::string::npos && n<libName.length()) {
1649 100 : base=libName.substr(0,n);
1650 : }
1651 :
1652 50 : if(extension=="cpp") {
1653 120 : libName="./"+base+"."+config::getVersionLong()+"."+config::getSoExt();
1654 : // full path command, including environment setup
1655 : // this will work even if plumed is not in the execution path or if it has been
1656 : // installed with a name different from "plumed"
1657 150 : std::string cmd=config::getEnvCommand()+" \""+config::getPlumedRoot()+"\"/scripts/mklib.sh -n -o "+libName+" "+fileName;
1658 :
1659 30 : if(std::getenv("PLUMED_LOAD_ACTION_DEBUG")) {
1660 0 : log<<"Executing: "<<cmd;
1661 : } else {
1662 30 : log<<"Compiling: "<<fileName<<" to "<<libName;
1663 : }
1664 :
1665 30 : if(comm.Get_size()>0) {
1666 30 : log<<" (only on master node)";
1667 : }
1668 30 : log<<"\n";
1669 :
1670 : // On MPI process (intracomm), we use Get_rank to make sure a single process does the compilation
1671 : // Processes from multiple replicas might simultaneously do the compilation.
1672 30 : if(comm.Get_rank()==0) {
1673 28 : static Tools::CriticalSectionWithKey<std::string> section;
1674 : // This is only locking commands that are running with identical arguments.
1675 : // It is not necessary for correctness (a second mklib would just result in a no op since
1676 : // the library is already there, even if running simultaneously).
1677 : // It however decreases the system load if many threads are used.
1678 : auto s=section.startStop(cmd);
1679 28 : int ret=std::system(cmd.c_str());
1680 28 : if(ret!=0) {
1681 2 : plumed_error() <<"An error happened while executing command "<<cmd<<"\n";
1682 : }
1683 28 : }
1684 29 : comm.Barrier();
1685 : } else {
1686 40 : libName=base+"."+config::getSoExt();
1687 : }
1688 :
1689 : // If we have multiple threads (each holding a Plumed object), each of them
1690 : // will load the library, but each of them will only see actions registered
1691 : // from the owned library
1692 49 : auto *p=dlloader.load(libName);
1693 50 : log<<"Loading shared library "<<libName.c_str()<<" at "<<p<<"\n";
1694 49 : log<<"Here is the list of new actions\n";
1695 49 : log<<"\n";
1696 98 : for(const auto & a : actionRegister().getKeysWithDLHandle(p)) {
1697 49 : log<<a<<"\n";
1698 49 : }
1699 49 : log<<"\n";
1700 : } else {
1701 0 : plumed_error()<<"While loading library "<< fileName << " loading was not enabled, please check if dlopen was found at configure time";
1702 : }
1703 49 : }
1704 :
1705 393572 : void PlumedMain::resetInputs() {
1706 1250418 : for(const auto & ip : inputs) {
1707 856846 : if( ip->isActive() && ip->hasBeenSet() ) {
1708 380649 : ip->reset();
1709 : }
1710 : }
1711 393572 : }
1712 :
1713 78716 : double PlumedMain::getBias() const {
1714 78716 : return bias;
1715 : }
1716 :
1717 474 : double PlumedMain::getWork() const {
1718 474 : return work;
1719 : }
1720 :
1721 142 : FILE* PlumedMain::fopen(const char *path, const char *mode) {
1722 142 : std::string mmode(mode);
1723 143 : std::string ppath(path);
1724 142 : std::string ppathsuf=ppath+getSuffix();
1725 142 : FILE*fp=std::fopen(const_cast<char*>(ppathsuf.c_str()),const_cast<char*>(mmode.c_str()));
1726 142 : if(!fp) {
1727 1 : fp=std::fopen(const_cast<char*>(ppath.c_str()),const_cast<char*>(mmode.c_str()));
1728 : }
1729 1 : if(!fp) {
1730 2 : plumed_error() << "file " << ppath << " cannot be found";
1731 : }
1732 141 : return fp;
1733 : }
1734 :
1735 187 : int PlumedMain::fclose(FILE*fp) {
1736 187 : return std::fclose(fp);
1737 : }
1738 :
1739 8078 : std::string PlumedMain::cite(const std::string&item) {
1740 8078 : return citations.cite(item);
1741 : }
1742 :
1743 4156 : void PlumedMain::fflush() {
1744 9081 : for(const auto & p : files) {
1745 4925 : p->flush();
1746 : }
1747 4156 : }
1748 :
1749 6765 : void PlumedMain::insertFile(FileBase&f) {
1750 6765 : files.insert(&f);
1751 6765 : }
1752 :
1753 7312 : void PlumedMain::eraseFile(FileBase&f) {
1754 7312 : files.erase(&f);
1755 7312 : }
1756 :
1757 66 : void PlumedMain::stop() {
1758 66 : stopNow=true;
1759 66 : }
1760 :
1761 1270 : void PlumedMain::runJobsAtEndOfCalculation() {
1762 52079 : for(const auto & p : actionSet) {
1763 50809 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(p.get());
1764 50809 : if( av && av->calculateOnUpdate() ) {
1765 427 : p->activate();
1766 : }
1767 : }
1768 52078 : for(const auto & p : actionSet) {
1769 50809 : ActionPilot* ap=dynamic_cast<ActionPilot*>(p.get());
1770 50809 : ActionWithValue* av=dynamic_cast<ActionWithValue*>(p.get());
1771 50809 : if( av && av->calculateOnUpdate() ) {
1772 427 : p->calculate();
1773 50382 : } else if( ap && !av && ap->getStride()==0 ) {
1774 87 : p->update();
1775 : } else {
1776 50295 : p->runFinalJobs();
1777 : }
1778 : }
1779 1269 : }
1780 :
1781 8810295 : unsigned PlumedMain::increaseReferenceCounter() noexcept {
1782 8810295 : return ++referenceCounter;
1783 : }
1784 :
1785 8808797 : unsigned PlumedMain::decreaseReferenceCounter() noexcept {
1786 8808797 : return --referenceCounter;
1787 : }
1788 :
1789 41 : unsigned PlumedMain::useCountReferenceCounter() const noexcept {
1790 41 : return referenceCounter;
1791 : }
1792 :
1793 60 : bool PlumedMain::valueExists( const std::string& name ) const {
1794 240 : for(const auto & p : inputs) {
1795 240 : if( p->getLabel()==name ) {
1796 : return true;
1797 : }
1798 : }
1799 : return false;
1800 : }
1801 :
1802 973772 : void PlumedMain::setInputValue( const std::string& name, const unsigned& start, const unsigned& stride, const TypesafePtr & val ) {
1803 : bool found=false;
1804 2304153 : for(const auto & pp : inputs) {
1805 2304153 : if( pp->setValuePointer( name, val ) ) {
1806 973768 : pp->setStart(name, start);
1807 973768 : pp->setStride(name, stride);
1808 : found=true;
1809 : break;
1810 : }
1811 : }
1812 0 : plumed_massert( found, "found no action to set named " + name );
1813 973768 : }
1814 :
1815 657478 : void PlumedMain::setInputForce( const std::string& name, const TypesafePtr & val ) {
1816 : bool found=false;
1817 1512080 : for(const auto & pp : inputs) {
1818 1512080 : if( pp->setForcePointer( name, val ) ) {
1819 : found=true;
1820 : break;
1821 : }
1822 : }
1823 657460 : plumed_massert( found, "found no action to set named " + name );
1824 657460 : }
1825 :
1826 2269 : void PlumedMain::setUnits( const bool& natural, const Units& u ) {
1827 2269 : passtools->usingNaturalUnits = natural;
1828 2269 : passtools->units=u;
1829 2269 : std::vector<ActionToPutData*> idata = actionSet.select<ActionToPutData*>();
1830 17574 : for(const auto & ip : idata) {
1831 15305 : ip->updateUnits( passtools.get() );
1832 : }
1833 61209 : for(const auto & p : actionSet ) {
1834 58940 : p->resetStoredTimestep();
1835 : }
1836 2269 : }
1837 :
1838 393515 : void PlumedMain::startStep() {
1839 1250134 : for(const auto & ip : inputs) {
1840 856619 : ip->resetForStepStart();
1841 : }
1842 393515 : }
1843 :
1844 126 : void PlumedMain::writeBinary(std::ostream&o)const {
1845 504 : for(const auto & ip : inputs) {
1846 378 : ip->writeBinary(o);
1847 : }
1848 126 : }
1849 :
1850 126 : void PlumedMain::readBinary(std::istream&i) {
1851 504 : for(const auto & ip : inputs) {
1852 378 : ip->readBinary(i);
1853 : }
1854 126 : }
1855 :
1856 41 : void PlumedMain::setEnergyValue( const std::string& name ) {
1857 41 : name_of_energy = name;
1858 41 : }
1859 :
1860 15586 : int PlumedMain::getRealPrecision() const {
1861 15586 : return passtools->getRealPrecision();
1862 : }
1863 :
1864 6310 : bool PlumedMain::usingNaturalUnits() const {
1865 6310 : return passtools->usingNaturalUnits;
1866 : }
1867 :
1868 15308954 : const Units& PlumedMain::getUnits() {
1869 15308954 : return passtools->units;
1870 : }
1871 :
1872 4 : PlumedMain::DeprecatedAtoms& PlumedMain::getAtoms() {
1873 4 : return datoms;
1874 : }
1875 :
1876 78228 : void PlumedMain::plumedQuantityToMD( const std::string& unit, const double& eng, const TypesafePtr & m) const {
1877 78228 : passtools->double2MD( eng/passtools->getUnitConversion(unit),m );
1878 78228 : }
1879 :
1880 0 : double PlumedMain::MDQuantityToPLUMED( const std::string& unit, const TypesafePtr & m) const {
1881 0 : double x=passtools->MD2double(m);
1882 0 : return x*passtools->getUnitConversion(unit);
1883 : }
1884 :
1885 1 : double PlumedMain::DeprecatedAtoms::getKBoltzmann() const {
1886 1 : if( plumed.usingNaturalUnits() ) {
1887 : return 1.0;
1888 : }
1889 1 : return kBoltzmann/plumed.getUnits().getEnergy();
1890 : }
1891 :
1892 1 : double PlumedMain::DeprecatedAtoms::getKbT() const {
1893 1 : ActionForInterface* kb=plumed.getActionSet().selectWithLabel<ActionForInterface*>("kBT");
1894 1 : if( kb ) {
1895 1 : return (kb->copyOutput(0))->get();
1896 : }
1897 : return 0.0;
1898 : }
1899 :
1900 1 : int PlumedMain::DeprecatedAtoms::getNatoms() const {
1901 1 : std::vector<ActionToPutData*> atpd=plumed.getActionSet().select<ActionToPutData*>();
1902 1 : for(const auto & pp : atpd ) {
1903 1 : if( pp->getRole()=="x" ) {
1904 1 : return (pp->copyOutput(0))->getShape()[0];
1905 : }
1906 : }
1907 : return 0;
1908 : }
1909 :
1910 1 : bool PlumedMain::DeprecatedAtoms::usingNaturalUnits() const {
1911 1 : return plumed.usingNaturalUnits();
1912 : }
1913 :
1914 0 : void PlumedMain::DeprecatedAtoms::setCollectEnergy(bool b) const {
1915 0 : plumed.readInputWords( Tools::getWords(plumed.MDEngine + "_energy: ENERGY"), false );
1916 0 : plumed.setEnergyValue( plumed.MDEngine + "_energy" );
1917 0 : }
1918 :
1919 0 : double PlumedMain::DeprecatedAtoms::getEnergy() const {
1920 0 : ActionToPutData* av = plumed.getActionSet().selectWithLabel<ActionToPutData*>( plumed.MDEngine + "_energy" );
1921 0 : return (av->copyOutput(0))->get();
1922 : }
1923 :
1924 17 : void PlumedMain::activateParseOnlyMode() {
1925 17 : doParseOnly=true;
1926 17 : }
1927 :
1928 3179 : bool PlumedMain::parseOnlyMode() const {
1929 3179 : return doParseOnly;
1930 : }
1931 :
1932 96 : void PlumedMain::getKeywordsForAction( const std::string& action, Keywords& keys ) const {
1933 96 : actionRegister().getKeywords( dlloader.getHandles(), action, keys );
1934 96 : }
1935 :
1936 5004 : bool PlumedMain::checkAction(const std::string& action) const {
1937 5004 : return actionRegister().check(dlloader.getHandles(),action);
1938 : }
1939 : #ifdef __PLUMED_HAS_PYTHON
1940 : // This is here to stop cppcheck throwing an error
1941 : #endif
1942 :
1943 : #ifdef __PLUMED_HAS_DLADDR
1944 : // This is here to stop cppcheck throwing an error
1945 : #endif
1946 :
1947 : }
1948 :
1949 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
