src/gromacs/modularsimulator/modularsimulator.cpp

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