Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2019-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 "CLTool.h"
23 : #include "core/CLToolRegister.h"
24 : #include "tools/Communicator.h"
25 : #include "tools/Tools.h"
26 : #include "config/Config.h"
27 : #include "tools/PlumedHandle.h"
28 : #include "tools/Stopwatch.h"
29 : #include "tools/Log.h"
30 : #include "tools/DLLoader.h"
31 : #include "tools/Random.h"
32 : #include "tools/TrajectoryParser.h"
33 : #include "tools/AtomDistribution.h"
34 : #include "tools/AtomDistributionFiles.h"
35 :
36 : #include <cstdio>
37 : #include <string>
38 : #include <vector>
39 : #include <fstream>
40 : #include <atomic>
41 : #include <csignal>
42 : #include <cstdio>
43 : #include <random>
44 : #include <algorithm>
45 : #include <chrono>
46 : #include <string_view>
47 : #include <optional>
48 :
49 : namespace PLMD {
50 : namespace cltools {
51 :
52 : //+PLUMEDOC TOOLS benchmark
53 : /*
54 : benchmark is a lightweight reimplementation of [driver](driver.md) that can be used to run benchmark calculations
55 :
56 : The main difference between [driver](driver.md) and benchmark is that benchmark generates a trajectory in memory rather than reading a
57 : trajectory from a file. This approach is better for timing the overhead of the plumed library. If you do similar benchmarking with driver
58 : the timings you get are dominated by the time spent doing the I/O operations that are required to read the trajectory.
59 :
60 : ## Basic usage
61 :
62 : If you want to use benchmark you first create a sample `plumed.dat` file for testing. For example:
63 :
64 : ```plumed
65 : WHOLEMOLECULES ENTITY0=1-10000
66 : p: POSITION ATOM=1
67 : RESTRAINT ARG=p.x KAPPA=1 AT=0
68 : ```
69 :
70 : You can then run this benchmark using the following command:
71 :
72 : ```plumed
73 : plumed benchmark
74 : ```
75 :
76 : Notice, that benchmark will read an input file called `plumed.dat` by default. You can specify a different name for you PLUMED input file
77 : by using the `--plumed` flag.
78 :
79 : ## Running with a different PLUMED version
80 :
81 : If you want to run a benchmark against a previous plumed version in a controlled setting you can do so by using the command:
82 :
83 : ```plumed
84 : plumed-runtime benchmark --kernel /path/to/lib/libplumedKernel.so
85 : ```
86 :
87 : If you use this command the version of PLUMED that is in your environment calls the version of the library that is specified using the
88 : `--kernel` flag. Running the benchmark in this way ensures that you are running in a controlled setting, where systematic errors
89 : in the comparison are minimized.
90 :
91 : !!! warning "using plumed-runtime"
92 :
93 : You use the `plumed-runtime` executable here to avoid conflicts between different
94 : plumed versions. You will find the `plumed-runtime` executable in your path if you are using the non installed version of plumed,
95 : and in `$prefix/lib/plumed` if you installed plumed in $prefix,.
96 :
97 : ## Comparing multiple versions
98 :
99 : The best way to compare two versions of plumed on the same input is to pass multiple colon-separated kernels as shown below:
100 :
101 : ```plumed
102 : plumed-runtime benchmark --kernel /path/to/lib/libplumedKernel.so:/path2/to/lib/libplumedKernel.so:this
103 : ```
104 :
105 : Here `this` means the kernel of the version with which you are running the benchmark. This comparison runs the three
106 : instances simultaneously (alternating them) so that systematic differences in the load of your machine will affect them
107 : to the same extent.
108 :
109 : In case the different versions require modified plumed.dat files, or if you simply want to compare
110 : two different plumed input files that compute the same thing, you can also use multiple plumed input files:
111 :
112 : ```plumed
113 : plumed-runtime benchmark --kernel /path/to/lib/libplumedKernel.so:this --plumed plumed1.dat:plumed2.dat
114 : ```
115 :
116 : Similarly, you might want to run two different inputs using the same kernel by using an input like this:
117 :
118 : ```plumed
119 : plumed-runtime benchmark --plumed plumed1.dat:plumed2.dat
120 : ```
121 :
122 : ## Profiling
123 :
124 : If you want to attach a profiler to the process on the fly, you might find it convenient to use `--nsteps -1`.
125 : This options ensures that the simulation runs forever unless interrupted with CTRL-C. When interrupted, the result of the timers
126 : should be displayed anyway.
127 : You can also set a maximum time for the calculating by using the `--maxtime` flag.
128 :
129 : If you run a profiler when testing multiple PLUMED versions it can be difficult to determine which function is from
130 : each version. We therefore recommended you recompile separate PLUMED instances with a separate C++ namespace (`-DPLMD=PLUMED_version_1`)
131 : so that you will be able to distinguish them. In addition, compiling with `CXXFLAGS="-g -O3"` will make the profiling
132 : report more complete and will likely highlight lines of code that are particularly computationally demanding.
133 :
134 : ## MPI runs
135 :
136 : You can also run a benchmark that emulates a domain decomposition if plumed has been compiled with MPI
137 : and you run with `mpirun` and a command like the one shown below:
138 :
139 : ```plumed
140 : mpirun -np 4 plumed-runtime benchmark
141 : ```
142 :
143 : If you load separate PLUMED instances as discussed above, they should all be compiled against the same MPI version.
144 : Notice that when using MPI signals (CTRL-C) might not work.
145 :
146 : Since some of the data transfer could happen asynchronously, you might want to use the `--sleep` option
147 : to simulate a lag between the `prepareCalc` and `performCalc` actions. This part of the calculation will not contribute
148 : to the output timings, but will obviously slow down your test.
149 :
150 : ## Output
151 :
152 : In the output you will see the usual reports about timings produced by the internal
153 : timers of the tested plumed instances.
154 :
155 : In addition, this tool monitors the timing externally, with some slightly different criterion:
156 :
157 : - First, the initialization (construction of the input) will be shown with a separate timer,
158 : as well as the timing for the first step.
159 : - Second, the timer corresponding to the calculation will be split in three parts, reporting
160 : execution of the first 20% (warm-up) and the next two blocks of 40% each.
161 : - Finally, you might notice some discrepancies because some of the actions that are usually
162 : not expensive are not included in the internal timers. The external timer will
163 : thus provide a better estimate of the total elapsed time that includes everything.
164 :
165 : The internal timers are still useful to monitor what happens at the different stages
166 : of the calculattion. If you want more detailed information you can also use a
167 : [DEBUG](DEBUG.md) action with the `DETAILED_TIMERS`, to determine how much time is spnt in each action.
168 :
169 : When you run multiple version, a comparative analisys of the time spent within PLUMED in the various
170 : instances will be done. For each PLUMED instance you run, this analysis shows the ratio between the total time each PLUMED instance ran for and the total time the first
171 : PLUMED instance ran for. In other words, the first time that the first PLUMED instance ran for is used as the basis for comparisons. Errors on these estimates of the timings
172 : are calculated using bootstrapping and the warm-up phase is discarded in the analysis.
173 :
174 : */
175 : //+ENDPLUMEDOC
176 :
177 : // We use an anonymous namespace here to avoid clashes with variables
178 : // declared in other parts of the code
179 : namespace {
180 :
181 : //this is a sugar for changing idea faster about the rng
182 : using generator = std::mt19937;
183 :
184 : std::atomic<bool> signalReceived(false);
185 :
186 : class SignalHandlerGuard {
187 : public:
188 4 : SignalHandlerGuard(int signal, void (*newHandler)(int)) : signal_(signal) {
189 : // Store the current handler before setting the new one
190 4 : prevHandler_ = std::signal(signal, newHandler);
191 4 : if (prevHandler_ == SIG_ERR) {
192 0 : throw std::runtime_error("Failed to set signal handler");
193 : }
194 4 : }
195 :
196 : ~SignalHandlerGuard() {
197 : // Restore the previous handler upon destruction
198 2 : std::signal(signal_, prevHandler_);
199 2 : }
200 :
201 : // Delete copy constructor and assignment operator to prevent copying
202 : SignalHandlerGuard(const SignalHandlerGuard&) = delete;
203 : SignalHandlerGuard& operator=(const SignalHandlerGuard&) = delete;
204 :
205 : private:
206 : int signal_;
207 : void (*prevHandler_)(int);
208 : };
209 :
210 0 : extern "C" void signalHandler(int signal) {
211 0 : if (signal == SIGINT) {
212 : signalReceived.store(true);
213 0 : fprintf(stderr, "Signal interrupt received\n");
214 : }
215 0 : if (signal == SIGTERM) {
216 : signalReceived.store(true);
217 0 : fprintf(stderr, "Signal termination received\n");
218 : }
219 0 : }
220 :
221 : /// This base class contains members that are movable with default operations
222 : struct KernelBase {
223 : std::string path;
224 : std::string plumed_dat;
225 : PlumedHandle handle;
226 : Stopwatch stopwatch;
227 : std::vector<long long int> timings;
228 : double comparative_timing=-1.0;
229 : double comparative_timing_error=-1.0;
230 4 : KernelBase(const std::string & path_,const std::string & plumed_dat_, Log* log_):
231 4 : path(path_),
232 4 : plumed_dat(plumed_dat_),
233 4 : handle([&]() {
234 4 : if(path_=="this") {
235 4 : return PlumedHandle();
236 : } else {
237 0 : return PlumedHandle::dlopen(path.c_str());
238 : }
239 : }()),
240 4 : stopwatch(*log_) {
241 4 : }
242 : };
243 :
244 : /// Local structure handling a kernel and the related timers.
245 : /// This structure specifically contain the Log, which needs special treatment
246 : /// in move semantics
247 : struct Kernel :
248 : public KernelBase {
249 : Log* log=nullptr;
250 4 : Kernel(const std::string & path_,const std::string & the_plumed_dat, Log* log_):
251 : KernelBase(path_,the_plumed_dat,log_),
252 4 : log(log_) {
253 : }
254 :
255 6 : ~Kernel() {
256 6 : if(log) {
257 4 : (*log)<<"\n";
258 4 : (*log) <<"Kernel: "<<path<<"\n";
259 4 : (*log) <<"Input: "<<plumed_dat<<"\n";
260 4 : if(comparative_timing>0.0) {
261 2 : (*log).printf("Comparative: %.3f +- %.3f\n",comparative_timing,comparative_timing_error);
262 : }
263 : }
264 6 : }
265 :
266 2 : Kernel(Kernel && other) noexcept:
267 : KernelBase(std::move(other)),
268 2 : log(other.log) {
269 1 : other.log=nullptr; // ensure no log is done in the moved away object
270 : }
271 :
272 2 : Kernel & operator=(Kernel && other) noexcept {
273 2 : if(this != &other) {
274 2 : KernelBase::operator=(std::move(other));
275 2 : log=other.log;
276 2 : other.log=nullptr; // ensure no log is done in the moved away object
277 : }
278 2 : return *this;
279 : }
280 : };
281 :
282 :
283 : class Benchmark:
284 : public CLTool {
285 : public:
286 : static void registerKeywords( Keywords& keys );
287 : explicit Benchmark(const CLToolOptions& co );
288 : int main(FILE* in, FILE*out,Communicator& pc) override;
289 :
290 7 : std::string description()const override {
291 7 : return "run a calculation with a fixed trajectory to find bottlenecks in PLUMED";
292 : }
293 :
294 : //this does the parsing
295 3 : std::optional<std::unique_ptr<AtomDistribution>> parseAtomDistribution(Log& log) {
296 : {
297 3 : std::string trajectoryFile="";
298 : int nn=0;
299 3 : std::string trajectory_fmt="";
300 18 : for (const auto & trj_type : TrajectoryParser::trajectoryOptions()) {
301 : std::string tmp;
302 30 : parse("--i"+trj_type, tmp);
303 15 : if (tmp.length()>0) {
304 1 : log << "Using --i"<<trj_type<<"=" << tmp << "\n";
305 : trajectory_fmt=trj_type;
306 1 : ++nn;
307 : trajectoryFile=tmp;
308 : }
309 3 : }
310 3 : bool use_molfile=false;
311 : #ifdef __PLUMED_HAS_MOLFILE_PLUGINS
312 : {
313 3 : auto plugins_names=TrajectoryParser::getMolfilePluginsnames() ;
314 30 : for(unsigned i=0; i<plugins_names.size(); i++) {
315 27 : std::string molfile_key="--mf_"+plugins_names[i];
316 : std::string traj_molfile;
317 27 : parse(molfile_key,traj_molfile);
318 27 : if(traj_molfile.length()>0) {
319 0 : ++nn;
320 0 : log << "Using --mf_"<<plugins_names[i]<<"=" << traj_molfile << "\n";
321 : trajectoryFile=traj_molfile;
322 : trajectory_fmt=plugins_names[i];
323 0 : use_molfile=true;
324 : }
325 : }
326 3 : }
327 : #endif
328 3 : if(nn>1) {
329 0 : std::fprintf(stderr,"ERROR: cannot provide more than one trajectory file\n");
330 : //let the "main"
331 0 : return std::nullopt;
332 : }
333 3 : if (nn==1) {
334 1 : return std::make_unique<fileTraj>(
335 : trajectory_fmt,
336 : trajectoryFile,
337 : use_molfile,
338 2 : -1
339 : );
340 : }
341 : }
342 : std::string atomicDistr;
343 2 : parse("--atom-distribution",atomicDistr);
344 2 : if(atomicDistr != "") {
345 2 : log << "Using --atom-distribution=" << atomicDistr << "\n";
346 2 : return AtomDistribution::getAtomDistribution(atomicDistr);
347 : }
348 0 : return std::nullopt;
349 : }
350 :
351 : //parse and evenually decorate the AtomDistribution
352 3 : std::optional<std::unique_ptr<AtomDistribution>> createAtomDistribution(Log& log, unsigned nat) {
353 3 : auto toret= parseAtomDistribution(log);
354 3 : if (!toret.has_value()) {
355 0 : return std::nullopt;
356 : }
357 :
358 3 : if(std::string trajectoryModifications="";
359 6 : parse("--modify-trajectory",trajectoryModifications)) {
360 :
361 1 : log << "Using --modify-trajectory=" << trajectoryModifications << "\n";
362 : std::unique_ptr<AtomDistribution> tmp = std::move(*toret);
363 1 : toret = AtomDistribution::decorateAtomDistribution(std::move(tmp),trajectoryModifications);
364 1 : }
365 :
366 : return toret;
367 : }
368 : };
369 :
370 21974 : PLUMED_REGISTER_CLTOOL(Benchmark,"benchmark")
371 :
372 7321 : void Benchmark::registerKeywords( Keywords& keys ) {
373 7321 : CLTool::registerKeywords( keys );
374 7321 : keys.add("compulsory","--plumed","plumed.dat","colon separated path(s) to the input file(s)");
375 7321 : keys.add("compulsory","--kernel","this","colon separated path(s) to kernel(s)");
376 7321 : keys.add("compulsory","--natoms","100000","the number of atoms to use for the simulation");
377 : // Maybe "--natoms" can be more clear when calling --help if we use reset_style to "atoms"
378 7321 : keys.add("compulsory","--nsteps","2000","number of steps of MD to perform (-1 means forever)");
379 7321 : keys.add("compulsory","--maxtime","-1","maximum number of seconds (-1 means forever)");
380 7321 : keys.add("compulsory","--sleep","0","number of seconds of sleep, mimicking MD calculation");
381 7321 : keys.add("compulsory","--atom-distribution","line|wiggle 0.5","the kind of possible atomic displacement at each step");
382 : // Maybe "--atom-distribution" can be more clear when calling --help if we use reset_style to "atoms"
383 7321 : keys.add("optional","--dump-trajectory","dump the trajectory to this file");
384 7321 : keys.addFlag("--ad-help",false,"print a small help text for the atom distributions");
385 7321 : keys.addFlag("--domain-decomposition",false,"simulate domain decomposition, implies --shuffle");
386 7321 : keys.addFlag("--shuffled",false,"reshuffle atoms");
387 7321 : TrajectoryParser::registerKeywords(keys);
388 7321 : keys.add("optional","--modify-trajectory","apply some modifications to the trajectory used in the benchmark, can modify a parsed file");
389 7321 : }
390 :
391 11 : Benchmark::Benchmark(const CLToolOptions& co ):
392 11 : CLTool(co) {
393 11 : inputdata=inputType::commandline;
394 11 : }
395 :
396 :
397 4 : int Benchmark::main(FILE* in, FILE*out,Communicator& pc) {
398 : // deterministic initializations to avoid issues with MPI
399 : generator rng;
400 4 : PLMD::Random atomicGenerator;
401 :
402 : struct FileDeleter {
403 : void operator()(FILE*f) const noexcept {
404 : if(f) {
405 4 : std::fclose(f);
406 : }
407 4 : }
408 : };
409 :
410 4 : std::unique_ptr<FILE,FileDeleter> log_dev_null{std::fopen("/dev/null","w")};
411 :
412 4 : Log log;
413 4 : if(pc.Get_rank()==0) {
414 4 : log.link(out);
415 : } else {
416 0 : log.link(log_dev_null.get());
417 : }
418 :
419 4 : bool printADHelp=false;
420 4 : parseFlag("--ad-help",printADHelp);
421 4 : if (printADHelp) {
422 1 : log << "Documentation for the atomic distributions:\n";
423 1 : log << "\nThese are the basic atomic distributions:\n";
424 1 : auto distrs = AtomDistribution::getDistributionDocumentation();
425 10 : for(const auto& d:distrs) {
426 9 : log.printf("%8s - %s\n\n",d.id.c_str(),d.doc.c_str());
427 : }
428 1 : auto modifiers = AtomDistribution::getDecoratorsDocumentation();
429 :
430 1 : log << "\nThese are the modifiers that can be applied to any atomic distributions:\n";
431 6 : for(const auto& d:modifiers) {
432 5 : log.printf("%8s - %s\n\n",d.id.c_str(),d.doc.c_str());
433 : }
434 1 : log << R"==(
435 : To use these distributions in the benchmark you can simple the cli option --atom-distribution
436 : like:
437 : --atom-distribution="sc", with no extra modifiers
438 : --atom-distribution="line|wiggle 0.5", use the pipe for modify the base distribution
439 : --atom-distribution="ifcc|scale 2.3|reply 2 1 1", do not add spaces before or after the pipes
440 :
441 : As an extra tool you can append more modificators with --modify-trajectory,
442 : it is not necessary since you can use --atom-distribution, but can be applied
443 : to a trajectory read from a file, for example "box" can be used for adding
444 : the pbcs to certain file parsers.
445 : Remember that --modify-trajectory accepts only modifiers, with the
446 : same syntax of --atom-distribution
447 : --modify-trajectory="box 5 5 5"
448 : --modify-trajectory="box 0 1 1 1 1 0 1 0 1|reply 4 2 6"
449 : )==";
450 : return 0;
451 1 : }
452 :
453 3 : log.setLinePrefix("BENCH: ");
454 3 : log <<"Welcome to PLUMED benchmark\n";
455 3 : std::vector<Kernel> kernels;
456 :
457 : // perform comparative analysis
458 : // ensure that kernels vector is destroyed from last to first element upon exit
459 3 : auto kernels_deleter=[&log](auto f) {
460 3 : if(!f) {
461 0 : return;
462 : }
463 3 : if(f->empty()) {
464 : return;
465 : }
466 : generator bootstrapRng;
467 :
468 : const auto size=f->back().timings.size();
469 : //B are the bootstrap iterations
470 : constexpr int B=200;
471 3 : const size_t numblocks=size;
472 : // for some reasons, blocked bootstrap underestimates error
473 : // For now I keep it off. If I remove it, the code below could be simplified
474 : // if(numblocks>20) numblocks=20;
475 : // size is zero if no timings were recorded (e.g. when only dumping the trajectory):
476 : // do not divide by zero (it is a SIGFPE in builds without optimization)
477 4 : const auto blocksize=numblocks>0 ? size/numblocks : 0;
478 :
479 3 : if(f->size()<2) {
480 2 : log<<"Single run, skipping comparative analysis\n";
481 1 : } else if(size<10) {
482 0 : log<<"Too small sample, skipping comparative analysis\n";
483 : } else
484 : try {
485 :
486 1 : log<<"Running comparative analysis, "<<numblocks<<" blocks with size "<<blocksize<<"\n";
487 :
488 1 : std::vector<std::size_t> choice(size);
489 1 : std::uniform_int_distribution<> distrib(0, numblocks-1);
490 1 : std::vector<std::vector<long long int>> blocks(f->size());
491 :
492 : {
493 : int i=0;
494 3 : for(auto it = f->rbegin(); it != f->rend(); ++it,++i) {
495 : size_t l=0;
496 2 : blocks[i].assign(numblocks,0);
497 66 : for(auto j=0ULL; j<numblocks; j++) {
498 128 : for(auto k=0ULL; k<blocksize; k++) {
499 64 : plumed_assert(l<it->timings.size());
500 64 : blocks[i][j]+=it->timings[l];
501 64 : l++;
502 : }
503 : }
504 : }
505 : }
506 :
507 1 : std::vector<std::vector<double>> ratios(f->size());
508 3 : for(auto & r : ratios) {
509 : //B are the bootstrap iterations
510 2 : r.resize(B);
511 : }
512 :
513 : //B are the bootstrap iterations
514 201 : for(unsigned b=0; b<B; b++) {
515 6600 : for(auto & c : choice) {
516 6400 : c=distrib(bootstrapRng);
517 : }
518 : long long int reference=0;
519 6600 : for(const auto & c : choice) {
520 6400 : reference+=blocks[0][c];
521 : }
522 600 : for(auto i=0ULL; i<blocks.size(); i++) {
523 : long long int estimate=0;
524 : // this would lead to separate bootstrap samples for each estimate:
525 : // for(auto & c : choice){c=distrib(bootstrapRng);}
526 13200 : for(const auto & c : choice) {
527 12800 : estimate+=blocks[i][c];
528 : }
529 400 : ratios[i][b]=double(estimate)/double(reference);
530 : }
531 : }
532 :
533 : {
534 : int i=0;
535 3 : for(auto it = f->rbegin(); it != f->rend(); ++it,++i) {
536 : double sum=0.0;
537 : double sum2=0.0;
538 402 : for(auto r : ratios[i]) {
539 400 : sum+=r;
540 400 : sum2+=r*r;
541 : }
542 : //B are the bootstrap iterations
543 2 : it->comparative_timing=sum/B;
544 2 : it->comparative_timing_error=std::sqrt(sum2/B-sum*sum/(B*B));
545 : }
546 : }
547 :
548 1 : } catch(std::exception & e) {
549 0 : log<<"Unexpected error during comparative analysis\n";
550 0 : log<<e.what()<<"\n";
551 : }
552 7 : while(!f->empty()) {
553 : f->pop_back();
554 : }
555 :
556 : };
557 : std::unique_ptr<decltype(kernels),decltype(kernels_deleter)> kernels_deleter_obj(&kernels,kernels_deleter);
558 :
559 :
560 : // construct the kernels vector:
561 : {
562 3 : std::vector<std::string> allpaths;
563 :
564 : {
565 : std::string paths;
566 3 : parse("--kernel",paths);
567 3 : log <<"Using --kernel=" << paths << "\n";
568 6 : allpaths=Tools::getWords(paths,":");
569 : }
570 :
571 3 : std::vector<std::string> allplumed;
572 : {
573 : std::string paths;
574 3 : parse("--plumed",paths);
575 3 : log <<"Using --plumed=" << paths << "\n";
576 6 : allplumed=Tools::getWords(paths,":");
577 : }
578 :
579 3 : plumed_assert(allplumed.size()>0 && allpaths.size()>0);
580 :
581 : // this check only works on MacOS
582 : #if defined(__APPLE__)
583 : // if any of the paths if different from "this", we check if libplumed was loaded locally to avoid conflicts.
584 : if(std::any_of(allpaths.begin(),allpaths.end(),[](auto value) {
585 : return value != "this";
586 : })) {
587 : if(DLLoader::isPlumedGlobal()) {
588 : plumed_error()<<"It looks like libplumed is loaded in the global namespace, you cannot load a different version of the kernel\n"
589 : <<"Please make sure you use the plumed-runtime executable and that the env var PLUMED_LOAD_NAMESPACE is not set to GLOBAL";
590 : }
591 : }
592 : #endif
593 :
594 3 : if(allplumed.size()>1 && allpaths.size()>1 && allplumed.size() != allpaths.size()) {
595 0 : plumed_error() << "--kernel and --plumed should have either one element or the same number of elements";
596 : }
597 :
598 3 : if(allplumed.size()>1 && allpaths.size()==1)
599 2 : for(unsigned i=1; i<allplumed.size(); i++) {
600 1 : allpaths.push_back(allpaths[0]);
601 : }
602 3 : if(allplumed.size()==1 && allpaths.size()>1)
603 0 : for(unsigned i=1; i<allpaths.size(); i++) {
604 0 : allplumed.push_back(allplumed[0]);
605 : }
606 :
607 7 : for(unsigned i=0; i<allpaths.size(); i++) {
608 4 : kernels.emplace_back(allpaths[i],allplumed[i],&log);
609 : }
610 3 : }
611 :
612 : // reverse order so that log happens in the forward order:
613 3 : std::reverse(kernels.begin(),kernels.end());
614 :
615 : // read other flags:
616 3 : bool shuffled=false;
617 3 : parseFlag("--shuffled",shuffled);
618 :
619 : int nf;
620 3 : parse("--nsteps",nf);
621 3 : log << "Using --nsteps=" << nf << "\n";
622 : unsigned natoms;
623 3 : parse("--natoms",natoms);
624 3 : log << "Using --natoms=" << natoms << "\n";
625 : double maxtime;
626 3 : parse("--maxtime",maxtime);
627 3 : log << "Using --maxtime=" << maxtime << "\n";
628 :
629 3 : bool domain_decomposition=false;
630 3 : parseFlag("--domain-decomposition",domain_decomposition);
631 :
632 3 : if(pc.Get_size()>1) {
633 0 : domain_decomposition=true;
634 : }
635 3 : if(domain_decomposition) {
636 1 : shuffled=true;
637 : }
638 :
639 3 : if (shuffled) {
640 2 : log << "Using --shuffled\n";
641 : }
642 3 : if (domain_decomposition) {
643 1 : log << "Using --domain-decomposition\n";
644 : }
645 :
646 : double timeToSleep;
647 3 : parse("--sleep",timeToSleep);
648 3 : log << "Using --sleep=" << timeToSleep << "\n";
649 :
650 3 : std::vector<int> shuffled_indexes;
651 3 : std::unique_ptr<AtomDistribution> distribution;
652 3 : if(auto checkDistr = createAtomDistribution(log,natoms);
653 : checkDistr.has_value()) {
654 : distribution = std::move (*checkDistr);
655 3 : if (distribution->overrideNat(natoms)) {
656 1 : log << "Distribution overrode --natoms, Using --natoms=" << natoms << "\n";
657 : }
658 : } else {
659 0 : std::fprintf(stderr,"ERROR: problem with setting up the trajectory for the benchmark\n");
660 : return 1;
661 : }
662 :
663 : {
664 : std::string fileToDump;
665 6 : if(parse("--dump-trajectory",fileToDump)) {
666 1 : log << "Saving the trajectory to \"" << fileToDump << "\" and exiting\n";
667 1 : std::vector<double> cell(9);
668 1 : std::vector<Vector> pos(natoms);
669 1 : std::ofstream ofile(fileToDump);
670 1 : if (nf<0) {
671 : //if the user accidentally sets infinite steps, we set it to print only one
672 0 : nf=1;
673 : }
674 3 : for(int step=0; step<nf; ++step) {
675 2 : auto sw=kernels[0].stopwatch.startStop("TrajectoryGeneration");
676 2 : distribution->frame(pos,cell,step,atomicGenerator);
677 :
678 2 : ofile << natoms << "\n"
679 2 : << cell[0] << " " << cell[1] << " " << cell[2] << " "
680 2 : << cell[3] << " " << cell[4] << " " << cell[5] << " "
681 2 : << cell[6] << " " << cell[7] << " " << cell[8]
682 2 : << "\n";
683 56 : for(unsigned i=0; i<natoms; ++i) {
684 54 : ofile << "X\t" << pos[i]<< "\n";
685 : }
686 2 : }
687 1 : ofile.close();
688 : return 0;
689 1 : }
690 : }
691 :
692 2 : log <<"Initializing the setup of the kernel(s)\n";
693 2 : const auto initial_time=std::chrono::high_resolution_clock::now();
694 :
695 5 : for(auto & k : kernels) {
696 : auto & p(k.handle);
697 3 : auto sw=k.stopwatch.startStop("A Initialization");
698 : // MPI should also be initialized (it is not with plumed --no-mpi)
699 3 : if(Communicator::initialized() && domain_decomposition) {
700 0 : p.cmd("setMPIComm",&pc.Get_comm());
701 : }
702 3 : p.cmd("setRealPrecision",(int)sizeof(double));
703 3 : p.cmd("setMDLengthUnits",1.0);
704 3 : p.cmd("setMDChargeUnits",1.0);
705 3 : p.cmd("setMDMassUnits",1.0);
706 3 : p.cmd("setMDEngine","benchmarks");
707 3 : p.cmd("setTimestep",1.0);
708 3 : p.cmd("setPlumedDat",k.plumed_dat.c_str());
709 3 : p.cmd("setLog",out);
710 3 : p.cmd("setNatoms",natoms);
711 3 : p.cmd("init");
712 3 : }
713 :
714 2 : std::vector<double> cell( 9 ), virial( 9 );
715 2 : std::vector<Vector> pos( natoms ), forces( natoms );
716 : // positions in the shuffled order, used with --shuffled
717 2 : std::vector<Vector> pos_shuffled;
718 2 : std::vector<double> masses( natoms, 1 ), charges( natoms, 0 );
719 :
720 2 : if(shuffled) {
721 2 : shuffled_indexes.resize(natoms);
722 93 : for(unsigned i=0; i<natoms; i++) {
723 91 : shuffled_indexes[i]=i;
724 : }
725 2 : std::shuffle(shuffled_indexes.begin(),shuffled_indexes.end(),rng);
726 2 : pos_shuffled.resize(natoms);
727 : }
728 :
729 : // non owning pointers, used for shuffling the execution order
730 2 : std::vector<Kernel*> kernels_ptr;
731 5 : for(unsigned i=0; i<kernels.size(); i++) {
732 3 : kernels_ptr.push_back(&kernels[i]);
733 : }
734 :
735 2 : int plumedStopCondition=0;
736 : bool fast_finish=false;
737 : int part=0;
738 :
739 2 : log<<"Starting MD loop\n";
740 2 : log<<"Use CTRL+C to stop at any time and collect timers (not working in MPI runs)\n";
741 : // trap signals:
742 2 : SignalHandlerGuard sigIntGuard(SIGINT, signalHandler);
743 2 : SignalHandlerGuard sigTermGuard(SIGTERM, signalHandler);
744 :
745 45 : for(int step=0; nf<0 || step<nf; ++step) {
746 43 : std::shuffle(kernels_ptr.begin(),kernels_ptr.end(),rng);
747 43 : distribution->frame(pos,cell,step,atomicGenerator);
748 43 : if(shuffled) {
749 : // atom shuffled_indexes[k] is passed in place k, so its position is copied there:
750 : // the structure seen by PLUMED does not depend on the permutation
751 : // (which is not portable, since std::shuffle is implementation defined)
752 2684 : for(unsigned k=0; k<natoms; k++) {
753 2641 : pos_shuffled[k]=pos[shuffled_indexes[k]];
754 : }
755 : }
756 43 : auto & pos_passed(shuffled?pos_shuffled:pos);
757 :
758 : double* pos_ptr;
759 : double* for_ptr;
760 : double* charges_ptr;
761 : double* masses_ptr;
762 : int* indexes_ptr=nullptr;
763 : int n_local_atoms;
764 :
765 43 : if(domain_decomposition) {
766 3 : const auto nproc=pc.Get_size();
767 3 : const auto nn=natoms/nproc;
768 : //using int to remove warning, MPI don't work with unsigned
769 3 : int excess=natoms%nproc;
770 3 : const auto myrank=pc.Get_rank();
771 : auto shift=0;
772 3 : n_local_atoms=nn;
773 3 : if(myrank<excess) {
774 0 : n_local_atoms+=1;
775 : }
776 3 : for(int i=0; i<myrank; i++) {
777 0 : shift+=nn;
778 0 : if(i<excess) {
779 0 : shift+=1;
780 : }
781 : }
782 3 : pos_ptr=&pos_passed[shift][0];
783 : for_ptr=&forces[shift][0];
784 : charges_ptr=&charges[shift];
785 : masses_ptr=&masses[shift];
786 3 : indexes_ptr=shuffled_indexes.data()+shift;
787 : } else {
788 : pos_ptr=&pos_passed[0][0];
789 : for_ptr=&forces[0][0];
790 : charges_ptr=&charges[0];
791 : masses_ptr=&masses[0];
792 40 : n_local_atoms=natoms;
793 : indexes_ptr=shuffled_indexes.data();
794 : }
795 :
796 : const char* sw_name;
797 : if(part==0) {
798 : sw_name="B0 First step";
799 : } else if(part==1) {
800 : sw_name="B1 Warm-up";
801 : } else if(part==2) {
802 : sw_name="B2 Calculation part 1";
803 : } else {
804 : sw_name="B3 Calculation part 2";
805 : }
806 :
807 :
808 126 : for(unsigned i=0; i<kernels_ptr.size(); i++) {
809 83 : auto & p(kernels_ptr[i]->handle);
810 :
811 : {
812 83 : auto sw=kernels_ptr[i]->stopwatch.startPause(sw_name);
813 83 : p.cmd("setStep",step);
814 83 : p.cmd("setStopFlag",&plumedStopCondition);
815 83 : p.cmd("setForces",for_ptr, {n_local_atoms,3});
816 83 : p.cmd("setBox",&cell[0], {3,3});
817 83 : p.cmd("setVirial",&virial[0], {3,3});
818 83 : p.cmd("setPositions",pos_ptr, {n_local_atoms,3});
819 83 : p.cmd("setMasses",masses_ptr, {n_local_atoms});
820 83 : p.cmd("setCharges",charges_ptr, {n_local_atoms});
821 83 : if(shuffled) {
822 83 : p.cmd("setAtomsNlocal",n_local_atoms);
823 83 : p.cmd("setAtomsGatindex",indexes_ptr, {n_local_atoms});
824 : }
825 83 : p.cmd("prepareCalc");
826 83 : }
827 :
828 : // mimick MD calculation here
829 : {
830 : unsigned k=0;
831 83 : auto start=std::chrono::high_resolution_clock::now();
832 83 : while(std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now()-start).count()<(long long int)1e9*timeToSleep) {
833 0 : k+=i*i;
834 : }
835 : std::fprintf(log_dev_null.get(),"%u",k);
836 : }
837 :
838 : {
839 83 : auto sw=kernels_ptr[i]->stopwatch.startStop(sw_name);
840 83 : p.cmd("performCalc");
841 83 : }
842 :
843 83 : if(kernels_ptr.size()>1 && part>1) {
844 64 : kernels_ptr[i]->timings.push_back(kernels_ptr[i]->stopwatch.getLastCycle(sw_name));
845 : }
846 83 : if(plumedStopCondition || signalReceived.load()) {
847 : fast_finish=true;
848 : }
849 : }
850 43 : auto elapsed=std::chrono::duration_cast<std::chrono::nanoseconds>(std::chrono::high_resolution_clock::now()-initial_time).count();
851 43 : if(part==0) {
852 : part=1;
853 : }
854 41 : if(part<2) {
855 9 : if((maxtime>0 && elapsed>(long long int)(0.2*1e9*maxtime)) || (nf>0 && step+1>=nf/5) || (maxtime<0 && nf<0 && step+1>=100)) {
856 : part=2;
857 2 : log<<"Warm-up completed\n";
858 : }
859 : }
860 41 : if(part<3) {
861 25 : if((maxtime>0 && elapsed>(long long int)(0.6*1e9*maxtime)) || (nf>0 && step+1>=3*nf/5)) {
862 : part=3;
863 2 : log<<"60% completed\n";
864 : }
865 : }
866 :
867 43 : if(maxtime>0 && elapsed>(long long int)(1e9*maxtime)) {
868 : fast_finish=true;
869 : }
870 :
871 : {
872 43 : unsigned tmp=fast_finish;
873 43 : pc.Bcast(tmp,0);
874 43 : fast_finish=tmp;
875 : }
876 43 : if(fast_finish) {
877 : break;
878 : }
879 : }
880 :
881 : return 0;
882 4 : }
883 :
884 : } // namespace unnamed
885 : } // namespace cltools
886 : } // namespace PLMD
|
