2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2019, 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_mdrun
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/gmxlib/network.h"
51 #include "gromacs/gmxlib/nrnb.h"
52 #include "gromacs/mdlib/checkpointhandler.h"
53 #include "gromacs/mdlib/constr.h"
54 #include "gromacs/mdlib/energyoutput.h"
55 #include "gromacs/mdlib/mdatoms.h"
56 #include "gromacs/mdlib/resethandler.h"
57 #include "gromacs/mdlib/stat.h"
58 #include "gromacs/mdlib/update.h"
59 #include "gromacs/mdrun/replicaexchange.h"
60 #include "gromacs/mdrun/shellfc.h"
61 #include "gromacs/mdrunutility/handlerestart.h"
62 #include "gromacs/mdrunutility/printtime.h"
63 #include "gromacs/mdtypes/commrec.h"
64 #include "gromacs/mdtypes/fcdata.h"
65 #include "gromacs/mdtypes/inputrec.h"
66 #include "gromacs/mdtypes/mdrunoptions.h"
67 #include "gromacs/mdtypes/observableshistory.h"
68 #include "gromacs/mdtypes/state.h"
69 #include "gromacs/nbnxm/nbnxm.h"
70 #include "gromacs/timing/walltime_accounting.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/utility/cstringutil.h"
73 #include "gromacs/utility/fatalerror.h"
75 #include "compositesimulatorelement.h"
76 #include "computeglobalselement.h"
77 #include "constraintelement.h"
78 #include "energyelement.h"
79 #include "forceelement.h"
80 #include "propagator.h"
81 #include "shellfcelement.h"
82 #include "signallers.h"
83 #include "statepropagatordata.h"
84 #include "trajectoryelement.h"
88 void ModularSimulator::run()
90 GMX_LOG(mdlog.info).asParagraph().
91 appendText("Using the modular simulator.");
92 constructElementsAndSignallers();
94 for (auto &signaller : signallerCallList_)
96 signaller->signallerSetup();
98 if (pmeLoadBalanceHelper_)
100 pmeLoadBalanceHelper_->setup();
104 domDecHelper_->setup();
107 for (auto &element : elementsOwnershipList_)
109 element->elementSetup();
112 while (step_ <= signalHelper_->lastStep_)
116 while (!taskQueue_.empty())
118 auto task = std::move(taskQueue_.front());
125 for (auto &element : elementsOwnershipList_)
127 element->elementTeardown();
129 if (pmeLoadBalanceHelper_)
131 pmeLoadBalanceHelper_->teardown();
136 void ModularSimulator::simulatorSetup()
138 if (!mdrunOptions.writeConfout)
140 // This is on by default, and the main known use case for
141 // turning it off is for convenience in benchmarking, which is
142 // something that should not show up in the general user
144 GMX_LOG(mdlog.info).asParagraph().
145 appendText("The -noconfout functionality is deprecated, and "
146 "may be removed in a future version.");
151 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
152 std::string timeString;
153 fprintf(stderr, "starting mdrun '%s'\n",
154 *(top_global->name));
155 if (inputrec->nsteps >= 0)
157 timeString = formatString(
158 "%8.1f", static_cast<double>(inputrec->init_step+inputrec->nsteps)*inputrec->delta_t);
162 timeString = "infinite";
164 if (inputrec->init_step > 0)
166 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
167 gmx_step_str(inputrec->init_step+inputrec->nsteps, sbuf),
169 gmx_step_str(inputrec->init_step, sbuf2),
170 inputrec->init_step*inputrec->delta_t);
174 fprintf(stderr, "%s steps, %s ps.\n",
175 gmx_step_str(inputrec->nsteps, sbuf), timeString.c_str());
177 fprintf(fplog, "\n");
180 walltime_accounting_start_time(walltime_accounting);
181 wallcycle_start(wcycle, ewcRUN);
182 print_start(fplog, cr, walltime_accounting, "mdrun");
184 step_ = inputrec->init_step;
187 void ModularSimulator::preStep(
188 Step step, Time gmx_unused time,
189 bool isNeighborSearchingStep)
191 if (stopHandler_->stoppingAfterCurrentStep(isNeighborSearchingStep) &&
192 step != signalHelper_->lastStep_)
195 * Stop handler wants to stop after the current step, which was
196 * not known when building the current task queue. This happens
197 * e.g. when a stop is signalled by OS. We therefore want to purge
198 * the task queue now, and re-schedule this step as last step.
201 std::queue<SimulatorRunFunctionPtr>().swap(taskQueue_);
207 resetHandler_->setSignal(walltime_accounting);
208 // This is a hack to avoid having to rewrite StopHandler to be a NeighborSearchSignaller
209 // and accept the step as input. Eventually, we want to do that, but currently this would
210 // require introducing NeighborSearchSignaller in the legacy do_md or a lot of code
212 stophandlerIsNSStep_ = isNeighborSearchingStep;
213 stophandlerCurrentStep_ = step;
214 stopHandler_->setSignal();
216 wallcycle_start(wcycle, ewcSTEP);
219 void ModularSimulator::postStep(Step step, Time gmx_unused time)
224 if (do_per_step(step, inputrec->nstlog))
226 if (fflush(fplog) != 0)
228 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
232 const bool do_verbose = mdrunOptions.verbose &&
233 (step % mdrunOptions.verboseStepPrintInterval == 0 ||
234 step == inputrec->init_step || step == signalHelper_->lastStep_);
235 // Print the remaining wall clock time for the run
237 (do_verbose || gmx_got_usr_signal()) &&
238 !(pmeLoadBalanceHelper_ && pmeLoadBalanceHelper_->pmePrinting()))
240 print_time(stderr, walltime_accounting, step, inputrec, cr);
243 double cycles = wallcycle_stop(wcycle, ewcSTEP);
244 if (DOMAINDECOMP(cr) && wcycle)
246 dd_cycles_add(cr->dd, static_cast<float>(cycles), ddCyclStep);
249 resetHandler_->resetCounters(
250 step, step - inputrec->init_step, mdlog, fplog, cr, fr->nbv.get(),
252 pmeLoadBalanceHelper_ ? pmeLoadBalanceHelper_->loadBalancingObject() : nullptr,
253 wcycle, walltime_accounting);
256 void ModularSimulator::simulatorTeardown()
259 // Stop measuring walltime
260 walltime_accounting_end_time(walltime_accounting);
262 if (!thisRankHasDuty(cr, DUTY_PME))
264 /* Tell the PME only node to finish */
265 gmx_pme_send_finish(cr);
268 walltime_accounting_set_nsteps_done(walltime_accounting, step_ - inputrec->init_step);
271 void ModularSimulator::populateTaskQueue()
273 auto registerRunFunction = std::make_unique<RegisterRunFunction>(
274 [this](SimulatorRunFunctionPtr ptr)
275 {taskQueue_.push(std::move(ptr)); });
277 Time startTime = inputrec->init_t;
278 Time timeStep = inputrec->delta_t;
279 Time time = startTime + step_*timeStep;
281 // Run an initial call to the signallers
282 for (auto &signaller : signallerCallList_)
284 signaller->signal(step_, time);
287 if (checkpointHelper_)
289 checkpointHelper_->run(step_, time);
292 if (pmeLoadBalanceHelper_)
294 pmeLoadBalanceHelper_->run(step_, time);
298 domDecHelper_->run(step_, time);
303 // local variables for lambda capturing
304 const int step = step_;
305 const bool isNSStep = step == signalHelper_->nextNSStep_;
308 (*registerRunFunction)(
309 std::make_unique<SimulatorRunFunction>(
310 [this, step, time, isNSStep](){preStep(step, time, isNSStep); }));
311 // register elements for step
312 for (auto &element : elementCallList_)
314 element->scheduleTask(step_, time, registerRunFunction);
316 // register post-step
317 (*registerRunFunction)(
318 std::make_unique<SimulatorRunFunction>(
319 [this, step, time](){postStep(step, time); }));
323 time = startTime + step_*timeStep;
324 for (auto &signaller : signallerCallList_)
326 signaller->signal(step_, time);
329 while (step_ != signalHelper_->nextNSStep_ && step_ <= signalHelper_->lastStep_);
332 void ModularSimulator::constructElementsAndSignallers()
334 /* When restarting from a checkpoint, it can be appropriate to
335 * initialize ekind from quantities in the checkpoint. Otherwise,
336 * compute_globals must initialize ekind before the simulation
337 * starts/restarts. However, only the master rank knows what was
338 * found in the checkpoint file, so we have to communicate in
339 * order to coordinate the restart.
341 * TODO (modular) This should become obsolete when checkpoint reading
342 * happens within the modular simulator framework: The energy
343 * element should read its data from the checkpoint file pointer,
344 * and signal to the compute globals element if it needs anything
347 * TODO (legacy) Consider removing this communication if/when checkpoint
348 * reading directly follows .tpr reading, because all ranks can
349 * agree on hasReadEkinState at that time.
351 bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
354 gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr);
356 if (hasReadEkinState)
358 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
362 * Build data structures
364 auto statePropagatorData = std::make_unique<StatePropagatorData>(
365 top_global->natoms, fplog, cr, state_global,
366 inputrec->nstxout, inputrec->nstvout,
367 inputrec->nstfout, inputrec->nstxout_compressed,
368 fr->nbv->useGpu(), inputrec, mdAtoms->mdatoms());
369 auto statePropagatorDataPtr = compat::make_not_null(statePropagatorData.get());
371 auto energyElement = std::make_unique<EnergyElement>(
372 statePropagatorDataPtr, top_global, inputrec, mdAtoms, enerd, ekind,
373 constr, fplog, fcd, mdModulesNotifier, MASTER(cr), observablesHistory, startingBehavior);
374 auto energyElementPtr = compat::make_not_null(energyElement.get());
376 topologyHolder_ = std::make_unique<TopologyHolder>(
377 *top_global, cr, inputrec, fr,
378 mdAtoms, constr, vsite);
383 const bool simulationsShareState = false;
384 stopHandler_ = stopHandlerBuilder->getStopHandlerMD(
385 compat::not_null<SimulationSignal*>(&signals_[eglsSTOPCOND]),
386 simulationsShareState, MASTER(cr), inputrec->nstlist, mdrunOptions.reproducible,
387 nstglobalcomm_, mdrunOptions.maximumHoursToRun, inputrec->nstlist == 0, fplog,
388 stophandlerCurrentStep_, stophandlerIsNSStep_, walltime_accounting);
391 * Create simulator builders
393 SignallerBuilder<NeighborSearchSignaller> neighborSearchSignallerBuilder;
394 SignallerBuilder<LastStepSignaller> lastStepSignallerBuilder;
395 SignallerBuilder<LoggingSignaller> loggingSignallerBuilder;
396 SignallerBuilder<EnergySignaller> energySignallerBuilder;
397 TrajectoryElementBuilder trajectoryElementBuilder;
400 * Register data structures to signallers
402 trajectoryElementBuilder.registerWriterClient(statePropagatorDataPtr);
403 trajectoryElementBuilder.registerSignallerClient(statePropagatorDataPtr);
405 trajectoryElementBuilder.registerWriterClient(energyElementPtr);
406 trajectoryElementBuilder.registerSignallerClient(energyElementPtr);
407 energySignallerBuilder.registerSignallerClient(energyElementPtr);
408 loggingSignallerBuilder.registerSignallerClient(energyElementPtr);
410 // Register the simulator itself to the neighbor search / last step signaller
411 neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
412 lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
415 * Build integrator - this takes care of force calculation, propagation,
416 * constraining, and of the place the statePropagatorData and the energy element
417 * have a full timestep state.
419 CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback = nullptr;
420 auto integrator = buildIntegrator(
421 &neighborSearchSignallerBuilder,
422 &energySignallerBuilder,
423 &loggingSignallerBuilder,
424 &trajectoryElementBuilder,
425 &checkBondedInteractionsCallback,
426 statePropagatorDataPtr,
431 * Build infrastructure elements
434 if (PmeLoadBalanceHelper::doPmeLoadBalancing(mdrunOptions, inputrec, fr))
436 pmeLoadBalanceHelper_ = std::make_unique<PmeLoadBalanceHelper>(
437 mdrunOptions.verbose, statePropagatorDataPtr, fplog,
438 cr, mdlog, inputrec, wcycle, fr);
439 neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(pmeLoadBalanceHelper_.get()));
442 if (DOMAINDECOMP(cr))
445 checkBondedInteractionsCallback,
446 "Domain decomposition needs a callback for check the number of bonded interactions.");
447 domDecHelper_ = std::make_unique<DomDecHelper>(
448 mdrunOptions.verbose, mdrunOptions.verboseStepPrintInterval,
449 statePropagatorDataPtr, topologyHolder_.get(), std::move(checkBondedInteractionsCallback),
450 nstglobalcomm_, fplog, cr, mdlog, constr, inputrec, mdAtoms,
451 nrnb, wcycle, fr, vsite, imdSession, pull_work);
452 neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(domDecHelper_.get()));
455 const bool simulationsShareResetCounters = false;
456 resetHandler_ = std::make_unique<ResetHandler>(
457 compat::make_not_null<SimulationSignal*>(&signals_[eglsRESETCOUNTERS]),
458 simulationsShareResetCounters, inputrec->nsteps, MASTER(cr),
459 mdrunOptions.timingOptions.resetHalfway, mdrunOptions.maximumHoursToRun,
460 mdlog, wcycle, walltime_accounting);
463 * Build signaller list
465 * Note that as signallers depend on each others, the order of calling the signallers
466 * matters. It is the responsibility of this builder to ensure that the order is
469 auto energySignaller = energySignallerBuilder.build(
470 inputrec->nstcalcenergy);
471 trajectoryElementBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
472 loggingSignallerBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
473 auto trajectoryElement = trajectoryElementBuilder.build(
474 fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier,
475 inputrec, top_global, oenv, wcycle, startingBehavior);
477 // Add checkpoint helper here since we need a pointer to the trajectory element and
478 // need to register it with the lastStepSignallerBuilder
479 auto checkpointHandler = std::make_unique<CheckpointHandler>(
480 compat::make_not_null<SimulationSignal*>(&signals_[eglsCHKPT]),
481 simulationsShareState, inputrec->nstlist == 0, MASTER(cr),
482 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
483 std::vector<ICheckpointHelperClient*> checkpointClients = {statePropagatorDataPtr, energyElementPtr};
484 checkpointHelper_ = std::make_unique<CheckpointHelper>(
485 std::move(checkpointClients),
486 std::move(checkpointHandler),
487 inputrec->init_step, trajectoryElement.get(),
488 top_global->natoms, fplog, cr,
489 observablesHistory, walltime_accounting, state_global,
490 mdrunOptions.writeConfout);
491 lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(checkpointHelper_.get()));
493 lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(trajectoryElement.get()));
494 auto loggingSignaller = loggingSignallerBuilder.build(
498 lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(loggingSignaller.get()));
499 auto lastStepSignaller = lastStepSignallerBuilder.build(
503 neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(lastStepSignaller.get()));
504 auto neighborSearchSignaller = neighborSearchSignallerBuilder.build(
509 addToCallListAndMove(std::move(neighborSearchSignaller), signallerCallList_, signallersOwnershipList_);
510 addToCallListAndMove(std::move(lastStepSignaller), signallerCallList_, signallersOwnershipList_);
511 addToCallListAndMove(std::move(loggingSignaller), signallerCallList_, signallersOwnershipList_);
512 addToCallList(trajectoryElement, signallerCallList_);
513 addToCallListAndMove(std::move(energySignaller), signallerCallList_, signallersOwnershipList_);
516 * Build the element list
518 * This is the actual sequence of (non-infrastructure) elements to be run.
519 * For NVE, the trajectory element is used outside of the integrator
520 * (composite) element, as well as the checkpoint helper. The checkpoint
521 * helper should be on top of the loop, and is only part of the simulator
522 * call list to be able to react to the last step being signalled.
524 addToCallList(checkpointHelper_, elementCallList_);
525 addToCallListAndMove(std::move(integrator), elementCallList_, elementsOwnershipList_);
526 addToCallListAndMove(std::move(trajectoryElement), elementCallList_, elementsOwnershipList_);
527 // for vv, we need to setup statePropagatorData after the compute
528 // globals so that we reset the right velocities
529 // TODO: Avoid this by getting rid of the need of resetting velocities in vv
530 elementsOwnershipList_.emplace_back(std::move(statePropagatorData));
531 elementsOwnershipList_.emplace_back(std::move(energyElement));
534 std::unique_ptr<ISimulatorElement> ModularSimulator::buildForces(
535 SignallerBuilder<NeighborSearchSignaller> *neighborSearchSignallerBuilder,
536 SignallerBuilder<EnergySignaller> *energySignallerBuilder,
537 StatePropagatorData *statePropagatorDataPtr,
538 EnergyElement *energyElementPtr)
540 const bool isVerbose = mdrunOptions.verbose;
541 const bool isDynamicBox = inputrecDynamicBox(inputrec);
542 // Check for polarizable models and flexible constraints
543 if (ShellFCElement::doShellsOrFlexConstraints(
544 &topologyHolder_->globalTopology(), constr ? constr->numFlexibleConstraints() : 0))
546 auto shellFCElement = std::make_unique<ShellFCElement>(
547 statePropagatorDataPtr, energyElementPtr, isVerbose, isDynamicBox, fplog,
548 cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle, mdScheduleWork,
549 vsite, imdSession, pull_work, constr, &topologyHolder_->globalTopology());
550 topologyHolder_->registerClient(shellFCElement.get());
551 neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(shellFCElement.get()));
552 energySignallerBuilder->registerSignallerClient(compat::make_not_null(shellFCElement.get()));
554 // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
555 return std::move(shellFCElement);
559 auto forceElement = std::make_unique<ForceElement>(
560 statePropagatorDataPtr, energyElementPtr, isDynamicBox, fplog,
561 cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle,
562 mdScheduleWork, vsite, imdSession, pull_work);
563 topologyHolder_->registerClient(forceElement.get());
564 neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
565 energySignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
567 // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
568 return std::move(forceElement);
572 std::unique_ptr<ISimulatorElement> ModularSimulator::buildIntegrator(
573 SignallerBuilder<NeighborSearchSignaller> *neighborSearchSignallerBuilder,
574 SignallerBuilder<EnergySignaller> *energySignallerBuilder,
575 SignallerBuilder<LoggingSignaller> *loggingSignallerBuilder,
576 TrajectoryElementBuilder *trajectoryElementBuilder,
577 CheckBondedInteractionsCallbackPtr *checkBondedInteractionsCallback,
578 compat::not_null<StatePropagatorData*> statePropagatorDataPtr,
579 compat::not_null<EnergyElement*> energyElementPtr,
580 bool hasReadEkinState)
582 auto forceElement = buildForces(
583 neighborSearchSignallerBuilder,
584 energySignallerBuilder,
585 statePropagatorDataPtr,
588 // list of elements owned by the simulator composite object
589 std::vector< std::unique_ptr<ISimulatorElement> > elementsOwnershipList;
590 // call list of the simulator composite object
591 std::vector< compat::not_null<ISimulatorElement*> > elementCallList;
593 std::function<void()> needToCheckNumberOfBondedInteractions;
594 if (inputrec->eI == eiMD)
596 auto computeGlobalsElement =
597 std::make_unique< ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog> >(
598 statePropagatorDataPtr, energyElementPtr, nstglobalcomm_, fplog, mdlog, cr,
599 inputrec, mdAtoms, nrnb, wcycle, fr,
600 &topologyHolder_->globalTopology(), constr, hasReadEkinState);
601 topologyHolder_->registerClient(computeGlobalsElement.get());
602 energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
603 trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
605 *checkBondedInteractionsCallback = computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
607 auto propagator = std::make_unique< Propagator<IntegrationStep::LeapFrog> >(
608 inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
610 addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
611 addToCallList(statePropagatorDataPtr, elementCallList); // we have a full microstate at time t here!
612 addToCallListAndMove(std::move(propagator), elementCallList, elementsOwnershipList);
615 auto constraintElement = std::make_unique< ConstraintsElement<ConstraintVariable::Positions> >(
616 constr, statePropagatorDataPtr, energyElementPtr, MASTER(cr),
617 fplog, inputrec, mdAtoms->mdatoms());
618 auto constraintElementPtr = compat::make_not_null(constraintElement.get());
619 energySignallerBuilder->registerSignallerClient(constraintElementPtr);
620 trajectoryElementBuilder->registerSignallerClient(constraintElementPtr);
621 loggingSignallerBuilder->registerSignallerClient(constraintElementPtr);
623 addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
626 addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
627 addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
629 else if (inputrec->eI == eiVV)
631 auto computeGlobalsElementAtFullTimeStep =
632 std::make_unique< ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerletAtFullTimeStep> >(
633 statePropagatorDataPtr, energyElementPtr, nstglobalcomm_, fplog, mdlog, cr,
634 inputrec, mdAtoms, nrnb, wcycle, fr,
635 &topologyHolder_->globalTopology(), constr, hasReadEkinState);
636 topologyHolder_->registerClient(computeGlobalsElementAtFullTimeStep.get());
637 energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElementAtFullTimeStep.get()));
638 trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElementAtFullTimeStep.get()));
640 auto computeGlobalsElementAfterCoordinateUpdate =
641 std::make_unique<ComputeGlobalsElement <ComputeGlobalsAlgorithm::VelocityVerletAfterCoordinateUpdate> >(
642 statePropagatorDataPtr, energyElementPtr, nstglobalcomm_, fplog, mdlog, cr,
643 inputrec, mdAtoms, nrnb, wcycle, fr,
644 &topologyHolder_->globalTopology(), constr, hasReadEkinState);
645 topologyHolder_->registerClient(computeGlobalsElementAfterCoordinateUpdate.get());
646 energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElementAfterCoordinateUpdate.get()));
647 trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElementAfterCoordinateUpdate.get()));
649 *checkBondedInteractionsCallback = computeGlobalsElementAfterCoordinateUpdate->getCheckNumberOfBondedInteractionsCallback();
651 auto propagatorVelocities = std::make_unique< Propagator <IntegrationStep::VelocitiesOnly> >(
652 inputrec->delta_t * 0.5, statePropagatorDataPtr, mdAtoms, wcycle);
653 auto propagatorVelocitiesAndPositions = std::make_unique< Propagator <IntegrationStep::VelocityVerletPositionsAndVelocities> >(
654 inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
656 addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
657 addToCallListAndMove(std::move(propagatorVelocities), elementCallList, elementsOwnershipList);
660 auto constraintElement = std::make_unique< ConstraintsElement<ConstraintVariable::Velocities> >(
661 constr, statePropagatorDataPtr, energyElementPtr, MASTER(cr),
662 fplog, inputrec, mdAtoms->mdatoms());
663 energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
664 trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
665 loggingSignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
667 addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
669 addToCallListAndMove(std::move(computeGlobalsElementAtFullTimeStep), elementCallList, elementsOwnershipList);
670 addToCallList(statePropagatorDataPtr, elementCallList); // we have a full microstate at time t here!
671 addToCallListAndMove(std::move(propagatorVelocitiesAndPositions), elementCallList, elementsOwnershipList);
674 auto constraintElement = std::make_unique< ConstraintsElement<ConstraintVariable::Positions> >(
675 constr, statePropagatorDataPtr, energyElementPtr, MASTER(cr),
676 fplog, inputrec, mdAtoms->mdatoms());
677 energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
678 trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
679 loggingSignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
681 addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
683 addToCallListAndMove(std::move(computeGlobalsElementAfterCoordinateUpdate), elementCallList, elementsOwnershipList);
684 addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
688 gmx_fatal(FARGS, "Integrator not implemented for the modular simulator.");
691 auto integrator = std::make_unique<CompositeSimulatorElement>(
692 std::move(elementCallList), std::move(elementsOwnershipList));
693 // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
694 return std::move(integrator);
697 bool ModularSimulator::isInputCompatible(
699 const t_inputrec *inputrec,
701 const gmx_vsite_t *vsite,
702 const gmx_multisim_t *ms,
703 const ReplicaExchangeParameters &replExParams,
707 ObservablesHistory *observablesHistory,
708 const gmx_membed_t *membed)
710 auto conditionalAssert =
711 [exitOnFailure](bool condition, const char* message)
715 GMX_RELEASE_ASSERT(condition, message);
720 bool isInputCompatible = true;
722 // GMX_USE_MODULAR_SIMULATOR allows to use modular simulator also for non-standard uses,
723 // such as the leap-frog integrator
724 const auto modularSimulatorExplicitlyTurnedOn =
725 (getenv("GMX_USE_MODULAR_SIMULATOR") != nullptr);
727 isInputCompatible = isInputCompatible && conditionalAssert(
728 inputrec->eI == eiMD || inputrec->eI == eiVV,
729 "Only integrators md and md-vv are supported by the modular simulator.");
730 isInputCompatible = isInputCompatible && conditionalAssert(
731 inputrec->eI != eiMD || modularSimulatorExplicitlyTurnedOn,
732 "Set GMX_USE_MODULAR_SIMULATOR=ON to use the modular simulator with integrator md.");
733 isInputCompatible = isInputCompatible && conditionalAssert(
735 "Rerun is not supported by the modular simulator.");
736 isInputCompatible = isInputCompatible && conditionalAssert(
737 inputrec->etc == etcNO,
738 "Temperature coupling is not supported by the modular simulator.");
739 isInputCompatible = isInputCompatible && conditionalAssert(
740 inputrec->epc == epcNO,
741 "Pressure coupling is not supported by the modular simulator.");
742 isInputCompatible = isInputCompatible && conditionalAssert(
743 !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec) || inputrecNvtTrotter(inputrec)),
744 "Legacy Trotter decomposition is not supported by the modular simulator.");
745 isInputCompatible = isInputCompatible && conditionalAssert(
746 inputrec->efep == efepNO,
747 "Free energy calculation is not supported by the modular simulator.");
748 isInputCompatible = isInputCompatible && conditionalAssert(
750 "Virtual sites are not supported by the modular simulator.");
751 isInputCompatible = isInputCompatible && conditionalAssert(
753 "AWH is not supported by the modular simulator.");
754 isInputCompatible = isInputCompatible && conditionalAssert(
756 "Multi-sim are not supported by the modular simulator.");
757 isInputCompatible = isInputCompatible && conditionalAssert(
758 replExParams.exchangeInterval == 0,
759 "Replica exchange is not supported by the modular simulator.");
760 isInputCompatible = isInputCompatible && conditionalAssert(
761 fcd->disres.nsystems <= 1,
762 "Ensemble restraints are not supported by the modular simulator.");
763 isInputCompatible = isInputCompatible && conditionalAssert(
764 !doSimulatedAnnealing(inputrec),
765 "Simulated annealing is not supported by the modular simulator.");
766 isInputCompatible = isInputCompatible && conditionalAssert(
768 "Simulated tempering is not supported by the modular simulator.");
769 isInputCompatible = isInputCompatible && conditionalAssert(
770 !inputrec->bExpanded,
771 "Expanded ensemble simulations are not supported by the modular simulator.");
772 isInputCompatible = isInputCompatible && conditionalAssert(
773 !(opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr),
774 "Essential dynamics is not supported by the modular simulator.");
775 isInputCompatible = isInputCompatible && conditionalAssert(
776 inputrec->eSwapCoords == eswapNO,
777 "Ion / water position swapping is not supported by the modular simulator.");
778 isInputCompatible = isInputCompatible && conditionalAssert(
780 "Interactive MD is not supported by the modular simulator.");
781 isInputCompatible = isInputCompatible && conditionalAssert(
783 "Membrane embedding is not supported by the modular simulator.");
784 // TODO: Change this to the boolean passed when we merge the user interface change for the GPU update.
785 isInputCompatible = isInputCompatible && conditionalAssert(
786 getenv("GMX_UPDATE_CONSTRAIN_GPU") == nullptr,
787 "Integration on the GPU is not supported by the modular simulator.");
788 // Modular simulator is centered around NS updates
789 // TODO: think how to handle nstlist == 0
790 isInputCompatible = isInputCompatible && conditionalAssert(
791 inputrec->nstlist != 0,
792 "Simulations without neighbor list update are not supported by the modular simulator.");
793 isInputCompatible = isInputCompatible && conditionalAssert(
795 "GMX_FAHCORE not supported by the modular simulator.");
797 return isInputCompatible;
800 void ModularSimulator::checkInputForDisabledFunctionality()
804 inputrec, doRerun, vsite, ms, replExParams,
805 fcd, nfile, fnm, observablesHistory, membed);
808 SignallerCallbackPtr ModularSimulator::SignalHelper::registerLastStepCallback()
810 return std::make_unique<SignallerCallback>(
811 [this](Step step, Time gmx_unused time){this->lastStep_ = step; });
814 SignallerCallbackPtr ModularSimulator::SignalHelper::registerNSCallback()
816 return std::make_unique<SignallerCallback>(
817 [this](Step step, Time gmx_unused time)
818 {this->nextNSStep_ = step; });