src/gromacs/modularsimulator/modularsimulator.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*! \internal \file
  36  * \brief Defines the modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "modularsimulator.h"
  45
  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"
  80
  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 "signallers.h"
  90 #include "statepropagatordata.h"
  91 #include "trajectoryelement.h"
  92 #include "vrescalethermostat.h"
  93
  94 namespace gmx
  95 {
  96 void ModularSimulator::run()
  97 {
  98     GMX_LOG(mdlog.info).asParagraph().appendText("Using the modular simulator.");
  99     constructElementsAndSignallers();
 100     simulatorSetup();
 101     for (auto& signaller : signallerCallList_)
 102     {
 103         signaller->signallerSetup();
 104     }
 105     if (domDecHelper_)
 106     {
 107         domDecHelper_->setup();
 108     }
 109
 110     for (auto& element : elementsOwnershipList_)
 111     {
 112         element->elementSetup();
 113     }
 114     if (pmeLoadBalanceHelper_)
 115     {
 116         // State must have been initialized so pmeLoadBalanceHelper_ gets a valid box
 117         pmeLoadBalanceHelper_->setup();
 118     }
 119
 120     while (step_ <= signalHelper_->lastStep_)
 121     {
 122         populateTaskQueue();
 123
 124         while (!taskQueue_.empty())
 125         {
 126             auto task = std::move(taskQueue_.front());
 127             taskQueue_.pop();
 128             // run function
 129             (*task)();
 130         }
 131     }
 132
 133     for (auto& element : elementsOwnershipList_)
 134     {
 135         element->elementTeardown();
 136     }
 137     if (pmeLoadBalanceHelper_)
 138     {
 139         pmeLoadBalanceHelper_->teardown();
 140     }
 141     simulatorTeardown();
 142 }
 143
 144 void ModularSimulator::simulatorSetup()
 145 {
 146     if (!mdrunOptions.writeConfout)
 147     {
 148         // This is on by default, and the main known use case for
 149         // turning it off is for convenience in benchmarking, which is
 150         // something that should not show up in the general user
 151         // interface.
 152         GMX_LOG(mdlog.info)
 153                 .asParagraph()
 154                 .appendText(
 155                         "The -noconfout functionality is deprecated, and "
 156                         "may be removed in a future version.");
 157     }
 158
 159     if (MASTER(cr))
 160     {
 161         char        sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
 162         std::string timeString;
 163         fprintf(stderr, "starting mdrun '%s'\n", *(top_global->name));
 164         if (inputrec->nsteps >= 0)
 165         {
 166             timeString = formatString("%8.1f", static_cast<double>(inputrec->init_step + inputrec->nsteps)
 167                                                        * inputrec->delta_t);
 168         }
 169         else
 170         {
 171             timeString = "infinite";
 172         }
 173         if (inputrec->init_step > 0)
 174         {
 175             fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 176                     gmx_step_str(inputrec->init_step + inputrec->nsteps, sbuf), timeString.c_str(),
 177                     gmx_step_str(inputrec->init_step, sbuf2), inputrec->init_step * inputrec->delta_t);
 178         }
 179         else
 180         {
 181             fprintf(stderr, "%s steps, %s ps.\n", gmx_step_str(inputrec->nsteps, sbuf),
 182                     timeString.c_str());
 183         }
 184         fprintf(fplog, "\n");
 185     }
 186
 187     walltime_accounting_start_time(walltime_accounting);
 188     wallcycle_start(wcycle, ewcRUN);
 189     print_start(fplog, cr, walltime_accounting, "mdrun");
 190
 191     step_ = inputrec->init_step;
 192 }
 193
 194 void ModularSimulator::preStep(Step step, Time gmx_unused time, bool isNeighborSearchingStep)
 195 {
 196     if (stopHandler_->stoppingAfterCurrentStep(isNeighborSearchingStep) && step != signalHelper_->lastStep_)
 197     {
 198         /*
 199          * Stop handler wants to stop after the current step, which was
 200          * not known when building the current task queue. This happens
 201          * e.g. when a stop is signalled by OS. We therefore want to purge
 202          * the task queue now, and re-schedule this step as last step.
 203          */
 204         // clear task queue
 205         std::queue<SimulatorRunFunctionPtr>().swap(taskQueue_);
 206         // rewind step
 207         step_ = step;
 208         return;
 209     }
 210
 211     resetHandler_->setSignal(walltime_accounting);
 212     // This is a hack to avoid having to rewrite StopHandler to be a NeighborSearchSignaller
 213     // and accept the step as input. Eventually, we want to do that, but currently this would
 214     // require introducing NeighborSearchSignaller in the legacy do_md or a lot of code
 215     // duplication.
 216     stophandlerIsNSStep_    = isNeighborSearchingStep;
 217     stophandlerCurrentStep_ = step;
 218     stopHandler_->setSignal();
 219
 220     wallcycle_start(wcycle, ewcSTEP);
 221 }
 222
 223 void ModularSimulator::postStep(Step step, Time gmx_unused time)
 224 {
 225     // Output stuff
 226     if (MASTER(cr))
 227     {
 228         if (do_per_step(step, inputrec->nstlog))
 229         {
 230             if (fflush(fplog) != 0)
 231             {
 232                 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
 233             }
 234         }
 235     }
 236     const bool do_verbose = mdrunOptions.verbose
 237                             && (step % mdrunOptions.verboseStepPrintInterval == 0
 238                                 || step == inputrec->init_step || step == signalHelper_->lastStep_);
 239     // Print the remaining wall clock time for the run
 240     if (MASTER(cr) && (do_verbose || gmx_got_usr_signal())
 241         && !(pmeLoadBalanceHelper_ && pmeLoadBalanceHelper_->pmePrinting()))
 242     {
 243         print_time(stderr, walltime_accounting, step, inputrec, cr);
 244     }
 245
 246     double cycles = wallcycle_stop(wcycle, ewcSTEP);
 247     if (DOMAINDECOMP(cr) && wcycle)
 248     {
 249         dd_cycles_add(cr->dd, static_cast<float>(cycles), ddCyclStep);
 250     }
 251
 252     resetHandler_->resetCounters(
 253             step, step - inputrec->init_step, mdlog, fplog, cr, fr->nbv.get(), nrnb, fr->pmedata,
 254             pmeLoadBalanceHelper_ ? pmeLoadBalanceHelper_->loadBalancingObject() : nullptr, wcycle,
 255             walltime_accounting);
 256 }
 257
 258 void ModularSimulator::simulatorTeardown()
 259 {
 260
 261     // Stop measuring walltime
 262     walltime_accounting_end_time(walltime_accounting);
 263
 264     if (!thisRankHasDuty(cr, DUTY_PME))
 265     {
 266         /* Tell the PME only node to finish */
 267         gmx_pme_send_finish(cr);
 268     }
 269
 270     walltime_accounting_set_nsteps_done(walltime_accounting, step_ - inputrec->init_step);
 271 }
 272
 273 void ModularSimulator::populateTaskQueue()
 274 {
 275     auto registerRunFunction = std::make_unique<RegisterRunFunction>(
 276             [this](SimulatorRunFunctionPtr ptr) { taskQueue_.push(std::move(ptr)); });
 277
 278     Time startTime = inputrec->init_t;
 279     Time timeStep  = inputrec->delta_t;
 280     Time time      = startTime + step_ * timeStep;
 281
 282     // Run an initial call to the signallers
 283     for (auto& signaller : signallerCallList_)
 284     {
 285         signaller->signal(step_, time);
 286     }
 287
 288     if (checkpointHelper_)
 289     {
 290         checkpointHelper_->run(step_, time);
 291     }
 292
 293     if (pmeLoadBalanceHelper_)
 294     {
 295         pmeLoadBalanceHelper_->run(step_, time);
 296     }
 297     if (domDecHelper_)
 298     {
 299         domDecHelper_->run(step_, time);
 300     }
 301
 302     do
 303     {
 304         // local variables for lambda capturing
 305         const int  step     = step_;
 306         const bool isNSStep = step == signalHelper_->nextNSStep_;
 307
 308         // register pre-step
 309         (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
 310                 [this, step, time, isNSStep]() { preStep(step, time, isNSStep); }));
 311         // register elements for step
 312         for (auto& element : elementCallList_)
 313         {
 314             element->scheduleTask(step_, time, registerRunFunction);
 315         }
 316         // register post-step
 317         (*registerRunFunction)(
 318                 std::make_unique<SimulatorRunFunction>([this, step, time]() { postStep(step, time); }));
 319
 320         // prepare next step
 321         step_++;
 322         time = startTime + step_ * timeStep;
 323         for (auto& signaller : signallerCallList_)
 324         {
 325             signaller->signal(step_, time);
 326         }
 327     } while (step_ != signalHelper_->nextNSStep_ && step_ <= signalHelper_->lastStep_);
 328 }
 329
 330 void ModularSimulator::constructElementsAndSignallers()
 331 {
 332     /* When restarting from a checkpoint, it can be appropriate to
 333      * initialize ekind from quantities in the checkpoint. Otherwise,
 334      * compute_globals must initialize ekind before the simulation
 335      * starts/restarts. However, only the master rank knows what was
 336      * found in the checkpoint file, so we have to communicate in
 337      * order to coordinate the restart.
 338      *
 339      * TODO (modular) This should become obsolete when checkpoint reading
 340      *      happens within the modular simulator framework: The energy
 341      *      element should read its data from the checkpoint file pointer,
 342      *      and signal to the compute globals element if it needs anything
 343      *      reduced.
 344      *
 345      * TODO (legacy) Consider removing this communication if/when checkpoint
 346      *      reading directly follows .tpr reading, because all ranks can
 347      *      agree on hasReadEkinState at that time.
 348      */
 349     bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
 350     if (PAR(cr))
 351     {
 352         gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr->mpi_comm_mygroup);
 353     }
 354     if (hasReadEkinState)
 355     {
 356         restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 357     }
 358
 359     /*
 360      * Build data structures
 361      */
 362     topologyHolder_ =
 363             std::make_unique<TopologyHolder>(*top_global, cr, inputrec, fr, mdAtoms, constr, vsite);
 364
 365     std::unique_ptr<FreeEnergyPerturbationElement> freeEnergyPerturbationElement    = nullptr;
 366     FreeEnergyPerturbationElement*                 freeEnergyPerturbationElementPtr = nullptr;
 367     if (inputrec->efep != efepNO)
 368     {
 369         freeEnergyPerturbationElement =
 370                 std::make_unique<FreeEnergyPerturbationElement>(fplog, inputrec, mdAtoms);
 371         freeEnergyPerturbationElementPtr = freeEnergyPerturbationElement.get();
 372     }
 373
 374     auto statePropagatorData = std::make_unique<StatePropagatorData>(
 375             top_global->natoms, fplog, cr, state_global, inputrec->nstxout, inputrec->nstvout,
 376             inputrec->nstfout, inputrec->nstxout_compressed, fr->nbv->useGpu(),
 377             freeEnergyPerturbationElementPtr, topologyHolder_.get(), fr->bMolPBC,
 378             mdrunOptions.writeConfout, opt2fn("-c", nfile, fnm), inputrec, mdAtoms->mdatoms());
 379     auto statePropagatorDataPtr = compat::make_not_null(statePropagatorData.get());
 380
 381     auto energyElement = std::make_unique<EnergyElement>(
 382             statePropagatorDataPtr, freeEnergyPerturbationElementPtr, top_global, inputrec, mdAtoms,
 383             enerd, ekind, constr, fplog, fcd, mdModulesNotifier, MASTER(cr), observablesHistory,
 384             startingBehavior);
 385     auto energyElementPtr = compat::make_not_null(energyElement.get());
 386
 387     /*
 388      * Build stop handler
 389      */
 390     const bool simulationsShareState = false;
 391     stopHandler_                     = stopHandlerBuilder->getStopHandlerMD(
 392             compat::not_null<SimulationSignal*>(&signals_[eglsSTOPCOND]), simulationsShareState,
 393             MASTER(cr), inputrec->nstlist, mdrunOptions.reproducible, nstglobalcomm_,
 394             mdrunOptions.maximumHoursToRun, inputrec->nstlist == 0, fplog, stophandlerCurrentStep_,
 395             stophandlerIsNSStep_, walltime_accounting);
 396
 397     /*
 398      * Create simulator builders
 399      */
 400     SignallerBuilder<NeighborSearchSignaller> neighborSearchSignallerBuilder;
 401     SignallerBuilder<LastStepSignaller>       lastStepSignallerBuilder;
 402     SignallerBuilder<LoggingSignaller>        loggingSignallerBuilder;
 403     SignallerBuilder<EnergySignaller>         energySignallerBuilder;
 404     TrajectoryElementBuilder                  trajectoryElementBuilder;
 405
 406     /*
 407      * Register data structures to signallers
 408      */
 409     trajectoryElementBuilder.registerWriterClient(statePropagatorDataPtr);
 410     trajectoryElementBuilder.registerSignallerClient(statePropagatorDataPtr);
 411     lastStepSignallerBuilder.registerSignallerClient(statePropagatorDataPtr);
 412
 413     trajectoryElementBuilder.registerWriterClient(energyElementPtr);
 414     trajectoryElementBuilder.registerSignallerClient(energyElementPtr);
 415     energySignallerBuilder.registerSignallerClient(energyElementPtr);
 416
 417     // Register the simulator itself to the neighbor search / last step signaller
 418     neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
 419     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
 420
 421     /*
 422      * Build integrator - this takes care of force calculation, propagation,
 423      * constraining, and of the place the statePropagatorData and the energy element
 424      * have a full timestep state.
 425      */
 426     // TODO: Make a CheckpointHelperBuilder
 427     std::vector<ICheckpointHelperClient*> checkpointClients = { statePropagatorDataPtr, energyElementPtr,
 428                                                                 freeEnergyPerturbationElementPtr };
 429     CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback = nullptr;
 430     auto                               integrator =
 431             buildIntegrator(&neighborSearchSignallerBuilder, &energySignallerBuilder,
 432                             &loggingSignallerBuilder, &trajectoryElementBuilder, &checkpointClients,
 433                             &checkBondedInteractionsCallback, statePropagatorDataPtr,
 434                             energyElementPtr, freeEnergyPerturbationElementPtr, hasReadEkinState);
 435
 436     /*
 437      * Build infrastructure elements
 438      */
 439
 440     if (PmeLoadBalanceHelper::doPmeLoadBalancing(mdrunOptions, inputrec, fr))
 441     {
 442         pmeLoadBalanceHelper_ = std::make_unique<PmeLoadBalanceHelper>(
 443                 mdrunOptions.verbose, statePropagatorDataPtr, fplog, cr, mdlog, inputrec, wcycle, fr);
 444         neighborSearchSignallerBuilder.registerSignallerClient(
 445                 compat::make_not_null(pmeLoadBalanceHelper_.get()));
 446     }
 447
 448     if (DOMAINDECOMP(cr))
 449     {
 450         GMX_ASSERT(checkBondedInteractionsCallback,
 451                    "Domain decomposition needs a callback for check the number of bonded "
 452                    "interactions.");
 453         domDecHelper_ = std::make_unique<DomDecHelper>(
 454                 mdrunOptions.verbose, mdrunOptions.verboseStepPrintInterval, statePropagatorDataPtr,
 455                 topologyHolder_.get(), std::move(checkBondedInteractionsCallback), nstglobalcomm_, fplog,
 456                 cr, mdlog, constr, inputrec, mdAtoms, nrnb, wcycle, fr, vsite, imdSession, pull_work);
 457         neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(domDecHelper_.get()));
 458     }
 459
 460     const bool simulationsShareResetCounters = false;
 461     resetHandler_                            = std::make_unique<ResetHandler>(
 462             compat::make_not_null<SimulationSignal*>(&signals_[eglsRESETCOUNTERS]),
 463             simulationsShareResetCounters, inputrec->nsteps, MASTER(cr),
 464             mdrunOptions.timingOptions.resetHalfway, mdrunOptions.maximumHoursToRun, mdlog, wcycle,
 465             walltime_accounting);
 466
 467     /*
 468      * Build signaller list
 469      *
 470      * Note that as signallers depend on each others, the order of calling the signallers
 471      * matters. It is the responsibility of this builder to ensure that the order is
 472      * maintained.
 473      */
 474     auto energySignaller = energySignallerBuilder.build(
 475             inputrec->nstcalcenergy, inputrec->fepvals->nstdhdl, inputrec->nstpcouple);
 476     trajectoryElementBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
 477     loggingSignallerBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
 478     auto trajectoryElement = trajectoryElementBuilder.build(
 479             fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, inputrec,
 480             top_global, oenv, wcycle, startingBehavior, simulationsShareState);
 481     loggingSignallerBuilder.registerSignallerClient(compat::make_not_null(trajectoryElement.get()));
 482
 483     // Add checkpoint helper here since we need a pointer to the trajectory element and
 484     // need to register it with the lastStepSignallerBuilder
 485     auto checkpointHandler = std::make_unique<CheckpointHandler>(
 486             compat::make_not_null<SimulationSignal*>(&signals_[eglsCHKPT]), simulationsShareState,
 487             inputrec->nstlist == 0, MASTER(cr), mdrunOptions.writeConfout,
 488             mdrunOptions.checkpointOptions.period);
 489     checkpointHelper_ = std::make_unique<CheckpointHelper>(
 490             std::move(checkpointClients), std::move(checkpointHandler), inputrec->init_step,
 491             trajectoryElement.get(), top_global->natoms, fplog, cr, observablesHistory,
 492             walltime_accounting, state_global, mdrunOptions.writeConfout);
 493     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(checkpointHelper_.get()));
 494
 495     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(trajectoryElement.get()));
 496     auto loggingSignaller =
 497             loggingSignallerBuilder.build(inputrec->nstlog, inputrec->init_step, inputrec->init_t);
 498     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(loggingSignaller.get()));
 499     auto lastStepSignaller =
 500             lastStepSignallerBuilder.build(inputrec->nsteps, inputrec->init_step, stopHandler_.get());
 501     neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(lastStepSignaller.get()));
 502     auto neighborSearchSignaller = neighborSearchSignallerBuilder.build(
 503             inputrec->nstlist, inputrec->init_step, inputrec->init_t);
 504
 505     addToCallListAndMove(std::move(neighborSearchSignaller), signallerCallList_, signallersOwnershipList_);
 506     addToCallListAndMove(std::move(lastStepSignaller), signallerCallList_, signallersOwnershipList_);
 507     addToCallListAndMove(std::move(loggingSignaller), signallerCallList_, signallersOwnershipList_);
 508     addToCallList(trajectoryElement, signallerCallList_);
 509     addToCallListAndMove(std::move(energySignaller), signallerCallList_, signallersOwnershipList_);
 510
 511     /*
 512      * Build the element list
 513      *
 514      * This is the actual sequence of (non-infrastructure) elements to be run.
 515      * For NVE, the trajectory element is used outside of the integrator
 516      * (composite) element, as well as the checkpoint helper. The checkpoint
 517      * helper should be on top of the loop, and is only part of the simulator
 518      * call list to be able to react to the last step being signalled.
 519      */
 520     addToCallList(checkpointHelper_, elementCallList_);
 521     if (freeEnergyPerturbationElement)
 522     {
 523         addToCallListAndMove(std::move(freeEnergyPerturbationElement), elementCallList_,
 524                              elementsOwnershipList_);
 525     }
 526     addToCallListAndMove(std::move(integrator), elementCallList_, elementsOwnershipList_);
 527     addToCallListAndMove(std::move(trajectoryElement), elementCallList_, elementsOwnershipList_);
 528     // for vv, we need to setup statePropagatorData after the compute
 529     // globals so that we reset the right velocities
 530     // TODO: Avoid this by getting rid of the need of resetting velocities in vv
 531     elementsOwnershipList_.emplace_back(std::move(statePropagatorData));
 532     elementsOwnershipList_.emplace_back(std::move(energyElement));
 533 }
 534
 535 std::unique_ptr<ISimulatorElement>
 536 ModularSimulator::buildForces(SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
 537                               SignallerBuilder<EnergySignaller>*         energySignallerBuilder,
 538                               StatePropagatorData*                       statePropagatorDataPtr,
 539                               EnergyElement*                             energyElementPtr,
 540                               FreeEnergyPerturbationElement* freeEnergyPerturbationElement)
 541 {
 542     const bool isVerbose    = mdrunOptions.verbose;
 543     const bool isDynamicBox = inputrecDynamicBox(inputrec);
 544
 545     auto forceElement = std::make_unique<ForceElement>(
 546             statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElement, isVerbose,
 547             isDynamicBox, fplog, cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle, runScheduleWork, vsite,
 548             imdSession, pull_work, constr, &topologyHolder_->globalTopology(), enforcedRotation);
 549     topologyHolder_->registerClient(forceElement.get());
 550     neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
 551     energySignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
 552
 553     // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 554     return std::move(forceElement);
 555 }
 556
 557 std::unique_ptr<ISimulatorElement> ModularSimulator::buildIntegrator(
 558         SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
 559         SignallerBuilder<EnergySignaller>*         energySignallerBuilder,
 560         SignallerBuilder<LoggingSignaller>*        loggingSignallerBuilder,
 561         TrajectoryElementBuilder*                  trajectoryElementBuilder,
 562         std::vector<ICheckpointHelperClient*>*     checkpointClients,
 563         CheckBondedInteractionsCallbackPtr*        checkBondedInteractionsCallback,
 564         compat::not_null<StatePropagatorData*>     statePropagatorDataPtr,
 565         compat::not_null<EnergyElement*>           energyElementPtr,
 566         FreeEnergyPerturbationElement*             freeEnergyPerturbationElementPtr,
 567         bool                                       hasReadEkinState)
 568 {
 569     auto forceElement =
 570             buildForces(neighborSearchSignallerBuilder, energySignallerBuilder,
 571                         statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr);
 572
 573     // list of elements owned by the simulator composite object
 574     std::vector<std::unique_ptr<ISimulatorElement>> elementsOwnershipList;
 575     // call list of the simulator composite object
 576     std::vector<compat::not_null<ISimulatorElement*>> elementCallList;
 577
 578     std::function<void()> needToCheckNumberOfBondedInteractions;
 579     if (inputrec->eI == eiMD)
 580     {
 581         auto computeGlobalsElement =
 582                 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog>>(
 583                         statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
 584                         &signals_, nstglobalcomm_, fplog, mdlog, cr, inputrec, mdAtoms, nrnb,
 585                         wcycle, fr, &topologyHolder_->globalTopology(), constr, hasReadEkinState);
 586         topologyHolder_->registerClient(computeGlobalsElement.get());
 587         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
 588         trajectoryElementBuilder->registerSignallerClient(
 589                 compat::make_not_null(computeGlobalsElement.get()));
 590
 591         *checkBondedInteractionsCallback =
 592                 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
 593
 594         auto propagator = std::make_unique<Propagator<IntegrationStep::LeapFrog>>(
 595                 inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
 596
 597         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
 598         addToCallList(statePropagatorDataPtr, elementCallList); // we have a full microstate at time t here!
 599         if (inputrec->etc == etcVRESCALE)
 600         {
 601             // TODO: With increased complexity of the propagator, this will need further development,
 602             //       e.g. using propagators templated for velocity propagation policies and a builder
 603             propagator->setNumVelocityScalingVariables(inputrec->opts.ngtc);
 604             auto thermostat = std::make_unique<VRescaleThermostat>(
 605                     inputrec->nsttcouple, -1, false, inputrec->ld_seed, inputrec->opts.ngtc,
 606                     inputrec->delta_t * inputrec->nsttcouple, inputrec->opts.ref_t, inputrec->opts.tau_t,
 607                     inputrec->opts.nrdf, energyElementPtr, propagator->viewOnVelocityScaling(),
 608                     propagator->velocityScalingCallback(), state_global, cr, inputrec->bContinuation);
 609             checkpointClients->emplace_back(thermostat.get());
 610             energyElementPtr->setVRescaleThermostat(thermostat.get());
 611             addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
 612         }
 613
 614         std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
 615         if (inputrec->epc == epcPARRINELLORAHMAN)
 616         {
 617             // Building the PR barostat here since it needs access to the propagator
 618             // and we want to be able to move the propagator object
 619             prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
 620                     inputrec->nstpcouple, -1, inputrec->delta_t * inputrec->nstpcouple,
 621                     inputrec->init_step, propagator->viewOnPRScalingMatrix(),
 622                     propagator->prScalingCallback(), statePropagatorDataPtr, energyElementPtr,
 623                     fplog, inputrec, mdAtoms, state_global, cr, inputrec->bContinuation);
 624             energyElementPtr->setParrinelloRahamnBarostat(prBarostat.get());
 625             checkpointClients->emplace_back(prBarostat.get());
 626         }
 627         addToCallListAndMove(std::move(propagator), elementCallList, elementsOwnershipList);
 628         if (constr)
 629         {
 630             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
 631                     constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
 632                     MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
 633             auto constraintElementPtr = compat::make_not_null(constraintElement.get());
 634             energySignallerBuilder->registerSignallerClient(constraintElementPtr);
 635             trajectoryElementBuilder->registerSignallerClient(constraintElementPtr);
 636             loggingSignallerBuilder->registerSignallerClient(constraintElementPtr);
 637
 638             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 639         }
 640
 641         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
 642         addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
 643         if (prBarostat)
 644         {
 645             addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
 646         }
 647     }
 648     else if (inputrec->eI == eiVV)
 649     {
 650         auto computeGlobalsElement =
 651                 std::make_unique<ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerlet>>(
 652                         statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
 653                         &signals_, nstglobalcomm_, fplog, mdlog, cr, inputrec, mdAtoms, nrnb,
 654                         wcycle, fr, &topologyHolder_->globalTopology(), constr, hasReadEkinState);
 655         topologyHolder_->registerClient(computeGlobalsElement.get());
 656         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
 657         trajectoryElementBuilder->registerSignallerClient(
 658                 compat::make_not_null(computeGlobalsElement.get()));
 659
 660         *checkBondedInteractionsCallback =
 661                 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
 662
 663         auto propagatorVelocities = std::make_unique<Propagator<IntegrationStep::VelocitiesOnly>>(
 664                 inputrec->delta_t * 0.5, statePropagatorDataPtr, mdAtoms, wcycle);
 665         auto propagatorVelocitiesAndPositions =
 666                 std::make_unique<Propagator<IntegrationStep::VelocityVerletPositionsAndVelocities>>(
 667                         inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
 668
 669         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
 670
 671         std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
 672         if (inputrec->epc == epcPARRINELLORAHMAN)
 673         {
 674             // Building the PR barostat here since it needs access to the propagator
 675             // and we want to be able to move the propagator object
 676             prBarostat = std::make_unique<ParrinelloRahmanBarostat>(
 677                     inputrec->nstpcouple, -1, inputrec->delta_t * inputrec->nstpcouple,
 678                     inputrec->init_step, propagatorVelocities->viewOnPRScalingMatrix(),
 679                     propagatorVelocities->prScalingCallback(), statePropagatorDataPtr, energyElementPtr,
 680                     fplog, inputrec, mdAtoms, state_global, cr, inputrec->bContinuation);
 681             energyElementPtr->setParrinelloRahamnBarostat(prBarostat.get());
 682             checkpointClients->emplace_back(prBarostat.get());
 683         }
 684         addToCallListAndMove(std::move(propagatorVelocities), elementCallList, elementsOwnershipList);
 685         if (constr)
 686         {
 687             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Velocities>>(
 688                     constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
 689                     MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
 690             energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 691             trajectoryElementBuilder->registerSignallerClient(
 692                     compat::make_not_null(constraintElement.get()));
 693             loggingSignallerBuilder->registerSignallerClient(
 694                     compat::make_not_null(constraintElement.get()));
 695
 696             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 697         }
 698         addToCallList(compat::make_not_null(computeGlobalsElement.get()), elementCallList);
 699         addToCallList(statePropagatorDataPtr, elementCallList); // we have a full microstate at time t here!
 700         if (inputrec->etc == etcVRESCALE)
 701         {
 702             // TODO: With increased complexity of the propagator, this will need further development,
 703             //       e.g. using propagators templated for velocity propagation policies and a builder
 704             propagatorVelocitiesAndPositions->setNumVelocityScalingVariables(inputrec->opts.ngtc);
 705             auto thermostat = std::make_unique<VRescaleThermostat>(
 706                     inputrec->nsttcouple, 0, true, inputrec->ld_seed, inputrec->opts.ngtc,
 707                     inputrec->delta_t * inputrec->nsttcouple, inputrec->opts.ref_t,
 708                     inputrec->opts.tau_t, inputrec->opts.nrdf, energyElementPtr,
 709                     propagatorVelocitiesAndPositions->viewOnVelocityScaling(),
 710                     propagatorVelocitiesAndPositions->velocityScalingCallback(), state_global, cr,
 711                     inputrec->bContinuation);
 712             checkpointClients->emplace_back(thermostat.get());
 713             energyElementPtr->setVRescaleThermostat(thermostat.get());
 714             addToCallListAndMove(std::move(thermostat), elementCallList, elementsOwnershipList);
 715         }
 716         addToCallListAndMove(std::move(propagatorVelocitiesAndPositions), elementCallList,
 717                              elementsOwnershipList);
 718         if (constr)
 719         {
 720             auto constraintElement = std::make_unique<ConstraintsElement<ConstraintVariable::Positions>>(
 721                     constr, statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr,
 722                     MASTER(cr), fplog, inputrec, mdAtoms->mdatoms());
 723             energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 724             trajectoryElementBuilder->registerSignallerClient(
 725                     compat::make_not_null(constraintElement.get()));
 726             loggingSignallerBuilder->registerSignallerClient(
 727                     compat::make_not_null(constraintElement.get()));
 728
 729             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 730         }
 731         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
 732         addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
 733         if (prBarostat)
 734         {
 735             addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
 736         }
 737     }
 738     else
 739     {
 740         gmx_fatal(FARGS, "Integrator not implemented for the modular simulator.");
 741     }
 742
 743     auto integrator = std::make_unique<CompositeSimulatorElement>(std::move(elementCallList),
 744                                                                   std::move(elementsOwnershipList));
 745     // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 746     return std::move(integrator);
 747 }
 748
 749 bool ModularSimulator::isInputCompatible(bool                             exitOnFailure,
 750                                          const t_inputrec*                inputrec,
 751                                          bool                             doRerun,
 752                                          const gmx_mtop_t&                globalTopology,
 753                                          const gmx_multisim_t*            ms,
 754                                          const ReplicaExchangeParameters& replExParams,
 755                                          const t_fcdata*                  fcd,
 756                                          bool                             doEssentialDynamics,
 757                                          bool                             doMembed)
 758 {
 759     auto conditionalAssert = [exitOnFailure](bool condition, const char* message) {
 760         if (exitOnFailure)
 761         {
 762             GMX_RELEASE_ASSERT(condition, message);
 763         }
 764         return condition;
 765     };
 766
 767     bool isInputCompatible = true;
 768
 769     // GMX_USE_MODULAR_SIMULATOR allows to use modular simulator also for non-standard uses,
 770     // such as the leap-frog integrator
 771     const auto modularSimulatorExplicitlyTurnedOn = (getenv("GMX_USE_MODULAR_SIMULATOR") != nullptr);
 772     // GMX_USE_MODULAR_SIMULATOR allows to use disable modular simulator for all uses,
 773     // including the velocity-verlet integrator used by default
 774     const auto modularSimulatorExplicitlyTurnedOff = (getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
 775
 776     GMX_RELEASE_ASSERT(
 777             !(modularSimulatorExplicitlyTurnedOn && modularSimulatorExplicitlyTurnedOff),
 778             "Cannot have both GMX_USE_MODULAR_SIMULATOR=ON and GMX_DISABLE_MODULAR_SIMULATOR=ON. "
 779             "Unset one of the two environment variables to explicitly chose which simulator to "
 780             "use, "
 781             "or unset both to recover default behavior.");
 782
 783     GMX_RELEASE_ASSERT(
 784             !(modularSimulatorExplicitlyTurnedOff && inputrec->eI == eiVV
 785               && inputrec->epc == epcPARRINELLORAHMAN),
 786             "Cannot use a Parrinello-Rahman barostat with md-vv and "
 787             "GMX_DISABLE_MODULAR_SIMULATOR=ON, "
 788             "as the Parrinello-Rahman barostat is not implemented in the legacy simulator. Unset "
 789             "GMX_DISABLE_MODULAR_SIMULATOR or use a different pressure control algorithm.");
 790
 791     isInputCompatible =
 792             isInputCompatible
 793             && conditionalAssert(
 794                        inputrec->eI == eiMD || inputrec->eI == eiVV,
 795                        "Only integrators md and md-vv are supported by the modular simulator.");
 796     isInputCompatible = isInputCompatible
 797                         && conditionalAssert(inputrec->eI != eiMD || modularSimulatorExplicitlyTurnedOn,
 798                                              "Set GMX_USE_MODULAR_SIMULATOR=ON to use the modular "
 799                                              "simulator with integrator md.");
 800     isInputCompatible =
 801             isInputCompatible
 802             && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
 803     isInputCompatible =
 804             isInputCompatible
 805             && conditionalAssert(
 806                        inputrec->etc == etcNO || inputrec->etc == etcVRESCALE,
 807                        "Only v-rescale thermostat is supported by the modular simulator.");
 808     isInputCompatible =
 809             isInputCompatible
 810             && conditionalAssert(
 811                        inputrec->epc == epcNO || inputrec->epc == epcPARRINELLORAHMAN,
 812                        "Only Parrinello-Rahman barostat is supported by the modular simulator.");
 813     isInputCompatible =
 814             isInputCompatible
 815             && conditionalAssert(
 816                        !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec)
 817                          || inputrecNvtTrotter(inputrec)),
 818                        "Legacy Trotter decomposition is not supported by the modular simulator.");
 819     isInputCompatible = isInputCompatible
 820                         && conditionalAssert(inputrec->efep == efepNO || inputrec->efep == efepYES
 821                                                      || inputrec->efep == efepSLOWGROWTH,
 822                                              "Expanded ensemble free energy calculation is not "
 823                                              "supported by the modular simulator.");
 824     isInputCompatible = isInputCompatible
 825                         && conditionalAssert(!inputrec->bPull,
 826                                              "Pulling is not supported by the modular simulator.");
 827     isInputCompatible =
 828             isInputCompatible
 829             && conditionalAssert(inputrec->opts.ngacc == 1 && inputrec->opts.acc[0][XX] == 0.0
 830                                          && inputrec->opts.acc[0][YY] == 0.0
 831                                          && inputrec->opts.acc[0][ZZ] == 0.0 && inputrec->cos_accel == 0.0,
 832                                  "Acceleration is not supported by the modular simulator.");
 833     isInputCompatible =
 834             isInputCompatible
 835             && conditionalAssert(inputrec->opts.ngfrz == 1 && inputrec->opts.nFreeze[0][XX] == 0
 836                                          && inputrec->opts.nFreeze[0][YY] == 0
 837                                          && inputrec->opts.nFreeze[0][ZZ] == 0,
 838                                  "Freeze groups are not supported by the modular simulator.");
 839     isInputCompatible =
 840             isInputCompatible
 841             && conditionalAssert(
 842                        inputrec->deform[XX][XX] == 0.0 && inputrec->deform[XX][YY] == 0.0
 843                                && inputrec->deform[XX][ZZ] == 0.0 && inputrec->deform[YY][XX] == 0.0
 844                                && inputrec->deform[YY][YY] == 0.0 && inputrec->deform[YY][ZZ] == 0.0
 845                                && inputrec->deform[ZZ][XX] == 0.0 && inputrec->deform[ZZ][YY] == 0.0
 846                                && inputrec->deform[ZZ][ZZ] == 0.0,
 847                        "Deformation is not supported by the modular simulator.");
 848     isInputCompatible =
 849             isInputCompatible
 850             && conditionalAssert(gmx_mtop_interaction_count(globalTopology, IF_VSITE) == 0,
 851                                  "Virtual sites are not supported by the modular simulator.");
 852     isInputCompatible = isInputCompatible
 853                         && conditionalAssert(!inputrec->bDoAwh,
 854                                              "AWH is not supported by the modular simulator.");
 855     isInputCompatible =
 856             isInputCompatible
 857             && conditionalAssert(gmx_mtop_ftype_count(globalTopology, F_DISRES) == 0,
 858                                  "Distance restraints are not supported by the modular simulator.");
 859     isInputCompatible =
 860             isInputCompatible
 861             && conditionalAssert(
 862                        gmx_mtop_ftype_count(globalTopology, F_ORIRES) == 0,
 863                        "Orientation restraints are not supported by the modular simulator.");
 864     isInputCompatible =
 865             isInputCompatible
 866             && conditionalAssert(ms == nullptr,
 867                                  "Multi-sim are not supported by the modular simulator.");
 868     isInputCompatible =
 869             isInputCompatible
 870             && conditionalAssert(replExParams.exchangeInterval == 0,
 871                                  "Replica exchange is not supported by the modular simulator.");
 872
 873     int numEnsembleRestraintSystems;
 874     if (fcd)
 875     {
 876         numEnsembleRestraintSystems = fcd->disres.nsystems;
 877     }
 878     else
 879     {
 880         auto distantRestraintEnsembleEnvVar = getenv("GMX_DISRE_ENSEMBLE_SIZE");
 881         numEnsembleRestraintSystems =
 882                 (ms != nullptr && distantRestraintEnsembleEnvVar != nullptr)
 883                         ? static_cast<int>(strtol(distantRestraintEnsembleEnvVar, nullptr, 10))
 884                         : 0;
 885     }
 886     isInputCompatible =
 887             isInputCompatible
 888             && conditionalAssert(numEnsembleRestraintSystems <= 1,
 889                                  "Ensemble restraints are not supported by the modular simulator.");
 890     isInputCompatible =
 891             isInputCompatible
 892             && conditionalAssert(!doSimulatedAnnealing(inputrec),
 893                                  "Simulated annealing is not supported by the modular simulator.");
 894     isInputCompatible =
 895             isInputCompatible
 896             && conditionalAssert(!inputrec->bSimTemp,
 897                                  "Simulated tempering is not supported by the modular simulator.");
 898     isInputCompatible = isInputCompatible
 899                         && conditionalAssert(!inputrec->bExpanded,
 900                                              "Expanded ensemble simulations are not supported by "
 901                                              "the modular simulator.");
 902     isInputCompatible =
 903             isInputCompatible
 904             && conditionalAssert(!doEssentialDynamics,
 905                                  "Essential dynamics is not supported by the modular simulator.");
 906     isInputCompatible = isInputCompatible
 907                         && conditionalAssert(inputrec->eSwapCoords == eswapNO,
 908                                              "Ion / water position swapping is not supported by "
 909                                              "the modular simulator.");
 910     isInputCompatible =
 911             isInputCompatible
 912             && conditionalAssert(!inputrec->bIMD,
 913                                  "Interactive MD is not supported by the modular simulator.");
 914     isInputCompatible =
 915             isInputCompatible
 916             && conditionalAssert(!doMembed,
 917                                  "Membrane embedding is not supported by the modular simulator.");
 918     // TODO: Change this to the boolean passed when we merge the user interface change for the GPU update.
 919     isInputCompatible =
 920             isInputCompatible
 921             && conditionalAssert(
 922                        getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") == nullptr,
 923                        "Integration on the GPU is not supported by the modular simulator.");
 924     // Modular simulator is centered around NS updates
 925     // TODO: think how to handle nstlist == 0
 926     isInputCompatible = isInputCompatible
 927                         && conditionalAssert(inputrec->nstlist != 0,
 928                                              "Simulations without neighbor list update are not "
 929                                              "supported by the modular simulator.");
 930     isInputCompatible = isInputCompatible
 931                         && conditionalAssert(!GMX_FAHCORE,
 932                                              "GMX_FAHCORE not supported by the modular simulator.");
 933
 934     return isInputCompatible;
 935 }
 936
 937 void ModularSimulator::checkInputForDisabledFunctionality()
 938 {
 939     isInputCompatible(true, inputrec, doRerun, *top_global, ms, replExParams, fcd,
 940                       opt2bSet("-ei", nfile, fnm), membed != nullptr);
 941     if (observablesHistory->edsamHistory)
 942     {
 943         gmx_fatal(FARGS,
 944                   "The checkpoint is from a run with essential dynamics sampling, "
 945                   "but the current run did not specify the -ei option. "
 946                   "Either specify the -ei option to mdrun, or do not use this checkpoint file.");
 947     }
 948 }
 949
 950 SignallerCallbackPtr ModularSimulator::SignalHelper::registerLastStepCallback()
 951 {
 952     return std::make_unique<SignallerCallback>(
 953             [this](Step step, Time gmx_unused time) { this->lastStep_ = step; });
 954 }
 955
 956 SignallerCallbackPtr ModularSimulator::SignalHelper::registerNSCallback()
 957 {
 958     return std::make_unique<SignallerCallback>(
 959             [this](Step step, Time gmx_unused time) { this->nextNSStep_ = step; });
 960 }
 961 } // namespace gmx