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