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 "ActionRegister.h"
26 : #include "ActionSet.h"
27 : #include "ActionWithValue.h"
28 : #include "ActionWithVirtualAtom.h"
29 : #include "Atoms.h"
30 : #include "CLToolMain.h"
31 : #include "ExchangePatterns.h"
32 : #include "GREX.h"
33 : #include "config/Config.h"
34 : #include "tools/Citations.h"
35 : #include "tools/Communicator.h"
36 : #include "tools/DLLoader.h"
37 : #include "tools/Exception.h"
38 : #include "tools/IFile.h"
39 : #include "tools/Log.h"
40 : #include "tools/OpenMP.h"
41 : #include "tools/Tools.h"
42 : #include "tools/Stopwatch.h"
43 : #include "tools/TypesafePtr.h"
44 : #include "lepton/Exception.h"
45 : #include "DataFetchingObject.h"
46 : #include <cstdlib>
47 : #include <cstdio>
48 : #include <cstring>
49 : #include <set>
50 : #include <unordered_map>
51 : #include <exception>
52 : #include <stdexcept>
53 : #include <ios>
54 : #include <new>
55 : #include <typeinfo>
56 : #include <iostream>
57 : #include <algorithm>
58 : #ifdef __PLUMED_LIBCXX11
59 : #include <system_error>
60 : #include <future>
61 : #include <memory>
62 : #include <functional>
63 : #endif
64 :
65 :
66 : namespace PLMD {
67 :
68 : /// Small utility just used in this file to throw arbitrary exceptions
69 26 : [[noreturn]] static void testThrow(const char* what) {
70 52 : auto words=Tools::getWords(what);
71 26 : plumed_assert(words.size()>0);
72 : #define __PLUMED_THROW_NOMSG(type) if(words[0]==#type) throw type()
73 : #define __PLUMED_THROW_MSG(type) if(words[0]==#type) throw type(what)
74 27 : __PLUMED_THROW_MSG(PLMD::ExceptionError);
75 26 : __PLUMED_THROW_MSG(PLMD::ExceptionDebug);
76 25 : __PLUMED_THROW_MSG(PLMD::Exception);
77 24 : __PLUMED_THROW_MSG(PLMD::lepton::Exception);
78 22 : __PLUMED_THROW_NOMSG(std::bad_exception);
79 : #ifdef __PLUMED_LIBCXX11
80 21 : __PLUMED_THROW_NOMSG(std::bad_array_new_length);
81 : #endif
82 20 : __PLUMED_THROW_NOMSG(std::bad_alloc);
83 : #ifdef __PLUMED_LIBCXX11
84 19 : __PLUMED_THROW_NOMSG(std::bad_function_call);
85 18 : __PLUMED_THROW_NOMSG(std::bad_weak_ptr);
86 : #endif
87 17 : __PLUMED_THROW_NOMSG(std::bad_cast);
88 16 : __PLUMED_THROW_NOMSG(std::bad_typeid);
89 15 : __PLUMED_THROW_MSG(std::underflow_error);
90 14 : __PLUMED_THROW_MSG(std::overflow_error);
91 13 : __PLUMED_THROW_MSG(std::range_error);
92 12 : __PLUMED_THROW_MSG(std::runtime_error);
93 11 : __PLUMED_THROW_MSG(std::out_of_range);
94 10 : __PLUMED_THROW_MSG(std::length_error);
95 9 : __PLUMED_THROW_MSG(std::domain_error);
96 8 : __PLUMED_THROW_MSG(std::invalid_argument);
97 7 : __PLUMED_THROW_MSG(std::logic_error);
98 :
99 : #ifdef __PLUMED_LIBCXX11
100 6 : if(words[0]=="std::system_error") {
101 4 : plumed_assert(words.size()>2);
102 : int error_code;
103 4 : Tools::convert(words[2],error_code);
104 4 : if(words[1]=="std::generic_category") throw std::system_error(error_code,std::generic_category(),what);
105 3 : if(words[1]=="std::system_category") throw std::system_error(error_code,std::system_category(),what);
106 2 : if(words[1]=="std::iostream_category") throw std::system_error(error_code,std::iostream_category(),what);
107 1 : if(words[1]=="std::future_category") throw std::system_error(error_code,std::future_category(),what);
108 : }
109 : #endif
110 :
111 2 : if(words[0]=="std::ios_base::failure") {
112 : #ifdef __PLUMED_LIBCXX11
113 1 : int error_code=0;
114 1 : if(words.size()>2) Tools::convert(words[2],error_code);
115 1 : if(words.size()>1 && words[1]=="std::generic_category") throw std::ios_base::failure(what,std::error_code(error_code,std::generic_category()));
116 1 : if(words.size()>1 && words[1]=="std::system_category") throw std::ios_base::failure(what,std::error_code(error_code,std::system_category()));
117 1 : if(words.size()>1 && words[1]=="std::iostream_category") throw std::ios_base::failure(what,std::error_code(error_code,std::iostream_category()));
118 1 : if(words.size()>1 && words[1]=="std::future_category") throw std::ios_base::failure(what,std::error_code(error_code,std::future_category()));
119 : #endif
120 2 : throw std::ios_base::failure(what);
121 : }
122 :
123 2 : plumed_error() << "unknown exception " << what;
124 26 : }
125 :
126 : #if 0
127 : // this is temporarily commented out since it makes
128 : // multithread calculations fail
129 : namespace {
130 : class Register {
131 : std::vector<PlumedMain*> instances;
132 : public:
133 : void add(PlumedMain* instance) {
134 : instances.push_back(instance);
135 : }
136 : void remove(PlumedMain* instance) {
137 : auto it = std::find(instances.begin(), instances.end(), instance);
138 : if(it==instances.end()) {
139 : std::cerr<<"WARNING: internal inconsistency in allocated PlumedMain instances\n";
140 : } else {
141 : instances.erase(it);
142 : }
143 : }
144 : ~Register() {
145 : if(instances.size()>0) std::cerr<<"PLUMED instances was not properly deallocated in your code: "<<instances.size()<<"\n";
146 : }
147 : };
148 :
149 : static Register myregister;
150 :
151 : }
152 : #endif
153 :
154 404395 : PlumedMain::PlumedMain():
155 404395 : initialized(false),
156 : // automatically write on log in destructor
157 404395 : stopwatch_fwd(log),
158 404395 : step(0),
159 404395 : active(false),
160 404395 : mydatafetcher(DataFetchingObject::create(sizeof(double),*this)),
161 404395 : endPlumed(false),
162 404395 : atoms_fwd(*this),
163 404395 : actionSet_fwd(*this),
164 404395 : bias(0.0),
165 404395 : work(0.0),
166 404395 : exchangeStep(false),
167 404395 : restart(false),
168 404395 : doCheckPoint(false),
169 404395 : stopNow(false),
170 404395 : novirial(false),
171 808790 : detailedTimers(false)
172 : {
173 404395 : log.link(comm);
174 404395 : log.setLinePrefix("PLUMED: ");
175 : //myregister.add(this);
176 404395 : }
177 :
178 : // destructor needed to delete forward declarated objects
179 808069 : PlumedMain::~PlumedMain() {
180 : //myregister.remove(this);
181 1212464 : }
182 :
183 : /////////////////////////////////////////////////////////////
184 : // MAIN INTERPRETER
185 :
186 : #define CHECK_INIT(ini,word) plumed_massert(ini,"cmd(\"" + word +"\") should be only used after plumed initialization")
187 : #define CHECK_NOTINIT(ini,word) plumed_massert(!(ini),"cmd(\"" + word +"\") should be only used before plumed initialization")
188 : #define CHECK_NOTNULL(val,word) plumed_massert(val,"NULL pointer received in cmd(\"" + word + "\")");
189 :
190 :
191 1040740 : void PlumedMain::cmd(const std::string & word,const TypesafePtr & val) {
192 :
193 : // Enumerate all possible commands:
194 : enum {
195 : #include "PlumedMainEnum.inc"
196 : };
197 :
198 : // Static object (initialized once) containing the map of commands:
199 : const static std::unordered_map<std::string, int> word_map = {
200 : #include "PlumedMainMap.inc"
201 1306698 : };
202 :
203 : try {
204 :
205 1040740 : auto ss=stopwatch.startPause();
206 :
207 1040740 : std::vector<std::string> words=Tools::getWords(word);
208 1040740 : unsigned nw=words.size();
209 1040740 : if(nw==0) {
210 : // do nothing
211 : } else {
212 : int iword=-1;
213 : double d;
214 : const auto it=word_map.find(words[0]);
215 1040740 : if(it!=word_map.end()) iword=it->second;
216 1040733 : switch(iword) {
217 43333 : case cmd_setBox:
218 43333 : CHECK_INIT(initialized,word);
219 43333 : CHECK_NOTNULL(val,word);
220 43333 : atoms.setBox(val);
221 : break;
222 47474 : case cmd_setPositions:
223 47474 : CHECK_INIT(initialized,word);
224 47474 : atoms.setPositions(val);
225 : break;
226 47478 : case cmd_setMasses:
227 47478 : CHECK_INIT(initialized,word);
228 47478 : atoms.setMasses(val);
229 : break;
230 37570 : case cmd_setCharges:
231 37570 : CHECK_INIT(initialized,word);
232 37570 : atoms.setCharges(val);
233 : break;
234 1 : case cmd_setPositionsX:
235 1 : CHECK_INIT(initialized,word);
236 1 : atoms.setPositions(val,0);
237 : break;
238 1 : case cmd_setPositionsY:
239 1 : CHECK_INIT(initialized,word);
240 1 : atoms.setPositions(val,1);
241 : break;
242 1 : case cmd_setPositionsZ:
243 1 : CHECK_INIT(initialized,word);
244 1 : atoms.setPositions(val,2);
245 : break;
246 37675 : case cmd_setVirial:
247 37675 : CHECK_INIT(initialized,word);
248 37675 : CHECK_NOTNULL(val,word);
249 37675 : atoms.setVirial(val);
250 : break;
251 9792 : case cmd_setEnergy:
252 9792 : CHECK_INIT(initialized,word);
253 9792 : CHECK_NOTNULL(val,word);
254 9792 : atoms.setEnergy(val);
255 : break;
256 47474 : case cmd_setForces:
257 47474 : CHECK_INIT(initialized,word);
258 47474 : atoms.setForces(val);
259 : break;
260 1 : case cmd_setForcesX:
261 1 : CHECK_INIT(initialized,word);
262 1 : atoms.setForces(val,0);
263 : break;
264 1 : case cmd_setForcesY:
265 1 : CHECK_INIT(initialized,word);
266 1 : atoms.setForces(val,1);
267 : break;
268 1 : case cmd_setForcesZ:
269 1 : CHECK_INIT(initialized,word);
270 1 : atoms.setForces(val,2);
271 : break;
272 245716 : case cmd_calc:
273 245716 : CHECK_INIT(initialized,word);
274 245716 : calc();
275 : break;
276 99 : case cmd_prepareDependencies:
277 99 : CHECK_INIT(initialized,word);
278 99 : prepareDependencies();
279 : break;
280 84 : case cmd_shareData:
281 84 : CHECK_INIT(initialized,word);
282 84 : shareData();
283 : break;
284 2070 : case cmd_prepareCalc:
285 2070 : CHECK_INIT(initialized,word);
286 2070 : prepareCalc();
287 : break;
288 10 : case cmd_performCalc:
289 10 : CHECK_INIT(initialized,word);
290 10 : performCalc();
291 : break;
292 2134 : case cmd_performCalcNoUpdate:
293 2134 : CHECK_INIT(initialized,word);
294 2134 : performCalcNoUpdate();
295 : break;
296 10 : case cmd_performCalcNoForces:
297 10 : CHECK_INIT(initialized,word);
298 10 : performCalcNoForces();
299 : break;
300 79 : case cmd_update:
301 79 : CHECK_INIT(initialized,word);
302 79 : update();
303 : break;
304 10028 : case cmd_setStep:
305 10028 : CHECK_INIT(initialized,word);
306 10034 : CHECK_NOTNULL(val,word);
307 10025 : step=val.get<int>();
308 10025 : atoms.startStep();
309 : break;
310 237805 : case cmd_setStepLong:
311 237805 : CHECK_INIT(initialized,word);
312 237805 : CHECK_NOTNULL(val,word);
313 237805 : step=val.get<long int>();
314 237805 : atoms.startStep();
315 : break;
316 : // words used less frequently:
317 1140 : case cmd_setAtomsNlocal:
318 1140 : CHECK_INIT(initialized,word);
319 1140 : CHECK_NOTNULL(val,word);
320 1140 : atoms.setAtomsNlocal(val.get<int>());
321 : break;
322 988 : case cmd_setAtomsGatindex:
323 988 : CHECK_INIT(initialized,word);
324 988 : atoms.setAtomsGatindex(val,false);
325 : break;
326 0 : case cmd_setAtomsFGatindex:
327 0 : CHECK_INIT(initialized,word);
328 0 : atoms.setAtomsGatindex(val,true);
329 : break;
330 152 : case cmd_setAtomsContiguous:
331 152 : CHECK_INIT(initialized,word);
332 152 : CHECK_NOTNULL(val,word);
333 152 : atoms.setAtomsContiguous(val.get<int>());
334 : break;
335 116 : case cmd_createFullList:
336 116 : CHECK_INIT(initialized,word);
337 116 : CHECK_NOTNULL(val,word);
338 116 : atoms.createFullList(val);
339 : break;
340 116 : case cmd_getFullList:
341 116 : CHECK_INIT(initialized,word);
342 116 : CHECK_NOTNULL(val,word);
343 116 : atoms.getFullList(val);
344 : break;
345 116 : case cmd_clearFullList:
346 116 : CHECK_INIT(initialized,word);
347 116 : atoms.clearFullList();
348 : break;
349 : /* ADDED WITH API==6 */
350 64 : case cmd_getDataRank:
351 64 : CHECK_INIT(initialized,words[0]); plumed_assert(nw==2 || nw==3);
352 128 : if( nw==2 ) DataFetchingObject::get_rank( actionSet, words[1], "", val);
353 0 : else DataFetchingObject::get_rank( actionSet, words[1], words[2], val);
354 : break;
355 : /* ADDED WITH API==6 */
356 0 : case cmd_getDataShape:
357 0 : CHECK_INIT(initialized,words[0]); plumed_assert(nw==2 || nw==3);
358 0 : if( nw==2 ) DataFetchingObject::get_shape( actionSet, words[1], "", val );
359 0 : else DataFetchingObject::get_shape( actionSet, words[1], words[2], val );
360 : break;
361 : /* ADDED WITH API==6 */
362 64 : case cmd_setMemoryForData:
363 64 : CHECK_INIT(initialized,words[0]); plumed_assert(nw==2 || nw==3);
364 128 : if( nw==2 ) mydatafetcher->setData( words[1], "", val );
365 0 : else mydatafetcher->setData( words[1], words[2], val );
366 : break;
367 : /* ADDED WITH API==6 */
368 : case cmd_setErrorHandler:
369 : {
370 0 : if(val) error_handler=*static_cast<const plumed_error_handler*>(val.get<const void*>());
371 0 : else error_handler.handler=NULL;
372 : }
373 : break;
374 0 : case cmd_read:
375 0 : CHECK_INIT(initialized,word);
376 0 : if(val)readInputFile(val.get<const char*>());
377 0 : else readInputFile("plumed.dat");
378 : break;
379 269 : case cmd_readInputLine:
380 269 : CHECK_INIT(initialized,word);
381 269 : CHECK_NOTNULL(val,word);
382 269 : readInputLine(val.get<const char*>());
383 230 : break;
384 1 : case cmd_readInputLines:
385 1 : CHECK_INIT(initialized,word);
386 1 : CHECK_NOTNULL(val,word);
387 1 : readInputLines(val.get<const char*>());
388 1 : break;
389 1 : case cmd_clear:
390 1 : CHECK_INIT(initialized,word);
391 1 : actionSet.clearDelete();
392 : break;
393 : case cmd_getApiVersion:
394 21 : CHECK_NOTNULL(val,word);
395 21 : val.set(int(9));
396 : break;
397 : // commands which can be used only before initialization:
398 914 : case cmd_init:
399 914 : CHECK_NOTINIT(initialized,word);
400 914 : init();
401 : break;
402 781 : case cmd_setRealPrecision:
403 781 : CHECK_NOTINIT(initialized,word);
404 781 : CHECK_NOTNULL(val,word);
405 781 : atoms.setRealPrecision(val.get<int>());
406 1560 : mydatafetcher=DataFetchingObject::create(val.get<int>(),*this);
407 780 : break;
408 718 : case cmd_setMDLengthUnits:
409 718 : CHECK_NOTINIT(initialized,word);
410 718 : CHECK_NOTNULL(val,word);
411 718 : atoms.MD2double(val,d);
412 718 : atoms.setMDLengthUnits(d);
413 : break;
414 718 : case cmd_setMDChargeUnits:
415 718 : CHECK_NOTINIT(initialized,word);
416 718 : CHECK_NOTNULL(val,word);
417 718 : atoms.MD2double(val,d);
418 718 : atoms.setMDChargeUnits(d);
419 : break;
420 718 : case cmd_setMDMassUnits:
421 718 : CHECK_NOTINIT(initialized,word);
422 718 : CHECK_NOTNULL(val,word);
423 718 : atoms.MD2double(val,d);
424 718 : atoms.setMDMassUnits(d);
425 : break;
426 45 : case cmd_setMDEnergyUnits:
427 45 : CHECK_NOTINIT(initialized,word);
428 45 : CHECK_NOTNULL(val,word);
429 45 : atoms.MD2double(val,d);
430 45 : atoms.setMDEnergyUnits(d);
431 : break;
432 6 : case cmd_setMDTimeUnits:
433 6 : CHECK_NOTINIT(initialized,word);
434 6 : CHECK_NOTNULL(val,word);
435 6 : atoms.MD2double(val,d);
436 6 : atoms.setMDTimeUnits(d);
437 : break;
438 0 : case cmd_setNaturalUnits:
439 : // set the boltzman constant for MD in natural units (kb=1)
440 : // only needed in LJ codes if the MD is passing temperatures to plumed (so, not yet...)
441 : // use as cmd("setNaturalUnits")
442 0 : CHECK_NOTINIT(initialized,word);
443 0 : atoms.setMDNaturalUnits(true);
444 : break;
445 52 : case cmd_setNoVirial:
446 52 : CHECK_NOTINIT(initialized,word);
447 52 : novirial=true;
448 52 : break;
449 773 : case cmd_setPlumedDat:
450 773 : CHECK_NOTINIT(initialized,word);
451 773 : CHECK_NOTNULL(val,word);
452 773 : plumedDat=val.get<const char*>();
453 : break;
454 337 : case cmd_setMPIComm:
455 337 : CHECK_NOTINIT(initialized,word);
456 337 : comm.Set_comm(val);
457 337 : atoms.setDomainDecomposition(comm);
458 : break;
459 0 : case cmd_setMPIFComm:
460 0 : CHECK_NOTINIT(initialized,word);
461 0 : comm.Set_fcomm(val);
462 0 : atoms.setDomainDecomposition(comm);
463 : break;
464 0 : case cmd_setMPImultiSimComm:
465 0 : CHECK_NOTINIT(initialized,word);
466 0 : multi_sim_comm.Set_comm(val);
467 : break;
468 897 : case cmd_setNatoms:
469 897 : CHECK_NOTINIT(initialized,word);
470 897 : CHECK_NOTNULL(val,word);
471 897 : atoms.setNatoms(val.get<int>());
472 : break;
473 774 : case cmd_setTimestep:
474 774 : CHECK_NOTINIT(initialized,word);
475 774 : CHECK_NOTNULL(val,word);
476 774 : atoms.setTimeStep(val);
477 : break;
478 : /* ADDED WITH API==2 */
479 59 : case cmd_setKbT:
480 59 : CHECK_NOTINIT(initialized,word);
481 59 : CHECK_NOTNULL(val,word);
482 59 : atoms.setKbT(val);
483 : break;
484 : /* ADDED WITH API==3 */
485 8 : case cmd_setRestart:
486 8 : CHECK_NOTINIT(initialized,word);
487 8 : CHECK_NOTNULL(val,word);
488 8 : if(val.get<int>()!=0) restart=true;
489 : break;
490 : /* ADDED WITH API==4 */
491 0 : case cmd_doCheckPoint:
492 0 : CHECK_INIT(initialized,word);
493 0 : CHECK_NOTNULL(val,word);
494 0 : doCheckPoint = false;
495 0 : if(val.get<int>()!=0) doCheckPoint = true;
496 : break;
497 : /* ADDED WITH API==6 */
498 : case cmd_setNumOMPthreads:
499 0 : CHECK_NOTNULL(val,word);
500 : {
501 0 : auto nt=val.get<unsigned>();
502 : if(nt==0) nt=1;
503 0 : OpenMP::setNumThreads(nt);
504 : }
505 : break;
506 : /* ADDED WITH API==6 */
507 : /* only used for testing */
508 : case cmd_throw:
509 26 : CHECK_NOTNULL(val,word);
510 26 : testThrow(val.get<const char*>());
511 : /* STOP API */
512 771 : case cmd_setMDEngine:
513 771 : CHECK_NOTINIT(initialized,word);
514 771 : CHECK_NOTNULL(val,word);
515 771 : MDEngine=val.get<const char*>();
516 : break;
517 758 : case cmd_setLog:
518 758 : CHECK_NOTINIT(initialized,word);
519 758 : log.link(val.get<FILE*>());
520 : break;
521 52 : case cmd_setLogFile:
522 52 : CHECK_NOTINIT(initialized,word);
523 52 : CHECK_NOTNULL(val,word);
524 138 : log.open(val.get<const char*>());
525 52 : break;
526 : // other commands that should be used after initialization:
527 245643 : case cmd_setStopFlag:
528 245643 : CHECK_INIT(initialized,word);
529 245643 : CHECK_NOTNULL(val,word);
530 245643 : val.get<int*>(); // just check type and discard pointer
531 245642 : stopFlag=val.copy();
532 245642 : break;
533 0 : case cmd_getExchangesFlag:
534 0 : CHECK_INIT(initialized,word);
535 0 : CHECK_NOTNULL(val,word);
536 0 : exchangePatterns.getFlag(*val.get<int*>()); // note: getFlag changes the value of the reference!
537 : break;
538 0 : case cmd_setExchangesSeed:
539 0 : CHECK_INIT(initialized,word);
540 0 : CHECK_NOTNULL(val,word);
541 0 : exchangePatterns.setSeed(val.get<int>());
542 : break;
543 0 : case cmd_setNumberOfReplicas:
544 0 : CHECK_INIT(initialized,word);
545 0 : CHECK_NOTNULL(val,word);
546 0 : exchangePatterns.setNofR(val.get<int>());
547 : break;
548 0 : case cmd_getExchangesList:
549 0 : CHECK_INIT(initialized,word);
550 0 : CHECK_NOTNULL(val,word);
551 0 : exchangePatterns.getList(val.get<int*>());
552 0 : break;
553 721 : case cmd_runFinalJobs:
554 721 : CHECK_INIT(initialized,word);
555 721 : runJobsAtEndOfCalculation();
556 : break;
557 96 : case cmd_isEnergyNeeded:
558 96 : CHECK_INIT(initialized,word);
559 96 : CHECK_NOTNULL(val,word);
560 96 : if(atoms.isEnergyNeeded()) val.set(int(1));
561 96 : else val.set(int(0));
562 : break;
563 2172 : case cmd_getBias:
564 2172 : CHECK_INIT(initialized,word);
565 2172 : CHECK_NOTNULL(val,word);
566 2172 : atoms.double2MD(getBias()/(atoms.getMDUnits().getEnergy()/atoms.getUnits().getEnergy()),val);
567 2172 : break;
568 : case cmd_checkAction:
569 2 : CHECK_NOTNULL(val,word);
570 2 : plumed_assert(nw==2);
571 3 : val.set(int(actionRegister().check(words[1]) ? 1:0));
572 2 : break;
573 : case cmd_setExtraCV:
574 20 : CHECK_NOTNULL(val,word);
575 20 : plumed_assert(nw==2);
576 20 : atoms.setExtraCV(words[1],val);
577 : break;
578 : case cmd_setExtraCVForce:
579 20 : CHECK_NOTNULL(val,word);
580 20 : plumed_assert(nw==2);
581 20 : atoms.setExtraCVForce(words[1],val);
582 : break;
583 1077 : case cmd_GREX:
584 1077 : if(!grex) grex=Tools::make_unique<GREX>(*this);
585 1077 : plumed_massert(grex,"error allocating grex");
586 : {
587 1077 : std::string kk=words[1];
588 1248 : for(unsigned i=2; i<words.size(); i++) kk+=" "+words[i];
589 2154 : grex->cmd(kk.c_str(),val);
590 : }
591 1077 : break;
592 10668 : case cmd_CLTool:
593 10668 : CHECK_NOTINIT(initialized,word);
594 10668 : if(!cltool) cltool=Tools::make_unique<CLToolMain>();
595 : {
596 10668 : std::string kk=words[1];
597 10668 : for(unsigned i=2; i<words.size(); i++) kk+=" "+words[i];
598 21336 : cltool->cmd(kk.c_str(),val);
599 : }
600 10668 : break;
601 : /* ADDED WITH API==7 */
602 : case cmd_convert:
603 : {
604 : double v;
605 22 : plumed_assert(words.size()==2);
606 22 : if(Tools::convertNoexcept(words[1],v)) atoms.double2MD(v,val);
607 : }
608 22 : break;
609 7 : default:
610 21 : plumed_merror("cannot interpret cmd(\"" + word + "\"). check plumed developers manual to see the available commands.");
611 : break;
612 : }
613 : }
614 :
615 1040912 : } catch (const std::exception &e) {
616 86 : if(log.isOpen()) {
617 51 : log<<"\n\n################################################################################\n\n";
618 51 : log<<e.what();
619 51 : log<<"\n\n################################################################################\n\n";
620 51 : log.flush();
621 : }
622 86 : throw;
623 86 : }
624 1040654 : }
625 :
626 : ////////////////////////////////////////////////////////////////////////
627 :
628 914 : void PlumedMain::init() {
629 : // check that initialization just happens once
630 914 : initialized=true;
631 914 : atoms.init();
632 914 : if(!log.isOpen()) log.link(stdout);
633 914 : log<<"PLUMED is starting\n";
634 2742 : log<<"Version: "<<config::getVersionLong()<<" (git: "<<config::getVersionGit()<<") "
635 3656 : <<"compiled on " <<config::getCompilationDate() << " at " << config::getCompilationTime() << "\n";
636 914 : log<<"Please cite these papers when using PLUMED ";
637 1828 : log<<cite("The PLUMED consortium, Nat. Methods 16, 670 (2019)");
638 1828 : log<<cite("Tribello, Bonomi, Branduardi, Camilloni, and Bussi, Comput. Phys. Commun. 185, 604 (2014)");
639 914 : log<<"\n";
640 914 : log<<"For further information see the PLUMED web page at http://www.plumed.org\n";
641 914 : log<<"Root: "<<config::getPlumedRoot()<<"\n";
642 1828 : log<<"For installed feature, see "<<config::getPlumedRoot() + "/src/config/config.txt\n";
643 914 : log.printf("Molecular dynamics engine: %s\n",MDEngine.c_str());
644 914 : log.printf("Precision of reals: %d\n",atoms.getRealPrecision());
645 1604 : log.printf("Running over %d %s\n",comm.Get_size(),(comm.Get_size()>1?"nodes":"node"));
646 914 : log<<"Number of threads: "<<OpenMP::getNumThreads()<<"\n";
647 914 : log<<"Cache line size: "<<OpenMP::getCachelineSize()<<"\n";
648 914 : log.printf("Number of atoms: %d\n",atoms.getNatoms());
649 914 : if(grex) log.printf("GROMACS-like replica exchange is on\n");
650 914 : log.printf("File suffix: %s\n",getSuffix().c_str());
651 914 : if(plumedDat.length()>0) {
652 773 : readInputFile(plumedDat);
653 : plumedDat="";
654 : }
655 914 : atoms.updateUnits();
656 914 : log.printf("Timestep: %f\n",atoms.getTimeStep());
657 914 : if(atoms.getKbT()>0.0)
658 53 : log.printf("KbT: %f\n",atoms.getKbT());
659 : else {
660 861 : log.printf("KbT has not been set by the MD engine\n");
661 861 : log.printf("It should be set by hand where needed\n");
662 : }
663 914 : log<<"Relevant bibliography:\n";
664 914 : log<<citations;
665 914 : log<<"Please read and cite where appropriate!\n";
666 914 : log<<"Finished setup\n";
667 914 : }
668 :
669 791 : void PlumedMain::readInputFile(const std::string & str) {
670 791 : plumed_assert(initialized);
671 791 : log<<"FILE: "<<str<<"\n";
672 791 : IFile ifile;
673 791 : ifile.link(*this);
674 791 : ifile.open(str);
675 791 : ifile.allowNoEOL();
676 791 : readInputFile(ifile);
677 791 : log<<"END FILE: "<<str<<"\n";
678 791 : log.flush();
679 :
680 791 : }
681 :
682 792 : void PlumedMain::readInputFile(IFile & ifile) {
683 : std::vector<std::string> words;
684 14523 : while(Tools::getParsedLine(ifile,words) && !endPlumed) readInputWords(words);
685 792 : endPlumed=false;
686 792 : pilots=actionSet.select<ActionPilot*>();
687 792 : }
688 :
689 353 : void PlumedMain::readInputLine(const std::string & str) {
690 353 : plumed_assert(initialized);
691 353 : if(str.empty()) return;
692 353 : std::vector<std::string> words=Tools::getWords(str);
693 353 : citations.clear();
694 353 : readInputWords(words);
695 314 : if(!citations.empty()) {
696 2 : log<<"Relevant bibliography:\n";
697 2 : log<<citations;
698 2 : log<<"Please read and cite where appropriate!\n";
699 : }
700 353 : }
701 :
702 1 : void PlumedMain::readInputLines(const std::string & str) {
703 1 : plumed_assert(initialized);
704 1 : if(str.empty()) return;
705 :
706 1 : log<<"FILE: (temporary)\n";
707 :
708 : // Open a temporary file
709 1 : auto fp=std::tmpfile();
710 1 : plumed_assert(fp);
711 :
712 : // make sure file is closed (and thus deleted) also if an exception occurs
713 1 : auto deleter=[](FILE* fp) { std::fclose(fp); };
714 : std::unique_ptr<FILE,decltype(deleter)> fp_deleter(fp,deleter);
715 :
716 1 : auto ret=std::fputs(str.c_str(),fp);
717 1 : plumed_assert(ret!=EOF);
718 :
719 1 : std::rewind(fp);
720 :
721 1 : IFile ifile;
722 1 : ifile.link(*this);
723 1 : ifile.link(fp);
724 1 : ifile.allowNoEOL();
725 :
726 1 : readInputFile(ifile);
727 1 : log<<"END FILE: (temporary)\n";
728 1 : }
729 :
730 14110 : void PlumedMain::readInputWords(const std::vector<std::string> & words) {
731 14110 : plumed_assert(initialized);
732 14110 : if(words.empty())return;
733 14110 : else if(words[0]=="_SET_SUFFIX") {
734 3 : plumed_assert(words.size()==2);
735 : setSuffix(words[1]);
736 : } else {
737 14107 : std::vector<std::string> interpreted(words);
738 14107 : Tools::interpretLabel(interpreted);
739 28177 : auto action=actionRegister().create(ActionOptions(*this,interpreted));
740 14070 : if(!action) {
741 : std::string msg;
742 : msg ="ERROR\nI cannot understand line:";
743 10 : for(unsigned i=0; i<interpreted.size(); ++i) msg+=" "+interpreted[i];
744 : msg+="\nMaybe a missing space or a typo?";
745 2 : log << msg;
746 2 : log.flush();
747 4 : plumed_merror(msg);
748 : }
749 14068 : action->checkRead();
750 14068 : actionSet.emplace_back(std::move(action));
751 14109 : };
752 :
753 28142 : pilots=actionSet.select<ActionPilot*>();
754 : }
755 :
756 : ////////////////////////////////////////////////////////////////////////
757 :
758 0 : void PlumedMain::exit(int c) {
759 0 : comm.Abort(c);
760 0 : }
761 :
762 14105 : Log& PlumedMain::getLog() {
763 14105 : return log;
764 : }
765 :
766 245716 : void PlumedMain::calc() {
767 245716 : prepareCalc();
768 245716 : performCalc();
769 245716 : }
770 :
771 247786 : void PlumedMain::prepareCalc() {
772 247786 : prepareDependencies();
773 247786 : shareData();
774 247786 : }
775 :
776 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
777 : // here we have the main steps in "calc()"
778 : // they can be called individually, but the standard thing is to
779 : // traverse them in this order:
780 248113 : void PlumedMain::prepareDependencies() {
781 :
782 : // Stopwatch is stopped when sw goes out of scope
783 248113 : auto sw=stopwatch.startStop("1 Prepare dependencies");
784 :
785 : // activate all the actions which are on step
786 : // activation is recursive and enables also the dependencies
787 : // before doing that, the prepare() method is called to see if there is some
788 : // new/changed dependency (up to now, only useful for dependences on virtual atoms,
789 : // which can be dynamically changed).
790 :
791 : // First switch off all actions
792 2071737 : for(const auto & p : actionSet) {
793 1823624 : p->deactivate();
794 : }
795 :
796 : // for optimization, an "active" flag remains false if no action at all is active
797 248113 : active=mydatafetcher->activate();
798 1587111 : for(unsigned i=0; i<pilots.size(); ++i) {
799 1338998 : if(pilots[i]->onStep()) {
800 1269461 : pilots[i]->activate();
801 1269461 : active=true;
802 : }
803 : };
804 :
805 : // also, if one of them is the total energy, tell to atoms that energy should be collected
806 2071737 : for(const auto & p : actionSet) {
807 1823624 : if(p->isActive()) {
808 1642537 : if(p->checkNeedsGradients()) p->setOption("GRADIENTS");
809 : }
810 : }
811 :
812 248113 : }
813 :
814 247870 : void PlumedMain::shareData() {
815 : // atom positions are shared (but only if there is something to do)
816 247870 : if(!active)return;
817 : // Stopwatch is stopped when sw goes out of scope
818 246382 : auto sw=stopwatch.startStop("2 Sharing data");
819 246382 : if(atoms.getNatoms()>0) atoms.share();
820 246382 : }
821 :
822 2134 : void PlumedMain::performCalcNoUpdate() {
823 2134 : waitData();
824 2134 : justCalculate();
825 2134 : backwardPropagate();
826 2134 : }
827 :
828 10 : void PlumedMain::performCalcNoForces() {
829 10 : waitData();
830 10 : justCalculate();
831 10 : }
832 :
833 245726 : void PlumedMain::performCalc() {
834 245726 : waitData();
835 245726 : justCalculate();
836 245726 : backwardPropagate();
837 245726 : update();
838 245726 : mydatafetcher->finishDataGrab();
839 245726 : }
840 :
841 247984 : void PlumedMain::waitData() {
842 247984 : if(!active)return;
843 : // Stopwatch is stopped when sw goes out of scope
844 246496 : auto sw=stopwatch.startStop("3 Waiting for data");
845 246496 : if(atoms.getNatoms()>0) atoms.wait();
846 246496 : }
847 :
848 247984 : void PlumedMain::justCalculate() {
849 247984 : if(!active)return;
850 : // Stopwatch is stopped when sw goes out of scope
851 246496 : auto sw=stopwatch.startStop("4 Calculating (forward loop)");
852 246496 : bias=0.0;
853 246496 : work=0.0;
854 :
855 246496 : int iaction=0;
856 : // calculate the active actions in order (assuming *backward* dependence)
857 2060664 : for(const auto & pp : actionSet) {
858 : Action* p(pp.get());
859 1814168 : if(p->isActive()) {
860 : // Stopwatch is stopped when sw goes out of scope.
861 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
862 1641587 : Stopwatch::Handler sw;
863 1641587 : if(detailedTimers) {
864 : std::string actionNumberLabel;
865 1690 : Tools::convert(iaction,actionNumberLabel);
866 1690 : const unsigned m=actionSet.size();
867 5070 : unsigned k=0; unsigned n=1; while(n<m) { n*=10; k++; }
868 1690 : const int pad=k-actionNumberLabel.length();
869 2815 : for(int i=0; i<pad; i++) actionNumberLabel=" "+actionNumberLabel;
870 1690 : sw=stopwatch.startStop("4A "+actionNumberLabel+" "+p->getLabel());
871 : }
872 1641587 : ActionWithValue*av=dynamic_cast<ActionWithValue*>(p);
873 1641587 : ActionAtomistic*aa=dynamic_cast<ActionAtomistic*>(p);
874 : {
875 1641587 : if(av) av->clearInputForces();
876 1641587 : if(av) av->clearDerivatives();
877 : }
878 : {
879 1641587 : if(aa) aa->clearOutputForces();
880 1641587 : if(aa) if(aa->isActive()) aa->retrieveAtoms();
881 : }
882 1641587 : if(p->checkNumericalDerivatives()) p->calculateNumericalDerivatives();
883 1641020 : else p->calculate();
884 : // This retrieves components called bias
885 1641587 : if(av) {
886 1389822 : bias+=av->getOutputQuantity("bias");
887 1389822 : work+=av->getOutputQuantity("work");
888 1389822 : av->setGradientsIfNeeded();
889 : }
890 1641587 : ActionWithVirtualAtom*avv=dynamic_cast<ActionWithVirtualAtom*>(p);
891 1641587 : if(avv)avv->setGradientsIfNeeded();
892 1641587 : }
893 1814168 : iaction++;
894 : }
895 246496 : }
896 :
897 0 : void PlumedMain::justApply() {
898 0 : backwardPropagate();
899 0 : update();
900 0 : }
901 :
902 247860 : void PlumedMain::backwardPropagate() {
903 247860 : if(!active)return;
904 246372 : int iaction=0;
905 : // Stopwatch is stopped when sw goes out of scope
906 246372 : auto sw=stopwatch.startStop("5 Applying (backward loop)");
907 : // apply them in reverse order
908 2059606 : for(auto pp=actionSet.rbegin(); pp!=actionSet.rend(); ++pp) {
909 : const auto & p(pp->get());
910 1813234 : if(p->isActive()) {
911 :
912 : // Stopwatch is stopped when sw goes out of scope.
913 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
914 1640833 : Stopwatch::Handler sw;
915 1640833 : if(detailedTimers) {
916 : std::string actionNumberLabel;
917 1690 : Tools::convert(iaction,actionNumberLabel);
918 1690 : const unsigned m=actionSet.size();
919 5070 : unsigned k=0; unsigned n=1; while(n<m) { n*=10; k++; }
920 1690 : const int pad=k-actionNumberLabel.length();
921 2815 : for(int i=0; i<pad; i++) actionNumberLabel=" "+actionNumberLabel;
922 1690 : sw=stopwatch.startStop("5A "+actionNumberLabel+" "+p->getLabel());
923 : }
924 :
925 1640833 : p->apply();
926 1640833 : ActionAtomistic*a=dynamic_cast<ActionAtomistic*>(p);
927 : // still ActionAtomistic has a special treatment, since they may need to add forces on atoms
928 1640833 : if(a) a->applyForces();
929 :
930 1640833 : }
931 1813234 : iaction++;
932 : }
933 :
934 : // Stopwatch is stopped when sw goes out of scope.
935 : // We explicitly declare a Stopwatch::Handler here to allow for conditional initialization.
936 246372 : Stopwatch::Handler sw1;
937 246485 : if(detailedTimers) sw1=stopwatch.startStop("5B Update forces");
938 : // this is updating the MD copy of the forces
939 246372 : if(atoms.getNatoms()>0) atoms.updateForces();
940 246372 : }
941 :
942 245805 : void PlumedMain::update() {
943 245805 : if(!active)return;
944 :
945 : // Stopwatch is stopped when sw goes out of scope
946 244317 : auto sw=stopwatch.startStop("6 Update");
947 :
948 : // update step (for statistics, etc)
949 244317 : updateFlags.push(true);
950 2051084 : for(const auto & p : actionSet) {
951 1806767 : p->beforeUpdate();
952 3440763 : if(p->isActive() && p->checkUpdate() && updateFlagsTop()) p->update();
953 : }
954 488638 : while(!updateFlags.empty()) updateFlags.pop();
955 : if(!updateFlags.empty()) plumed_merror("non matching changes in the update flags");
956 : // Check that no action has told the calculation to stop
957 244317 : if(stopNow) {
958 65 : if(stopFlag) stopFlag.set(int(1));
959 0 : else plumed_merror("your md code cannot handle plumed stop events - add a call to plumed.comm(stopFlag,stopCondition)");
960 : }
961 :
962 : // flush by default every 10000 steps
963 : // hopefully will not affect performance
964 : // also if receive checkpointing signal
965 244317 : if(step%10000==0||doCheckPoint) {
966 886 : fflush();
967 886 : log.flush();
968 15368 : for(const auto & p : actionSet) p->fflush();
969 : }
970 244317 : }
971 :
972 2 : void PlumedMain::load(const std::string& ss) {
973 2 : if(DLLoader::installed()) {
974 2 : std::string s=ss;
975 2 : size_t n=s.find_last_of(".");
976 3 : std::string extension="";
977 2 : std::string base=s;
978 4 : if(n!=std::string::npos && n<s.length()-1) extension=s.substr(n+1);
979 4 : if(n!=std::string::npos && n<s.length()) base=s.substr(0,n);
980 2 : if(extension=="cpp") {
981 : // full path command, including environment setup
982 : // this will work even if plumed is not in the execution path or if it has been
983 : // installed with a name different from "plumed"
984 4 : std::string cmd=config::getEnvCommand()+" \""+config::getPlumedRoot()+"\"/scripts/mklib.sh "+s;
985 2 : log<<"Executing: "<<cmd;
986 2 : if(comm.Get_size()>0) log<<" (only on master node)";
987 2 : log<<"\n";
988 2 : if(comm.Get_rank()==0) {
989 2 : int ret=std::system(cmd.c_str());
990 3 : if(ret!=0) plumed_error() <<"An error happened while executing command "<<cmd<<"\n";
991 : }
992 1 : comm.Barrier();
993 2 : base="./"+base;
994 : }
995 2 : s=base+"."+config::getSoExt();
996 1 : void *p=dlloader.load(s);
997 1 : if(!p) {
998 0 : plumed_error()<<"I cannot load library " << ss << " " << dlloader.error();
999 : }
1000 2 : log<<"Loading shared library "<<s.c_str()<<"\n";
1001 1 : log<<"Here is the new list of available actions\n";
1002 1 : log<<actionRegister();
1003 0 : } else plumed_error()<<"While loading library "<< ss << " loading was not enabled, please check if dlopen was found at configure time";
1004 1 : }
1005 :
1006 2805 : double PlumedMain::getBias() const {
1007 2805 : return bias;
1008 : }
1009 :
1010 450 : double PlumedMain::getWork() const {
1011 450 : return work;
1012 : }
1013 :
1014 73 : FILE* PlumedMain::fopen(const char *path, const char *mode) {
1015 73 : std::string mmode(mode);
1016 73 : std::string ppath(path);
1017 73 : std::string suffix(getSuffix());
1018 73 : std::string ppathsuf=ppath+suffix;
1019 73 : FILE*fp=std::fopen(const_cast<char*>(ppathsuf.c_str()),const_cast<char*>(mmode.c_str()));
1020 73 : if(!fp) fp=std::fopen(const_cast<char*>(ppath.c_str()),const_cast<char*>(mmode.c_str()));
1021 73 : plumed_massert(fp,"file " + ppath + " cannot be found");
1022 73 : return fp;
1023 : }
1024 :
1025 91 : int PlumedMain::fclose(FILE*fp) {
1026 91 : return std::fclose(fp);
1027 : }
1028 :
1029 3576 : std::string PlumedMain::cite(const std::string&item) {
1030 3576 : return citations.cite(item);
1031 : }
1032 :
1033 1447 : void PlumedMain::fflush() {
1034 4578 : for(const auto & p : files) {
1035 3131 : p->flush();
1036 : }
1037 1447 : }
1038 :
1039 4317 : void PlumedMain::insertFile(FileBase&f) {
1040 4317 : files.insert(&f);
1041 4317 : }
1042 :
1043 4607 : void PlumedMain::eraseFile(FileBase&f) {
1044 4607 : files.erase(&f);
1045 4607 : }
1046 :
1047 65 : void PlumedMain::stop() {
1048 65 : stopNow=true;
1049 65 : }
1050 :
1051 721 : void PlumedMain::runJobsAtEndOfCalculation() {
1052 14277 : for(const auto & p : actionSet) {
1053 13556 : p->runFinalJobs();
1054 : }
1055 721 : }
1056 :
1057 : #ifdef __PLUMED_HAS_PYTHON
1058 : // This is here to stop cppcheck throwing an error
1059 : #endif
1060 :
1061 : }
1062 :
1063 : //////////////////////////////////////////////////////////////////////////////////////////////////////////////////
|
