2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019,2020, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS 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 GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
36 * \brief Defines the modular simulator
38 * \author Pascal Merz <pascal.merz@me.com>
39 * \ingroup module_modularsimulator
44 #include "modularsimulator.h"
46 #include "gromacs/commandline/filenm.h"
47 #include "gromacs/domdec/domdec.h"
48 #include "gromacs/ewald/pme.h"
49 #include "gromacs/ewald/pme_load_balancing.h"
50 #include "gromacs/ewald/pme_pp.h"
51 #include "gromacs/gmxlib/network.h"
52 #include "gromacs/gmxlib/nrnb.h"
53 #include "gromacs/math/vec.h"
54 #include "gromacs/mdlib/checkpointhandler.h"
55 #include "gromacs/mdlib/constr.h"
56 #include "gromacs/mdlib/energyoutput.h"
57 #include "gromacs/mdlib/forcerec.h"
58 #include "gromacs/mdlib/mdatoms.h"
59 #include "gromacs/mdlib/resethandler.h"
60 #include "gromacs/mdlib/stat.h"
61 #include "gromacs/mdlib/update.h"
62 #include "gromacs/mdrun/replicaexchange.h"
63 #include "gromacs/mdrun/shellfc.h"
64 #include "gromacs/mdrunutility/handlerestart.h"
65 #include "gromacs/mdrunutility/printtime.h"
66 #include "gromacs/mdtypes/commrec.h"
67 #include "gromacs/mdtypes/fcdata.h"
68 #include "gromacs/mdtypes/forcerec.h"
69 #include "gromacs/mdtypes/inputrec.h"
70 #include "gromacs/mdtypes/mdatom.h"
71 #include "gromacs/mdtypes/mdrunoptions.h"
72 #include "gromacs/mdtypes/observableshistory.h"
73 #include "gromacs/mdtypes/state.h"
74 #include "gromacs/nbnxm/nbnxm.h"
75 #include "gromacs/timing/walltime_accounting.h"
76 #include "gromacs/topology/mtop_util.h"
77 #include "gromacs/topology/topology.h"
78 #include "gromacs/utility/cstringutil.h"
79 #include "gromacs/utility/fatalerror.h"
81 #include "compositesimulatorelement.h"
82 #include "computeglobalselement.h"
83 #include "constraintelement.h"
84 #include "energyelement.h"
85 #include "forceelement.h"
86 #include "freeenergyperturbationelement.h"
87 #include "parrinellorahmanbarostat.h"
88 #include "propagator.h"
89 #include "shellfcelement.h"
90 #include "signallers.h"
91 #include "statepropagatordata.h"
92 #include "trajectoryelement.h"
93 #include "vrescalethermostat.h"
97 void ModularSimulator::run()
99 GMX_LOG(mdlog.info).asParagraph().appendText("Using the modular simulator.");
100 constructElementsAndSignallers();
102 for (auto& signaller : signallerCallList_)
104 signaller->signallerSetup();
108 domDecHelper_->setup();
111 for (auto& element : elementsOwnershipList_)
113 element->elementSetup();
115 if (pmeLoadBalanceHelper_)
117 // State must have been initialized so pmeLoadBalanceHelper_ gets a valid box
118 pmeLoadBalanceHelper_->setup();
121 while (step_ <= signalHelper_->lastStep_)
125 while (!taskQueue_.empty())
127 auto task = std::move(taskQueue_.front());
134 for (auto& element : elementsOwnershipList_)
136 element->elementTeardown();
138 if (pmeLoadBalanceHelper_)
140 pmeLoadBalanceHelper_->teardown();
145 void ModularSimulator::simulatorSetup()
147 if (!mdrunOptions.writeConfout)
149 // This is on by default, and the main known use case for
150 // turning it off is for convenience in benchmarking, which is
151 // something that should not show up in the general user
156 "The -noconfout functionality is deprecated, and "
157 "may be removed in a future version.");
162 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
163 std::string timeString;
164 fprintf(stderr, "starting mdrun '%s'\n", *(top_global->name));
165 if (inputrec->nsteps >= 0)
167 timeString = formatString("%8.1f", static_cast<double>(inputrec->init_step + inputrec->nsteps)
168 * inputrec->delta_t);
172 timeString = "infinite";
174 if (inputrec->init_step > 0)
176 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
177 gmx_step_str(inputrec->init_step + inputrec->nsteps, sbuf), timeString.c_str(),
178 gmx_step_str(inputrec->init_step, sbuf2), inputrec->init_step * inputrec->delta_t);
182 fprintf(stderr, "%s steps, %s ps.\n", gmx_step_str(inputrec->nsteps, sbuf),
185 fprintf(fplog, "\n");
188 walltime_accounting_start_time(walltime_accounting);
189 wallcycle_start(wcycle, ewcRUN);
190 print_start(fplog, cr, walltime_accounting, "mdrun");
192 step_ = inputrec->init_step;
195 void ModularSimulator::preStep(Step step, Time gmx_unused time, bool isNeighborSearchingStep)
197 if (stopHandler_->stoppingAfterCurrentStep(isNeighborSearchingStep) && step != signalHelper_->lastStep_)
200 * Stop handler wants to stop after the current step, which was
201 * not known when building the current task queue. This happens
202 * e.g. when a stop is signalled by OS. We therefore want to purge
203 * the task queue now, and re-schedule this step as last step.
206 std::queue<SimulatorRunFunctionPtr>().swap(taskQueue_);
212 resetHandler_->setSignal(walltime_accounting);
213 // This is a hack to avoid having to rewrite StopHandler to be a NeighborSearchSignaller
214 // and accept the step as input. Eventually, we want to do that, but currently this would
215 // require introducing NeighborSearchSignaller in the legacy do_md or a lot of code
217 stophandlerIsNSStep_ = isNeighborSearchingStep;
218 stophandlerCurrentStep_ = step;
219 stopHandler_->setSignal();
221 wallcycle_start(wcycle, ewcSTEP);
224 void ModularSimulator::postStep(Step step, Time gmx_unused time)
229 if (do_per_step(step, inputrec->nstlog))
231 if (fflush(fplog) != 0)
233 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
237 const bool do_verbose = mdrunOptions.verbose
238 && (step % mdrunOptions.verboseStepPrintInterval == 0
239 || step == inputrec->init_step || step == signalHelper_->lastStep_);
240 // Print the remaining wall clock time for the run
241 if (MASTER(cr) && (do_verbose || gmx_got_usr_signal())
242 && !(pmeLoadBalanceHelper_ && pmeLoadBalanceHelper_->pmePrinting()))
244 print_time(stderr, walltime_accounting, step, inputrec, cr);
247 double cycles = wallcycle_stop(wcycle, ewcSTEP);
248 if (DOMAINDECOMP(cr) && wcycle)
250 dd_cycles_add(cr->dd, static_cast<float>(cycles), ddCyclStep);
253 resetHandler_->resetCounters(
254 step, step - inputrec->init_step, mdlog, fplog, cr, fr->nbv.get(), nrnb, fr->pmedata,
255 pmeLoadBalanceHelper_ ? pmeLoadBalanceHelper_->loadBalancingObject() : nullptr, wcycle,
256 walltime_accounting);
259 void ModularSimulator::simulatorTeardown()
262 // Stop measuring walltime
263 walltime_accounting_end_time(walltime_accounting);
265 if (!thisRankHasDuty(cr, DUTY_PME))
267 /* Tell the PME only node to finish */
268 gmx_pme_send_finish(cr);
271 walltime_accounting_set_nsteps_done(walltime_accounting, step_ - inputrec->init_step);
274 void ModularSimulator::populateTaskQueue()
276 auto registerRunFunction = std::make_unique<RegisterRunFunction>(
277 [this](SimulatorRunFunctionPtr ptr) { taskQueue_.push(std::move(ptr)); });
279 Time startTime = inputrec->init_t;
280 Time timeStep = inputrec->delta_t;
281 Time time = startTime + step_ * timeStep;
283 // Run an initial call to the signallers
284 for (auto& signaller : signallerCallList_)
286 signaller->signal(step_, time);
289 if (checkpointHelper_)
291 checkpointHelper_->run(step_, time);
294 if (pmeLoadBalanceHelper_)
296 pmeLoadBalanceHelper_->run(step_, time);
300 domDecHelper_->run(step_, time);
305 // local variables for lambda capturing
306 const int step = step_;
307 const bool isNSStep = step == signalHelper_->nextNSStep_;
310 (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
311 [this, step, time, isNSStep]() { preStep(step, time, isNSStep); }));
312 // register elements for step
313 for (auto& element : elementCallList_)
315 element->scheduleTask(step_, time, registerRunFunction);
317 // register post-step
318 (*registerRunFunction)(
319 std::make_unique<SimulatorRunFunction>([this, step, time]() { postStep(step, time); }));
323 time = startTime + step_ * timeStep;
324 for (auto& signaller : signallerCallList_)
326 signaller->signal(step_, time);
328 } while (step_ != signalHelper_->nextNSStep_ && step_ <= signalHelper_->lastStep_);
331 void ModularSimulator::constructElementsAndSignallers()
333 /* When restarting from a checkpoint, it can be appropriate to
334 * initialize ekind from quantities in the checkpoint. Otherwise,
335 * compute_globals must initialize ekind before the simulation
336 * starts/restarts. However, only the master rank knows what was
337 * found in the checkpoint file, so we have to communicate in
338 * order to coordinate the restart.
340 * TODO (modular) This should become obsolete when checkpoint reading
341 * happens within the modular simulator framework: The energy
342 * element should read its data from the checkpoint file pointer,
343 * and signal to the compute globals element if it needs anything
346 * TODO (legacy) Consider removing this communication if/when checkpoint
347 * reading directly follows .tpr reading, because all ranks can
348 * agree on hasReadEkinState at that time.
350 bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
353 gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr);
355 if (hasReadEkinState)
357 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
361 * Build data structures
364 std::make_unique<TopologyHolder>(*top_global, cr, inputrec, fr, mdAtoms, constr, vsite);
366 std::unique_ptr<FreeEnergyPerturbationElement> freeEnergyPerturbationElement = nullptr;
367 FreeEnergyPerturbationElement* freeEnergyPerturbationElementPtr = nullptr;
368 if (inputrec->efep != efepNO)
370 freeEnergyPerturbationElement =
371 std::make_unique<FreeEnergyPerturbationElement>(fplog, inputrec, mdAtoms);
372 freeEnergyPerturbationElementPtr = freeEnergyPerturbationElement.get();
375 auto statePropagatorData = std::make_unique<StatePropagatorData>(
376 top_global->natoms, fplog, cr, state_global, inputrec->nstxout, inputrec->nstvout,
377 inputrec->nstfout, inputrec->nstxout_compressed, fr->nbv->useGpu(),
378 freeEnergyPerturbationElementPtr, topologyHolder_.get(), fr->bMolPBC,
379 mdrunOptions.writeConfout, opt2fn("-c", nfile, fnm), inputrec, mdAtoms->mdatoms());
380 auto statePropagatorDataPtr = compat::make_not_null(statePropagatorData.get());
382 auto energyElement = std::make_unique<EnergyElement>(
383 statePropagatorDataPtr, freeEnergyPerturbationElementPtr, top_global, inputrec, mdAtoms,
384 enerd, ekind, constr, fplog, fcd, mdModulesNotifier, MASTER(cr), observablesHistory,
386 auto energyElementPtr = compat::make_not_null(energyElement.get());
391 const bool simulationsShareState = false;
392 stopHandler_ = stopHandlerBuilder->getStopHandlerMD(
393 compat::not_null<SimulationSignal*>(&signals_[eglsSTOPCOND]), simulationsShareState,
394 MASTER(cr), inputrec->nstlist, mdrunOptions.reproducible, nstglobalcomm_,
395 mdrunOptions.maximumHoursToRun, inputrec->nstlist == 0, fplog, stophandlerCurrentStep_,
396 stophandlerIsNSStep_, walltime_accounting);
399 * Create simulator builders
401 SignallerBuilder<NeighborSearchSignaller> neighborSearchSignallerBuilder;
402 SignallerBuilder<LastStepSignaller> lastStepSignallerBuilder;
403 SignallerBuilder<LoggingSignaller> loggingSignallerBuilder;
404 SignallerBuilder<EnergySignaller> energySignallerBuilder;
405 TrajectoryElementBuilder trajectoryElementBuilder;
408 * Register data structures to signallers
410 trajectoryElementBuilder.registerWriterClient(statePropagatorDataPtr);
411 trajectoryElementBuilder.registerSignallerClient(statePropagatorDataPtr);
412 lastStepSignallerBuilder.registerSignallerClient(statePropagatorDataPtr);
414 trajectoryElementBuilder.registerWriterClient(energyElementPtr);
415 trajectoryElementBuilder.registerSignallerClient(energyElementPtr);
416 energySignallerBuilder.registerSignallerClient(energyElementPtr);
418 // Register the simulator itself to the neighbor search / last step signaller
419 neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
420 lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
423 * Build integrator - this takes care of force calculation, propagation,
424 * constraining, and of the place the statePropagatorData and the energy element
425 * have a full timestep state.
427 // TODO: Make a CheckpointHelperBuilder
428 std::vector<ICheckpointHelperClient*> checkpointClients = { statePropagatorDataPtr, energyElementPtr,
429 freeEnergyPerturbationElementPtr };
430 CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback = nullptr;
432 buildIntegrator(&neighborSearchSignallerBuilder, &energySignallerBuilder,
433 &loggingSignallerBuilder, &trajectoryElementBuilder, &checkpointClients,
434 &checkBondedInteractionsCallback, statePropagatorDataPtr,
435 energyElementPtr, freeEnergyPerturbationElementPtr, hasReadEkinState);
438 * Build infrastructure elements
441 if (PmeLoadBalanceHelper::doPmeLoadBalancing(mdrunOptions, inputrec, fr))
443 pmeLoadBalanceHelper_ = std::make_unique<PmeLoadBalanceHelper>(
444 mdrunOptions.verbose, statePropagatorDataPtr, fplog, cr, mdlog, inputrec, wcycle, fr);
445 neighborSearchSignallerBuilder.registerSignallerClient(
446 compat::make_not_null(pmeLoadBalanceHelper_.get()));
449 if (DOMAINDECOMP(cr))
451 GMX_ASSERT(checkBondedInteractionsCallback,
452 "Domain decomposition needs a callback for check the number of bonded "
454 domDecHelper_ = std::make_unique<DomDecHelper>(
455 mdrunOptions.verbose, mdrunOptions.verboseStepPrintInterval, statePropagatorDataPtr,
456 topologyHolder_.get(), std::move(checkBondedInteractionsCallback), nstglobalcomm_, fplog,
457 cr, mdlog, constr, inputrec, mdAtoms, nrnb, wcycle, fr, vsite, imdSession, pull_work);
458 neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(domDecHelper_.get()));
461 const bool simulationsShareResetCounters = false;
462 resetHandler_ = std::make_unique<ResetHandler>(
463 compat::make_not_null<SimulationSignal*>(&signals_[eglsRESETCOUNTERS]),
464 simulationsShareResetCounters, inputrec->nsteps, MASTER(cr),
465 mdrunOptions.timingOptions.resetHalfway, mdrunOptions.maximumHoursToRun, mdlog, wcycle,
466 walltime_accounting);
469 * Build signaller list
471 * Note that as signallers depend on each others, the order of calling the signallers
472 * matters. It is the responsibility of this builder to ensure that the order is
475 auto energySignaller = energySignallerBuilder.build(
476 inputrec->nstcalcenergy, inputrec->fepvals->nstdhdl, inputrec->nstpcouple);
477 trajectoryElementBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
478 loggingSignallerBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
479 auto trajectoryElement = trajectoryElementBuilder.build(
480 fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec,
481 top_global, oenv, wcycle, startingBehavior, simulationsShareState);
482 loggingSignallerBuilder.registerSignallerClient(compat::make_not_null(trajectoryElement.get()));
484 // Add checkpoint helper here since we need a pointer to the trajectory element and
485 // need to register it with the lastStepSignallerBuilder
486 auto checkpointHandler = std::make_unique<CheckpointHandler>(
487 compat::make_not_null<SimulationSignal*>(&signals_[eglsCHKPT]), simulationsShareState,
488 inputrec->nstlist == 0, MASTER(cr), mdrunOptions.writeConfout,
489 mdrunOptions.checkpointOptions.period);
490 checkpointHelper_ = std::make_unique<CheckpointHelper>(
491 std::move(checkpointClients), std::move(checkpointHandler), inputrec->init_step,
492 trajectoryElement.get(), top_global->natoms, fplog, cr, observablesHistory,
493 walltime_accounting, state_global, mdrunOptions.writeConfout);
494 lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(checkpointHelper_.get()));
496 lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(trajectoryElement.get()));
497 auto loggingSignaller =
498 loggingSignallerBuilder.build(inputrec->nstlog, inputrec->init_step, inputrec->init_t);
499 lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(loggingSignaller.get()));
500 auto lastStepSignaller =
501 lastStepSignallerBuilder.build(inputrec->nsteps, inputrec->init_step, stopHandler_.get());
502 neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(lastStepSignaller.get()));
503 auto neighborSearchSignaller = neighborSearchSignallerBuilder.build(
504 inputrec->nstlist, inputrec->init_step, inputrec->init_t);
506 addToCallListAndMove(std::move(neighborSearchSignaller), signallerCallList_, signallersOwnershipList_);
507 addToCallListAndMove(std::move(lastStepSignaller), signallerCallList_, signallersOwnershipList_);
508 addToCallListAndMove(std::move(loggingSignaller), signallerCallList_, signallersOwnershipList_);
509 addToCallList(trajectoryElement, signallerCallList_);
510 addToCallListAndMove(std::move(energySignaller), signallerCallList_, signallersOwnershipList_);
513 * Build the element list
515 * This is the actual sequence of (non-infrastructure) elements to be run.
516 * For NVE, the trajectory element is used outside of the integrator
517 * (composite) element, as well as the checkpoint helper. The checkpoint
518 * helper should be on top of the loop, and is only part of the simulator
519 * call list to be able to react to the last step being signalled.
521 addToCallList(checkpointHelper_, elementCallList_);
522 if (freeEnergyPerturbationElement)
524 addToCallListAndMove(std::move(freeEnergyPerturbationElement), elementCallList_,
525 elementsOwnershipList_);
527 addToCallListAndMove(std::move(integrator), elementCallList_, elementsOwnershipList_);
528 addToCallListAndMove(std::move(trajectoryElement), elementCallList_, elementsOwnershipList_);
529 // for vv, we need to setup statePropagatorData after the compute
530 // globals so that we reset the right velocities
531 // TODO: Avoid this by getting rid of the need of resetting velocities in vv
532 elementsOwnershipList_.emplace_back(std::move(statePropagatorData));
533 elementsOwnershipList_.emplace_back(std::move(energyElement));
536 std::unique_ptr<ISimulatorElement>
537 ModularSimulator::buildForces(SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
538 SignallerBuilder<EnergySignaller>* energySignallerBuilder,
539 StatePropagatorData* statePropagatorDataPtr,
540 EnergyElement* energyElementPtr,
541 FreeEnergyPerturbationElement* freeEnergyPerturbationElement)
543 const bool isVerbose = mdrunOptions.verbose;
544 const bool isDynamicBox = inputrecDynamicBox(inputrec);
545 // Check for polarizable models and flexible constraints
546 if (ShellFCElement::doShellsOrFlexConstraints(topologyHolder_->globalTopology(),
547 constr ? constr->numFlexibleConstraints() : 0))
549 auto shellFCElement = std::make_unique<ShellFCElement>(
550 statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElement, isVerbose,
551 isDynamicBox, fplog, cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle, runScheduleWork, vsite,
552 imdSession, pull_work, constr, &topologyHolder_->globalTopology(), enforcedRotation);
553 topologyHolder_->registerClient(shellFCElement.get());
554 neighborSearchSignallerBuilder->registerSignallerClient(
555 compat::make_not_null(shellFCElement.get()));
556 energySignallerBuilder->registerSignallerClient(compat::make_not_null(shellFCElement.get()));
558 // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
559 return std::move(shellFCElement);
563 auto forceElement = std::make_unique<ForceElement>(
564 statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElement,
565 isDynamicBox, fplog, cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle, runScheduleWork,
566 vsite, imdSession, pull_work, enforcedRotation);
567 topologyHolder_->registerClient(forceElement.get());
568 neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
569 energySignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
571 // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
572 return std::move(forceElement);
576 std::unique_ptr<ISimulatorElement> ModularSimulator::buildIntegrator(
577 SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
578 SignallerBuilder<EnergySignaller>* energySignallerBuilder,
579 SignallerBuilder<LoggingSignaller>* loggingSignallerBuilder,
580 TrajectoryElementBuilder* trajectoryElementBuilder,
581 std::vector<ICheckpointHelperClient*>* checkpointClients,
582 CheckBondedInteractionsCallbackPtr* checkBondedInteractionsCallback,
583 compat::not_null<StatePropagatorData*> statePropagatorDataPtr,
584 compat::not_null<EnergyElement*> energyElementPtr,
585 FreeEnergyPerturbationElement* freeEnergyPerturbationElementPtr,
586 bool hasReadEkinState)
589 buildForces(neighborSearchSignallerBuilder, energySignallerBuilder,
590 statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr);
592 // list of elements owned by the simulator composite object
593 std::vector<std::unique_ptr<ISimulatorElement>> elementsOwnershipList;
594 // call list of the simulator composite object
595 std::vector<compat::not_null<ISimulatorElement*>> elementCallList;
597 std::function<void()> needToCheckNumberOfBondedInteractions;
598 if (inputrec->eI == eiMD)
600 auto computeGlobalsElement =
601 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>(
602 statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
603 &signals_, nstglobalcomm_, fplog, mdlog, cr, inputrec, mdAtoms, nrnb,
604 wcycle, fr, &topologyHolder_->globalTopology(), constr, hasReadEkinState);
605 topologyHolder_->registerClient(computeGlobalsElement.get());
606 energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
607 trajectoryElementBuilder->registerSignallerClient(
608 compat::make_not_null(computeGlobalsElement.get()));
610 *checkBondedInteractionsCallback =
611 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
613 auto propagator = std::make_unique<Propagator<IntegrationStep::LeapFrog>>(
614 inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
616 addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
617 addToCallList(statePropagatorDataPtr, elementCallList); // we have a full microstate at time t here!
618 if (inputrec->etc == etcVRESCALE)
620 // TODO: With increased complexity of the propagator, this will need further development,
621 // e.g. using propagators templated for velocity propagation policies and a builder
622 propagator->setNumVelocityScalingVariables(inputrec->opts.ngtc);
623 auto thermostat = std::make_unique<VRescaleThermostat>(
624 inputrec->nsttcouple, -1, false, inputrec->ld_seed, inputrec->opts.ngtc,
625 inputrec->delta_t * inputrec->nsttcouple, inputrec->opts.ref_t, inputrec->opts.tau_t,
626 inputrec->opts.nrdf, energyElementPtr, propagator->viewOnVelocityScaling(),
627 propagator->velocityScalingCallback(), state_global, cr, inputrec->bContinuation);
628 checkpointClients->emplace_back(thermostat.get());
629 energyElementPtr->setVRescaleThermostat(thermostat.get());
630 addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
633 std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
634 if (inputrec->epc == epcPARRINELLORAHMAN)
636 // Building the PR barostat here since it needs access to the propagator
637 // and we want to be able to move the propagator object
638 prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
639 inputrec->nstpcouple, -1, inputrec->delta_t * inputrec->nstpcouple,
640 inputrec->init_step, propagator->viewOnPRScalingMatrix(),
641 propagator->prScalingCallback(), statePropagatorDataPtr, energyElementPtr,
642 fplog, inputrec, mdAtoms, state_global, cr, inputrec->bContinuation);
643 energyElementPtr->setParrinelloRahamnBarostat(prBarostat.get());
644 checkpointClients->emplace_back(prBarostat.get());
646 addToCallListAndMove(std::move(propagator), elementCallList, elementsOwnershipList);
649 auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
650 constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
651 MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
652 auto constraintElementPtr = compat::make_not_null(constraintElement.get());
653 energySignallerBuilder->registerSignallerClient(constraintElementPtr);
654 trajectoryElementBuilder->registerSignallerClient(constraintElementPtr);
655 loggingSignallerBuilder->registerSignallerClient(constraintElementPtr);
657 addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
660 addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
661 addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
664 addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
667 else if (inputrec->eI == eiVV)
669 auto computeGlobalsElement =
670 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>(
671 statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
672 &signals_, nstglobalcomm_, fplog, mdlog, cr, inputrec, mdAtoms, nrnb,
673 wcycle, fr, &topologyHolder_->globalTopology(), constr, hasReadEkinState);
674 topologyHolder_->registerClient(computeGlobalsElement.get());
675 energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
676 trajectoryElementBuilder->registerSignallerClient(
677 compat::make_not_null(computeGlobalsElement.get()));
679 *checkBondedInteractionsCallback =
680 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
682 auto propagatorVelocities = std::make_unique<Propagator<IntegrationStep::VelocitiesOnly>>(
683 inputrec->delta_t * 0.5, statePropagatorDataPtr, mdAtoms, wcycle);
684 auto propagatorVelocitiesAndPositions =
685 std::make_unique<Propagator<IntegrationStep::VelocityVerletPositionsAndVelocities>>(
686 inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
688 addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
690 std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
691 if (inputrec->epc == epcPARRINELLORAHMAN)
693 // Building the PR barostat here since it needs access to the propagator
694 // and we want to be able to move the propagator object
695 prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
696 inputrec->nstpcouple, -1, inputrec->delta_t * inputrec->nstpcouple,
697 inputrec->init_step, propagatorVelocities->viewOnPRScalingMatrix(),
698 propagatorVelocities->prScalingCallback(), statePropagatorDataPtr, energyElementPtr,
699 fplog, inputrec, mdAtoms, state_global, cr, inputrec->bContinuation);
700 energyElementPtr->setParrinelloRahamnBarostat(prBarostat.get());
701 checkpointClients->emplace_back(prBarostat.get());
703 addToCallListAndMove(std::move(propagatorVelocities), elementCallList, elementsOwnershipList);
706 auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Velocities>>(
707 constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
708 MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
709 energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
710 trajectoryElementBuilder->registerSignallerClient(
711 compat::make_not_null(constraintElement.get()));
712 loggingSignallerBuilder->registerSignallerClient(
713 compat::make_not_null(constraintElement.get()));
715 addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
717 addToCallList(compat::make_not_null(computeGlobalsElement.get()), elementCallList);
718 addToCallList(statePropagatorDataPtr, elementCallList); // we have a full microstate at time t here!
719 if (inputrec->etc == etcVRESCALE)
721 // TODO: With increased complexity of the propagator, this will need further development,
722 // e.g. using propagators templated for velocity propagation policies and a builder
723 propagatorVelocitiesAndPositions->setNumVelocityScalingVariables(inputrec->opts.ngtc);
724 auto thermostat = std::make_unique<VRescaleThermostat>(
725 inputrec->nsttcouple, 0, true, inputrec->ld_seed, inputrec->opts.ngtc,
726 inputrec->delta_t * inputrec->nsttcouple, inputrec->opts.ref_t,
727 inputrec->opts.tau_t, inputrec->opts.nrdf, energyElementPtr,
728 propagatorVelocitiesAndPositions->viewOnVelocityScaling(),
729 propagatorVelocitiesAndPositions->velocityScalingCallback(), state_global, cr,
730 inputrec->bContinuation);
731 checkpointClients->emplace_back(thermostat.get());
732 energyElementPtr->setVRescaleThermostat(thermostat.get());
733 addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
735 addToCallListAndMove(std::move(propagatorVelocitiesAndPositions), elementCallList,
736 elementsOwnershipList);
739 auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
740 constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
741 MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
742 energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
743 trajectoryElementBuilder->registerSignallerClient(
744 compat::make_not_null(constraintElement.get()));
745 loggingSignallerBuilder->registerSignallerClient(
746 compat::make_not_null(constraintElement.get()));
748 addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
750 addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
751 addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
754 addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
759 gmx_fatal(FARGS, "Integrator not implemented for the modular simulator.");
762 auto integrator = std::make_unique<CompositeSimulatorElement>(std::move(elementCallList),
763 std::move(elementsOwnershipList));
764 // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
765 return std::move(integrator);
768 bool ModularSimulator::isInputCompatible(bool exitOnFailure,
769 const t_inputrec* inputrec,
771 const gmx_mtop_t& globalTopology,
772 const gmx_multisim_t* ms,
773 const ReplicaExchangeParameters& replExParams,
775 bool doEssentialDynamics,
778 auto conditionalAssert = [exitOnFailure](bool condition, const char* message) {
781 GMX_RELEASE_ASSERT(condition, message);
786 bool isInputCompatible = true;
788 // GMX_USE_MODULAR_SIMULATOR allows to use modular simulator also for non-standard uses,
789 // such as the leap-frog integrator
790 const auto modularSimulatorExplicitlyTurnedOn = (getenv("GMX_USE_MODULAR_SIMULATOR") != nullptr);
791 // GMX_USE_MODULAR_SIMULATOR allows to use disable modular simulator for all uses,
792 // including the velocity-verlet integrator used by default
793 const auto modularSimulatorExplicitlyTurnedOff = (getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
796 !(modularSimulatorExplicitlyTurnedOn && modularSimulatorExplicitlyTurnedOff),
797 "Cannot have both GMX_USE_MODULAR_SIMULATOR=ON and GMX_DISABLE_MODULAR_SIMULATOR=ON. "
798 "Unset one of the two environment variables to explicitly chose which simulator to "
800 "or unset both to recover default behavior.");
803 !(modularSimulatorExplicitlyTurnedOff && inputrec->eI == eiVV
804 && inputrec->epc == epcPARRINELLORAHMAN),
805 "Cannot use a Parrinello-Rahman barostat with md-vv and "
806 "GMX_DISABLE_MODULAR_SIMULATOR=ON, "
807 "as the Parrinello-Rahman barostat is not implemented in the legacy simulator. Unset "
808 "GMX_DISABLE_MODULAR_SIMULATOR or use a different pressure control algorithm.");
812 && conditionalAssert(
813 inputrec->eI == eiMD || inputrec->eI == eiVV,
814 "Only integrators md and md-vv are supported by the modular simulator.");
815 isInputCompatible = isInputCompatible
816 && conditionalAssert(inputrec->eI != eiMD || modularSimulatorExplicitlyTurnedOn,
817 "Set GMX_USE_MODULAR_SIMULATOR=ON to use the modular "
818 "simulator with integrator md.");
821 && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
824 && conditionalAssert(
825 inputrec->etc == etcNO || inputrec->etc == etcVRESCALE,
826 "Only v-rescale thermostat is supported by the modular simulator.");
829 && conditionalAssert(
830 inputrec->epc == epcNO || inputrec->epc == epcPARRINELLORAHMAN,
831 "Only Parrinello-Rahman barostat is supported by the modular simulator.");
834 && conditionalAssert(
835 !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec)
836 || inputrecNvtTrotter(inputrec)),
837 "Legacy Trotter decomposition is not supported by the modular simulator.");
838 isInputCompatible = isInputCompatible
839 && conditionalAssert(inputrec->efep == efepNO || inputrec->efep == efepYES
840 || inputrec->efep == efepSLOWGROWTH,
841 "Expanded ensemble free energy calculation is not "
842 "supported by the modular simulator.");
843 isInputCompatible = isInputCompatible
844 && conditionalAssert(!inputrec->bPull,
845 "Pulling is not supported by the modular simulator.");
848 && conditionalAssert(inputrec->opts.ngacc == 1 && inputrec->opts.acc[0][XX] == 0.0
849 && inputrec->opts.acc[0][YY] == 0.0
850 && inputrec->opts.acc[0][ZZ] == 0.0 && inputrec->cos_accel == 0.0,
851 "Acceleration is not supported by the modular simulator.");
854 && conditionalAssert(inputrec->opts.ngfrz == 1 && inputrec->opts.nFreeze[0][XX] == 0
855 && inputrec->opts.nFreeze[0][YY] == 0
856 && inputrec->opts.nFreeze[0][ZZ] == 0,
857 "Freeze groups are not supported by the modular simulator.");
860 && conditionalAssert(
861 inputrec->deform[XX][XX] == 0.0 && inputrec->deform[XX][YY] == 0.0
862 && inputrec->deform[XX][ZZ] == 0.0 && inputrec->deform[YY][XX] == 0.0
863 && inputrec->deform[YY][YY] == 0.0 && inputrec->deform[YY][ZZ] == 0.0
864 && inputrec->deform[ZZ][XX] == 0.0 && inputrec->deform[ZZ][YY] == 0.0
865 && inputrec->deform[ZZ][ZZ] == 0.0,
866 "Deformation is not supported by the modular simulator.");
869 && conditionalAssert(gmx_mtop_interaction_count(globalTopology, IF_VSITE) == 0,
870 "Virtual sites are not supported by the modular simulator.");
871 isInputCompatible = isInputCompatible
872 && conditionalAssert(!inputrec->bDoAwh,
873 "AWH is not supported by the modular simulator.");
876 && conditionalAssert(gmx_mtop_ftype_count(globalTopology, F_DISRES) == 0,
877 "Distance restraints are not supported by the modular simulator.");
880 && conditionalAssert(
881 gmx_mtop_ftype_count(globalTopology, F_ORIRES) == 0,
882 "Orientation restraints are not supported by the modular simulator.");
885 && conditionalAssert(ms == nullptr,
886 "Multi-sim are not supported by the modular simulator.");
889 && conditionalAssert(replExParams.exchangeInterval == 0,
890 "Replica exchange is not supported by the modular simulator.");
892 int numEnsembleRestraintSystems;
895 numEnsembleRestraintSystems = fcd->disres.nsystems;
899 auto distantRestraintEnsembleEnvVar = getenv("GMX_DISRE_ENSEMBLE_SIZE");
900 numEnsembleRestraintSystems =
901 (ms != nullptr && distantRestraintEnsembleEnvVar != nullptr)
902 ? static_cast<int>(strtol(distantRestraintEnsembleEnvVar, nullptr, 10))
907 && conditionalAssert(numEnsembleRestraintSystems <= 1,
908 "Ensemble restraints are not supported by the modular simulator.");
911 && conditionalAssert(!doSimulatedAnnealing(inputrec),
912 "Simulated annealing is not supported by the modular simulator.");
915 && conditionalAssert(!inputrec->bSimTemp,
916 "Simulated tempering is not supported by the modular simulator.");
917 isInputCompatible = isInputCompatible
918 && conditionalAssert(!inputrec->bExpanded,
919 "Expanded ensemble simulations are not supported by "
920 "the modular simulator.");
923 && conditionalAssert(!doEssentialDynamics,
924 "Essential dynamics is not supported by the modular simulator.");
925 isInputCompatible = isInputCompatible
926 && conditionalAssert(inputrec->eSwapCoords == eswapNO,
927 "Ion / water position swapping is not supported by "
928 "the modular simulator.");
931 && conditionalAssert(!inputrec->bIMD,
932 "Interactive MD is not supported by the modular simulator.");
935 && conditionalAssert(!doMembed,
936 "Membrane embedding is not supported by the modular simulator.");
937 const bool useGraph = !areMoleculesDistributedOverPbc(*inputrec, globalTopology, MDLogger());
940 && conditionalAssert(!useGraph, "Graph is not supported by the modular simulator.");
941 // TODO: Change this to the boolean passed when we merge the user interface change for the GPU update.
944 && conditionalAssert(
945 getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
946 "Integration on the GPU is not supported by the modular simulator.");
947 // Modular simulator is centered around NS updates
948 // TODO: think how to handle nstlist == 0
949 isInputCompatible = isInputCompatible
950 && conditionalAssert(inputrec->nstlist != 0,
951 "Simulations without neighbor list update are not "
952 "supported by the modular simulator.");
953 isInputCompatible = isInputCompatible
954 && conditionalAssert(!GMX_FAHCORE,
955 "GMX_FAHCORE not supported by the modular simulator.");
957 return isInputCompatible;
960 void ModularSimulator::checkInputForDisabledFunctionality()
962 isInputCompatible(true, inputrec, doRerun, *top_global, ms, replExParams, fcd,
963 opt2bSet("-ei", nfile, fnm), membed != nullptr);
964 if (observablesHistory->edsamHistory)
967 "The checkpoint is from a run with essential dynamics sampling, "
968 "but the current run did not specify the -ei option. "
969 "Either specify the -ei option to mdrun, or do not use this checkpoint file.");
973 SignallerCallbackPtr ModularSimulator::SignalHelper::registerLastStepCallback()
975 return std::make_unique<SignallerCallback>(
976 [this](Step step, Time gmx_unused time) { this->lastStep_ = step; });
979 SignallerCallbackPtr ModularSimulator::SignalHelper::registerNSCallback()
981 return std::make_unique<SignallerCallback>(
982 [this](Step step, Time gmx_unused time) { this->nextNSStep_ = step; });