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