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