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 "shellfcelement.h"
  90 #include "signallers.h"
  91 #include "statepropagatordata.h"
  92 #include "trajectoryelement.h"
  93 #include "vrescalethermostat.h"
  94
  95 namespace gmx
  96 {
  97 void ModularSimulator::run()
  98 {
  99     GMX_LOG(mdlog.info).asParagraph().appendText("Using the modular simulator.");
 100     constructElementsAndSignallers();
 101     simulatorSetup();
 102     for (auto& signaller : signallerCallList_)
 103     {
 104         signaller->signallerSetup();
 105     }
 106     if (domDecHelper_)
 107     {
 108         domDecHelper_->setup();
 109     }
 110
 111     for (auto& element : elementsOwnershipList_)
 112     {
 113         element->elementSetup();
 114     }
 115     if (pmeLoadBalanceHelper_)
 116     {
 117         // State must have been initialized so pmeLoadBalanceHelper_ gets a valid box
 118         pmeLoadBalanceHelper_->setup();
 119     }
 120
 121     while (step_ <= signalHelper_->lastStep_)
 122     {
 123         populateTaskQueue();
 124
 125         while (!taskQueue_.empty())
 126         {
 127             auto task = std::move(taskQueue_.front());
 128             taskQueue_.pop();
 129             // run function
 130             (*task)();
 131         }
 132     }
 133
 134     for (auto& element : elementsOwnershipList_)
 135     {
 136         element->elementTeardown();
 137     }
 138     if (pmeLoadBalanceHelper_)
 139     {
 140         pmeLoadBalanceHelper_->teardown();
 141     }
 142     simulatorTeardown();
 143 }
 144
 145 void ModularSimulator::simulatorSetup()
 146 {
 147     if (!mdrunOptions.writeConfout)
 148     {
 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
 152         // interface.
 153         GMX_LOG(mdlog.info)
 154                 .asParagraph()
 155                 .appendText(
 156                         "The -noconfout functionality is deprecated, and "
 157                         "may be removed in a future version.");
 158     }
 159
 160     if (MASTER(cr))
 161     {
 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)
 166         {
 167             timeString = formatString("%8.1f", static_cast<double>(inputrec->init_step + inputrec->nsteps)
 168                                                        * inputrec->delta_t);
 169         }
 170         else
 171         {
 172             timeString = "infinite";
 173         }
 174         if (inputrec->init_step > 0)
 175         {
 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);
 179         }
 180         else
 181         {
 182             fprintf(stderr, "%s steps, %s ps.\n", gmx_step_str(inputrec->nsteps, sbuf),
 183                     timeString.c_str());
 184         }
 185         fprintf(fplog, "\n");
 186     }
 187
 188     walltime_accounting_start_time(walltime_accounting);
 189     wallcycle_start(wcycle, ewcRUN);
 190     print_start(fplog, cr, walltime_accounting, "mdrun");
 191
 192     step_ = inputrec->init_step;
 193 }
 194
 195 void ModularSimulator::preStep(Step step, Time gmx_unused time, bool isNeighborSearchingStep)
 196 {
 197     if (stopHandler_->stoppingAfterCurrentStep(isNeighborSearchingStep) && step != signalHelper_->lastStep_)
 198     {
 199         /*
 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.
 204          */
 205         // clear task queue
 206         std::queue<SimulatorRunFunctionPtr>().swap(taskQueue_);
 207         // rewind step
 208         step_ = step;
 209         return;
 210     }
 211
 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
 216     // duplication.
 217     stophandlerIsNSStep_    = isNeighborSearchingStep;
 218     stophandlerCurrentStep_ = step;
 219     stopHandler_->setSignal();
 220
 221     wallcycle_start(wcycle, ewcSTEP);
 222 }
 223
 224 void ModularSimulator::postStep(Step step, Time gmx_unused time)
 225 {
 226     // Output stuff
 227     if (MASTER(cr))
 228     {
 229         if (do_per_step(step, inputrec->nstlog))
 230         {
 231             if (fflush(fplog) != 0)
 232             {
 233                 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
 234             }
 235         }
 236     }
 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()))
 243     {
 244         print_time(stderr, walltime_accounting, step, inputrec, cr);
 245     }
 246
 247     double cycles = wallcycle_stop(wcycle, ewcSTEP);
 248     if (DOMAINDECOMP(cr) && wcycle)
 249     {
 250         dd_cycles_add(cr->dd, static_cast<float>(cycles), ddCyclStep);
 251     }
 252
 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);
 257 }
 258
 259 void ModularSimulator::simulatorTeardown()
 260 {
 261
 262     // Stop measuring walltime
 263     walltime_accounting_end_time(walltime_accounting);
 264
 265     if (!thisRankHasDuty(cr, DUTY_PME))
 266     {
 267         /* Tell the PME only node to finish */
 268         gmx_pme_send_finish(cr);
 269     }
 270
 271     walltime_accounting_set_nsteps_done(walltime_accounting, step_ - inputrec->init_step);
 272 }
 273
 274 void ModularSimulator::populateTaskQueue()
 275 {
 276     auto registerRunFunction = std::make_unique<RegisterRunFunction>(
 277             [this](SimulatorRunFunctionPtr ptr) { taskQueue_.push(std::move(ptr)); });
 278
 279     Time startTime = inputrec->init_t;
 280     Time timeStep  = inputrec->delta_t;
 281     Time time      = startTime + step_ * timeStep;
 282
 283     // Run an initial call to the signallers
 284     for (auto& signaller : signallerCallList_)
 285     {
 286         signaller->signal(step_, time);
 287     }
 288
 289     if (checkpointHelper_)
 290     {
 291         checkpointHelper_->run(step_, time);
 292     }
 293
 294     if (pmeLoadBalanceHelper_)
 295     {
 296         pmeLoadBalanceHelper_->run(step_, time);
 297     }
 298     if (domDecHelper_)
 299     {
 300         domDecHelper_->run(step_, time);
 301     }
 302
 303     do
 304     {
 305         // local variables for lambda capturing
 306         const int  step     = step_;
 307         const bool isNSStep = step == signalHelper_->nextNSStep_;
 308
 309         // register pre-step
 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_)
 314         {
 315             element->scheduleTask(step_, time, registerRunFunction);
 316         }
 317         // register post-step
 318         (*registerRunFunction)(
 319                 std::make_unique<SimulatorRunFunction>([this, step, time]() { postStep(step, time); }));
 320
 321         // prepare next step
 322         step_++;
 323         time = startTime + step_ * timeStep;
 324         for (auto& signaller : signallerCallList_)
 325         {
 326             signaller->signal(step_, time);
 327         }
 328     } while (step_ != signalHelper_->nextNSStep_ && step_ <= signalHelper_->lastStep_);
 329 }
 330
 331 void ModularSimulator::constructElementsAndSignallers()
 332 {
 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.
 339      *
 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
 344      *      reduced.
 345      *
 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.
 349      */
 350     bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
 351     if (PAR(cr))
 352     {
 353         gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr);
 354     }
 355     if (hasReadEkinState)
 356     {
 357         restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 358     }
 359
 360     /*
 361      * Build data structures
 362      */
 363     topologyHolder_ =
 364             std::make_unique<TopologyHolder>(*top_global, cr, inputrec, fr, mdAtoms, constr, vsite);
 365
 366     std::unique_ptr<FreeEnergyPerturbationElement> freeEnergyPerturbationElement    = nullptr;
 367     FreeEnergyPerturbationElement*                 freeEnergyPerturbationElementPtr = nullptr;
 368     if (inputrec->efep != efepNO)
 369     {
 370         freeEnergyPerturbationElement =
 371                 std::make_unique<FreeEnergyPerturbationElement>(fplog, inputrec, mdAtoms);
 372         freeEnergyPerturbationElementPtr = freeEnergyPerturbationElement.get();
 373     }
 374
 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());
 381
 382     auto energyElement = std::make_unique<EnergyElement>(
 383             statePropagatorDataPtr, freeEnergyPerturbationElementPtr, top_global, inputrec, mdAtoms,
 384             enerd, ekind, constr, fplog, fcd, mdModulesNotifier, MASTER(cr), observablesHistory,
 385             startingBehavior);
 386     auto energyElementPtr = compat::make_not_null(energyElement.get());
 387
 388     /*
 389      * Build stop handler
 390      */
 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);
 397
 398     /*
 399      * Create simulator builders
 400      */
 401     SignallerBuilder<NeighborSearchSignaller> neighborSearchSignallerBuilder;
 402     SignallerBuilder<LastStepSignaller>       lastStepSignallerBuilder;
 403     SignallerBuilder<LoggingSignaller>        loggingSignallerBuilder;
 404     SignallerBuilder<EnergySignaller>         energySignallerBuilder;
 405     TrajectoryElementBuilder                  trajectoryElementBuilder;
 406
 407     /*
 408      * Register data structures to signallers
 409      */
 410     trajectoryElementBuilder.registerWriterClient(statePropagatorDataPtr);
 411     trajectoryElementBuilder.registerSignallerClient(statePropagatorDataPtr);
 412     lastStepSignallerBuilder.registerSignallerClient(statePropagatorDataPtr);
 413
 414     trajectoryElementBuilder.registerWriterClient(energyElementPtr);
 415     trajectoryElementBuilder.registerSignallerClient(energyElementPtr);
 416     energySignallerBuilder.registerSignallerClient(energyElementPtr);
 417
 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()));
 421
 422     /*
 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.
 426      */
 427     // TODO: Make a CheckpointHelperBuilder
 428     std::vector<ICheckpointHelperClient*> checkpointClients = { statePropagatorDataPtr, energyElementPtr,
 429                                                                 freeEnergyPerturbationElementPtr };
 430     CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback = nullptr;
 431     auto                               integrator =
 432             buildIntegrator(&neighborSearchSignallerBuilder, &energySignallerBuilder,
 433                             &loggingSignallerBuilder, &trajectoryElementBuilder, &checkpointClients,
 434                             &checkBondedInteractionsCallback, statePropagatorDataPtr,
 435                             energyElementPtr, freeEnergyPerturbationElementPtr, hasReadEkinState);
 436
 437     /*
 438      * Build infrastructure elements
 439      */
 440
 441     if (PmeLoadBalanceHelper::doPmeLoadBalancing(mdrunOptions, inputrec, fr))
 442     {
 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()));
 447     }
 448
 449     if (DOMAINDECOMP(cr))
 450     {
 451         GMX_ASSERT(checkBondedInteractionsCallback,
 452                    "Domain decomposition needs a callback for check the number of bonded "
 453                    "interactions.");
 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()));
 459     }
 460
 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);
 467
 468     /*
 469      * Build signaller list
 470      *
 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
 473      * maintained.
 474      */
 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()));
 483
 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()));
 495
 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);
 505
 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_);
 511
 512     /*
 513      * Build the element list
 514      *
 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.
 520      */
 521     addToCallList(checkpointHelper_, elementCallList_);
 522     if (freeEnergyPerturbationElement)
 523     {
 524         addToCallListAndMove(std::move(freeEnergyPerturbationElement), elementCallList_,
 525                              elementsOwnershipList_);
 526     }
 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));
 534 }
 535
 536 std::unique_ptr<ISimulatorElement>
 537 ModularSimulator::buildForces(SignallerBuilder<NeighborSearchSignaller>* neighborSearchSignallerBuilder,
 538                               SignallerBuilder<EnergySignaller>*         energySignallerBuilder,
 539                               StatePropagatorData*                       statePropagatorDataPtr,
 540                               EnergyElement*                             energyElementPtr,
 541                               FreeEnergyPerturbationElement* freeEnergyPerturbationElement)
 542 {
 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))
 548     {
 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()));
 557
 558         // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 559         return std::move(shellFCElement);
 560     }
 561     else
 562     {
 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()));
 570
 571         // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 572         return std::move(forceElement);
 573     }
 574 }
 575
 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)
 587 {
 588     auto forceElement =
 589             buildForces(neighborSearchSignallerBuilder, energySignallerBuilder,
 590                         statePropagatorDataPtr, energyElementPtr, freeEnergyPerturbationElementPtr);
 591
 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;
 596
 597     std::function<void()> needToCheckNumberOfBondedInteractions;
 598     if (inputrec->eI == eiMD)
 599     {
 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()));
 609
 610         *checkBondedInteractionsCallback =
 611                 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
 612
 613         auto propagator = std::make_unique<Propagator<IntegrationStep::LeapFrog>>(
 614                 inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
 615
 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)
 619         {
 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);
 631         }
 632
 633         std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
 634         if (inputrec->epc == epcPARRINELLORAHMAN)
 635         {
 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());
 645         }
 646         addToCallListAndMove(std::move(propagator), elementCallList, elementsOwnershipList);
 647         if (constr)
 648         {
 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);
 656
 657             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 658         }
 659
 660         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
 661         addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
 662         if (prBarostat)
 663         {
 664             addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
 665         }
 666     }
 667     else if (inputrec->eI == eiVV)
 668     {
 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()));
 678
 679         *checkBondedInteractionsCallback =
 680                 computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
 681
 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);
 687
 688         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
 689
 690         std::unique_ptr<ParrinelloRahmanBarostat> prBarostat = nullptr;
 691         if (inputrec->epc == epcPARRINELLORAHMAN)
 692         {
 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());
 702         }
 703         addToCallListAndMove(std::move(propagatorVelocities), elementCallList, elementsOwnershipList);
 704         if (constr)
 705         {
 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()));
 714
 715             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 716         }
 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)
 720         {
 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);
 734         }
 735         addToCallListAndMove(std::move(propagatorVelocitiesAndPositions), elementCallList,
 736                              elementsOwnershipList);
 737         if (constr)
 738         {
 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()));
 747
 748             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 749         }
 750         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
 751         addToCallList(energyElementPtr, elementCallList); // we have the energies at time t here!
 752         if (prBarostat)
 753         {
 754             addToCallListAndMove(std::move(prBarostat), elementCallList, elementsOwnershipList);
 755         }
 756     }
 757     else
 758     {
 759         gmx_fatal(FARGS, "Integrator not implemented for the modular simulator.");
 760     }
 761
 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);
 766 }
 767
 768 bool ModularSimulator::isInputCompatible(bool                             exitOnFailure,
 769                                          const t_inputrec*                inputrec,
 770                                          bool                             doRerun,
 771                                          const gmx_mtop_t&                globalTopology,
 772                                          const gmx_multisim_t*            ms,
 773                                          const ReplicaExchangeParameters& replExParams,
 774                                          const t_fcdata*                  fcd,
 775                                          bool                             doEssentialDynamics,
 776                                          bool                             doMembed)
 777 {
 778     auto conditionalAssert = [exitOnFailure](bool condition, const char* message) {
 779         if (exitOnFailure)
 780         {
 781             GMX_RELEASE_ASSERT(condition, message);
 782         }
 783         return condition;
 784     };
 785
 786     bool isInputCompatible = true;
 787
 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);
 794
 795     GMX_RELEASE_ASSERT(
 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 "
 799             "use, "
 800             "or unset both to recover default behavior.");
 801
 802     GMX_RELEASE_ASSERT(
 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.");
 809
 810     isInputCompatible =
 811             isInputCompatible
 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.");
 819     isInputCompatible =
 820             isInputCompatible
 821             && conditionalAssert(!doRerun, "Rerun is not supported by the modular simulator.");
 822     isInputCompatible =
 823             isInputCompatible
 824             && conditionalAssert(
 825                        inputrec->etc == etcNO || inputrec->etc == etcVRESCALE,
 826                        "Only v-rescale thermostat is supported by the modular simulator.");
 827     isInputCompatible =
 828             isInputCompatible
 829             && conditionalAssert(
 830                        inputrec->epc == epcNO || inputrec->epc == epcPARRINELLORAHMAN,
 831                        "Only Parrinello-Rahman barostat is supported by the modular simulator.");
 832     isInputCompatible =
 833             isInputCompatible
 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.");
 846     isInputCompatible =
 847             isInputCompatible
 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.");
 852     isInputCompatible =
 853             isInputCompatible
 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.");
 858     isInputCompatible =
 859             isInputCompatible
 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.");
 867     isInputCompatible =
 868             isInputCompatible
 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.");
 874     isInputCompatible =
 875             isInputCompatible
 876             && conditionalAssert(gmx_mtop_ftype_count(globalTopology, F_DISRES) == 0,
 877                                  "Distance restraints are not supported by the modular simulator.");
 878     isInputCompatible =
 879             isInputCompatible
 880             && conditionalAssert(
 881                        gmx_mtop_ftype_count(globalTopology, F_ORIRES) == 0,
 882                        "Orientation restraints are not supported by the modular simulator.");
 883     isInputCompatible =
 884             isInputCompatible
 885             && conditionalAssert(ms == nullptr,
 886                                  "Multi-sim are not supported by the modular simulator.");
 887     isInputCompatible =
 888             isInputCompatible
 889             && conditionalAssert(replExParams.exchangeInterval == 0,
 890                                  "Replica exchange is not supported by the modular simulator.");
 891
 892     int numEnsembleRestraintSystems;
 893     if (fcd)
 894     {
 895         numEnsembleRestraintSystems = fcd->disres.nsystems;
 896     }
 897     else
 898     {
 899         auto distantRestraintEnsembleEnvVar = getenv("GMX_DISRE_ENSEMBLE_SIZE");
 900         numEnsembleRestraintSystems =
 901                 (ms != nullptr && distantRestraintEnsembleEnvVar != nullptr)
 902                         ? static_cast<int>(strtol(distantRestraintEnsembleEnvVar, nullptr, 10))
 903                         : 0;
 904     }
 905     isInputCompatible =
 906             isInputCompatible
 907             && conditionalAssert(numEnsembleRestraintSystems <= 1,
 908                                  "Ensemble restraints are not supported by the modular simulator.");
 909     isInputCompatible =
 910             isInputCompatible
 911             && conditionalAssert(!doSimulatedAnnealing(inputrec),
 912                                  "Simulated annealing is not supported by the modular simulator.");
 913     isInputCompatible =
 914             isInputCompatible
 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.");
 921     isInputCompatible =
 922             isInputCompatible
 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.");
 929     isInputCompatible =
 930             isInputCompatible
 931             && conditionalAssert(!inputrec->bIMD,
 932                                  "Interactive MD is not supported by the modular simulator.");
 933     isInputCompatible =
 934             isInputCompatible
 935             && conditionalAssert(!doMembed,
 936                                  "Membrane embedding is not supported by the modular simulator.");
 937     const bool useGraph = !areMoleculesDistributedOverPbc(*inputrec, globalTopology, MDLogger());
 938     isInputCompatible =
 939             isInputCompatible
 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.
 942     isInputCompatible =
 943             isInputCompatible
 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.");
 956
 957     return isInputCompatible;
 958 }
 959
 960 void ModularSimulator::checkInputForDisabledFunctionality()
 961 {
 962     isInputCompatible(true, inputrec, doRerun, *top_global, ms, replExParams, fcd,
 963                       opt2bSet("-ei", nfile, fnm), membed != nullptr);
 964     if (observablesHistory->edsamHistory)
 965     {
 966         gmx_fatal(FARGS,
 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.");
 970     }
 971 }
 972
 973 SignallerCallbackPtr ModularSimulator::SignalHelper::registerLastStepCallback()
 974 {
 975     return std::make_unique<SignallerCallback>(
 976             [this](Step step, Time gmx_unused time) { this->lastStep_ = step; });
 977 }
 978
 979 SignallerCallbackPtr ModularSimulator::SignalHelper::registerNSCallback()
 980 {
 981     return std::make_unique<SignallerCallback>(
 982             [this](Step step, Time gmx_unused time) { this->nextNSStep_ = step; });
 983 }
 984 } // namespace gmx