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
  36  * \brief Defines the modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_mdrun
  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/mdlib/checkpointhandler.h"
  53 #include "gromacs/mdlib/constr.h"
  54 #include "gromacs/mdlib/energyoutput.h"
  55 #include "gromacs/mdlib/mdatoms.h"
  56 #include "gromacs/mdlib/resethandler.h"
  57 #include "gromacs/mdlib/stat.h"
  58 #include "gromacs/mdlib/update.h"
  59 #include "gromacs/mdrun/replicaexchange.h"
  60 #include "gromacs/mdrun/shellfc.h"
  61 #include "gromacs/mdrunutility/handlerestart.h"
  62 #include "gromacs/mdrunutility/printtime.h"
  63 #include "gromacs/mdtypes/commrec.h"
  64 #include "gromacs/mdtypes/fcdata.h"
  65 #include "gromacs/mdtypes/inputrec.h"
  66 #include "gromacs/mdtypes/mdrunoptions.h"
  67 #include "gromacs/mdtypes/observableshistory.h"
  68 #include "gromacs/mdtypes/state.h"
  69 #include "gromacs/nbnxm/nbnxm.h"
  70 #include "gromacs/timing/walltime_accounting.h"
  71 #include "gromacs/topology/topology.h"
  72 #include "gromacs/utility/cstringutil.h"
  73 #include "gromacs/utility/fatalerror.h"
  74
  75 #include "compositesimulatorelement.h"
  76 #include "computeglobalselement.h"
  77 #include "constraintelement.h"
  78 #include "energyelement.h"
  79 #include "forceelement.h"
  80 #include "propagator.h"
  81 #include "shellfcelement.h"
  82 #include "signallers.h"
  83 #include "statepropagatordata.h"
  84 #include "trajectoryelement.h"
  85
  86 namespace gmx
  87 {
  88 void ModularSimulator::run()
  89 {
  90     GMX_LOG(mdlog.info).asParagraph().
  91         appendText("Using the modular simulator.");
  92     constructElementsAndSignallers();
  93     simulatorSetup();
  94     for (auto &signaller : signallerCallList_)
  95     {
  96         signaller->signallerSetup();
  97     }
  98     if (pmeLoadBalanceHelper_)
  99     {
 100         pmeLoadBalanceHelper_->setup();
 101     }
 102     if (domDecHelper_)
 103     {
 104         domDecHelper_->setup();
 105     }
 106
 107     for (auto &element : elementsOwnershipList_)
 108     {
 109         element->elementSetup();
 110     }
 111
 112     while (step_ <= signalHelper_->lastStep_)
 113     {
 114         populateTaskQueue();
 115
 116         while (!taskQueue_.empty())
 117         {
 118             auto task = std::move(taskQueue_.front());
 119             taskQueue_.pop();
 120             // run function
 121             (*task)();
 122         }
 123     }
 124
 125     for (auto &element : elementsOwnershipList_)
 126     {
 127         element->elementTeardown();
 128     }
 129     if (pmeLoadBalanceHelper_)
 130     {
 131         pmeLoadBalanceHelper_->teardown();
 132     }
 133     simulatorTeardown();
 134 }
 135
 136 void ModularSimulator::simulatorSetup()
 137 {
 138     if (!mdrunOptions.writeConfout)
 139     {
 140         // This is on by default, and the main known use case for
 141         // turning it off is for convenience in benchmarking, which is
 142         // something that should not show up in the general user
 143         // interface.
 144         GMX_LOG(mdlog.info).asParagraph().
 145             appendText("The -noconfout functionality is deprecated, and "
 146                        "may be removed in a future version.");
 147     }
 148
 149     if (MASTER(cr))
 150     {
 151         char        sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
 152         std::string timeString;
 153         fprintf(stderr, "starting mdrun '%s'\n",
 154                 *(top_global->name));
 155         if (inputrec->nsteps >= 0)
 156         {
 157             timeString = formatString(
 158                         "%8.1f", static_cast<double>(inputrec->init_step+inputrec->nsteps)*inputrec->delta_t);
 159         }
 160         else
 161         {
 162             timeString = "infinite";
 163         }
 164         if (inputrec->init_step > 0)
 165         {
 166             fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 167                     gmx_step_str(inputrec->init_step+inputrec->nsteps, sbuf),
 168                     timeString.c_str(),
 169                     gmx_step_str(inputrec->init_step, sbuf2),
 170                     inputrec->init_step*inputrec->delta_t);
 171         }
 172         else
 173         {
 174             fprintf(stderr, "%s steps, %s ps.\n",
 175                     gmx_step_str(inputrec->nsteps, sbuf), timeString.c_str());
 176         }
 177         fprintf(fplog, "\n");
 178     }
 179
 180     walltime_accounting_start_time(walltime_accounting);
 181     wallcycle_start(wcycle, ewcRUN);
 182     print_start(fplog, cr, walltime_accounting, "mdrun");
 183
 184     step_ = inputrec->init_step;
 185 }
 186
 187 void ModularSimulator::preStep(
 188         Step step, Time gmx_unused time,
 189         bool isNeighborSearchingStep)
 190 {
 191     if (stopHandler_->stoppingAfterCurrentStep(isNeighborSearchingStep) &&
 192         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) &&
 237         (do_verbose || gmx_got_usr_signal()) &&
 238         !(pmeLoadBalanceHelper_ && pmeLoadBalanceHelper_->pmePrinting()))
 239     {
 240         print_time(stderr, walltime_accounting, step, inputrec, cr);
 241     }
 242
 243     double cycles = wallcycle_stop(wcycle, ewcSTEP);
 244     if (DOMAINDECOMP(cr) && wcycle)
 245     {
 246         dd_cycles_add(cr->dd, static_cast<float>(cycles), ddCyclStep);
 247     }
 248
 249     resetHandler_->resetCounters(
 250             step, step - inputrec->init_step, mdlog, fplog, cr, fr->nbv.get(),
 251             nrnb, fr->pmedata,
 252             pmeLoadBalanceHelper_ ? pmeLoadBalanceHelper_->loadBalancingObject() : nullptr,
 253             wcycle, walltime_accounting);
 254 }
 255
 256 void ModularSimulator::simulatorTeardown()
 257 {
 258
 259     // Stop measuring walltime
 260     walltime_accounting_end_time(walltime_accounting);
 261
 262     if (!thisRankHasDuty(cr, DUTY_PME))
 263     {
 264         /* Tell the PME only node to finish */
 265         gmx_pme_send_finish(cr);
 266     }
 267
 268     walltime_accounting_set_nsteps_done(walltime_accounting, step_ - inputrec->init_step);
 269 }
 270
 271 void ModularSimulator::populateTaskQueue()
 272 {
 273     auto registerRunFunction = std::make_unique<RegisterRunFunction>(
 274                 [this](SimulatorRunFunctionPtr ptr)
 275                 {taskQueue_.push(std::move(ptr)); });
 276
 277     Time startTime = inputrec->init_t;
 278     Time timeStep  = inputrec->delta_t;
 279     Time time      = startTime + step_*timeStep;
 280
 281     // Run an initial call to the signallers
 282     for (auto &signaller : signallerCallList_)
 283     {
 284         signaller->signal(step_, time);
 285     }
 286
 287     if (checkpointHelper_)
 288     {
 289         checkpointHelper_->run(step_, time);
 290     }
 291
 292     if (pmeLoadBalanceHelper_)
 293     {
 294         pmeLoadBalanceHelper_->run(step_, time);
 295     }
 296     if (domDecHelper_)
 297     {
 298         domDecHelper_->run(step_, time);
 299     }
 300
 301     do
 302     {
 303         // local variables for lambda capturing
 304         const int  step     = step_;
 305         const bool isNSStep = step == signalHelper_->nextNSStep_;
 306
 307         // register pre-step
 308         (*registerRunFunction)(
 309                 std::make_unique<SimulatorRunFunction>(
 310                         [this, step, time, isNSStep](){preStep(step, time, isNSStep); }));
 311         // register elements for step
 312         for (auto &element : elementCallList_)
 313         {
 314             element->scheduleTask(step_, time, registerRunFunction);
 315         }
 316         // register post-step
 317         (*registerRunFunction)(
 318                 std::make_unique<SimulatorRunFunction>(
 319                         [this, step, time](){postStep(step, time); }));
 320
 321         // prepare next step
 322         step_++;
 323         time = startTime + step_*timeStep;
 324         for (auto &signaller : signallerCallList_)
 325         {
 326             signaller->signal(step_, time);
 327         }
 328     }
 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);
 355     }
 356     if (hasReadEkinState)
 357     {
 358         restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 359     }
 360
 361     /*
 362      * Build data structures
 363      */
 364     auto statePropagatorData = std::make_unique<StatePropagatorData>(
 365                 top_global->natoms, fplog, cr, state_global,
 366                 inputrec->nstxout, inputrec->nstvout,
 367                 inputrec->nstfout, inputrec->nstxout_compressed,
 368                 fr->nbv->useGpu(), inputrec, mdAtoms->mdatoms());
 369     auto statePropagatorDataPtr = compat::make_not_null(statePropagatorData.get());
 370
 371     auto energyElement = std::make_unique<EnergyElement>(
 372                 statePropagatorDataPtr, top_global, inputrec, mdAtoms, enerd, ekind,
 373                 constr, fplog, fcd, mdModulesNotifier, MASTER(cr), observablesHistory, startingBehavior);
 374     auto energyElementPtr = compat::make_not_null(energyElement.get());
 375
 376     topologyHolder_ = std::make_unique<TopologyHolder>(
 377                 *top_global, cr, inputrec, fr,
 378                 mdAtoms, constr, vsite);
 379
 380     /*
 381      * Build stop handler
 382      */
 383     const bool simulationsShareState = false;
 384     stopHandler_ = stopHandlerBuilder->getStopHandlerMD(
 385                 compat::not_null<SimulationSignal*>(&signals_[eglsSTOPCOND]),
 386                 simulationsShareState, MASTER(cr), inputrec->nstlist, mdrunOptions.reproducible,
 387                 nstglobalcomm_, mdrunOptions.maximumHoursToRun, inputrec->nstlist == 0, fplog,
 388                 stophandlerCurrentStep_, stophandlerIsNSStep_, walltime_accounting);
 389
 390     /*
 391      * Create simulator builders
 392      */
 393     SignallerBuilder<NeighborSearchSignaller> neighborSearchSignallerBuilder;
 394     SignallerBuilder<LastStepSignaller>       lastStepSignallerBuilder;
 395     SignallerBuilder<LoggingSignaller>        loggingSignallerBuilder;
 396     SignallerBuilder<EnergySignaller>         energySignallerBuilder;
 397     TrajectoryElementBuilder                  trajectoryElementBuilder;
 398
 399     /*
 400      * Register data structures to signallers
 401      */
 402     trajectoryElementBuilder.registerWriterClient(statePropagatorDataPtr);
 403     trajectoryElementBuilder.registerSignallerClient(statePropagatorDataPtr);
 404
 405     trajectoryElementBuilder.registerWriterClient(energyElementPtr);
 406     trajectoryElementBuilder.registerSignallerClient(energyElementPtr);
 407     energySignallerBuilder.registerSignallerClient(energyElementPtr);
 408     loggingSignallerBuilder.registerSignallerClient(energyElementPtr);
 409
 410     // Register the simulator itself to the neighbor search / last step signaller
 411     neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
 412     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(signalHelper_.get()));
 413
 414     /*
 415      * Build integrator - this takes care of force calculation, propagation,
 416      * constraining, and of the place the statePropagatorData and the energy element
 417      * have a full timestep state.
 418      */
 419     CheckBondedInteractionsCallbackPtr checkBondedInteractionsCallback = nullptr;
 420     auto integrator = buildIntegrator(
 421                 &neighborSearchSignallerBuilder,
 422                 &energySignallerBuilder,
 423                 &loggingSignallerBuilder,
 424                 &trajectoryElementBuilder,
 425                 &checkBondedInteractionsCallback,
 426                 statePropagatorDataPtr,
 427                 energyElementPtr,
 428                 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,
 438                     cr, mdlog, inputrec, wcycle, fr);
 439         neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(pmeLoadBalanceHelper_.get()));
 440     }
 441
 442     if (DOMAINDECOMP(cr))
 443     {
 444         GMX_ASSERT(
 445                 checkBondedInteractionsCallback,
 446                 "Domain decomposition needs a callback for check the number of bonded interactions.");
 447         domDecHelper_ = std::make_unique<DomDecHelper>(
 448                     mdrunOptions.verbose, mdrunOptions.verboseStepPrintInterval,
 449                     statePropagatorDataPtr, topologyHolder_.get(), std::move(checkBondedInteractionsCallback),
 450                     nstglobalcomm_, fplog, cr, mdlog, constr, inputrec, mdAtoms,
 451                     nrnb, wcycle, fr, vsite, imdSession, pull_work);
 452         neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(domDecHelper_.get()));
 453     }
 454
 455     const bool simulationsShareResetCounters = false;
 456     resetHandler_ = std::make_unique<ResetHandler>(
 457                 compat::make_not_null<SimulationSignal*>(&signals_[eglsRESETCOUNTERS]),
 458                 simulationsShareResetCounters, inputrec->nsteps, MASTER(cr),
 459                 mdrunOptions.timingOptions.resetHalfway, mdrunOptions.maximumHoursToRun,
 460                 mdlog, wcycle, walltime_accounting);
 461
 462     /*
 463      * Build signaller list
 464      *
 465      * Note that as signallers depend on each others, the order of calling the signallers
 466      * matters. It is the responsibility of this builder to ensure that the order is
 467      * maintained.
 468      */
 469     auto energySignaller = energySignallerBuilder.build(
 470                 inputrec->nstcalcenergy);
 471     trajectoryElementBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
 472     loggingSignallerBuilder.registerSignallerClient(compat::make_not_null(energySignaller.get()));
 473     auto trajectoryElement = trajectoryElementBuilder.build(
 474                 fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier,
 475                 inputrec, top_global, oenv, wcycle, startingBehavior);
 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]),
 481                 simulationsShareState, inputrec->nstlist == 0, MASTER(cr),
 482                 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
 483     std::vector<ICheckpointHelperClient*> checkpointClients = {statePropagatorDataPtr, energyElementPtr};
 484     checkpointHelper_ = std::make_unique<CheckpointHelper>(
 485                 std::move(checkpointClients),
 486                 std::move(checkpointHandler),
 487                 inputrec->init_step, trajectoryElement.get(),
 488                 top_global->natoms, fplog, cr,
 489                 observablesHistory, walltime_accounting, state_global,
 490                 mdrunOptions.writeConfout);
 491     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(checkpointHelper_.get()));
 492
 493     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(trajectoryElement.get()));
 494     auto loggingSignaller = loggingSignallerBuilder.build(
 495                 inputrec->nstlog,
 496                 inputrec->init_step,
 497                 inputrec->init_t);
 498     lastStepSignallerBuilder.registerSignallerClient(compat::make_not_null(loggingSignaller.get()));
 499     auto lastStepSignaller = lastStepSignallerBuilder.build(
 500                 inputrec->nsteps,
 501                 inputrec->init_step,
 502                 stopHandler_.get());
 503     neighborSearchSignallerBuilder.registerSignallerClient(compat::make_not_null(lastStepSignaller.get()));
 504     auto neighborSearchSignaller = neighborSearchSignallerBuilder.build(
 505                 inputrec->nstlist,
 506                 inputrec->init_step,
 507                 inputrec->init_t);
 508
 509     addToCallListAndMove(std::move(neighborSearchSignaller), signallerCallList_, signallersOwnershipList_);
 510     addToCallListAndMove(std::move(lastStepSignaller), signallerCallList_, signallersOwnershipList_);
 511     addToCallListAndMove(std::move(loggingSignaller), signallerCallList_, signallersOwnershipList_);
 512     addToCallList(trajectoryElement, signallerCallList_);
 513     addToCallListAndMove(std::move(energySignaller), signallerCallList_, signallersOwnershipList_);
 514
 515     /*
 516      * Build the element list
 517      *
 518      * This is the actual sequence of (non-infrastructure) elements to be run.
 519      * For NVE, the trajectory element is used outside of the integrator
 520      * (composite) element, as well as the checkpoint helper. The checkpoint
 521      * helper should be on top of the loop, and is only part of the simulator
 522      * call list to be able to react to the last step being signalled.
 523      */
 524     addToCallList(checkpointHelper_, elementCallList_);
 525     addToCallListAndMove(std::move(integrator), elementCallList_, elementsOwnershipList_);
 526     addToCallListAndMove(std::move(trajectoryElement), elementCallList_, elementsOwnershipList_);
 527     // for vv, we need to setup statePropagatorData after the compute
 528     // globals so that we reset the right velocities
 529     // TODO: Avoid this by getting rid of the need of resetting velocities in vv
 530     elementsOwnershipList_.emplace_back(std::move(statePropagatorData));
 531     elementsOwnershipList_.emplace_back(std::move(energyElement));
 532 }
 533
 534 std::unique_ptr<ISimulatorElement> ModularSimulator::buildForces(
 535         SignallerBuilder<NeighborSearchSignaller> *neighborSearchSignallerBuilder,
 536         SignallerBuilder<EnergySignaller>         *energySignallerBuilder,
 537         StatePropagatorData                       *statePropagatorDataPtr,
 538         EnergyElement                             *energyElementPtr)
 539 {
 540     const bool isVerbose    = mdrunOptions.verbose;
 541     const bool isDynamicBox = inputrecDynamicBox(inputrec);
 542     // Check for polarizable models and flexible constraints
 543     if (ShellFCElement::doShellsOrFlexConstraints(
 544                 &topologyHolder_->globalTopology(), constr ? constr->numFlexibleConstraints() : 0))
 545     {
 546         auto shellFCElement = std::make_unique<ShellFCElement>(
 547                     statePropagatorDataPtr, energyElementPtr, isVerbose, isDynamicBox, fplog,
 548                     cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle, runScheduleWork,
 549                     vsite, imdSession, pull_work, constr, &topologyHolder_->globalTopology());
 550         topologyHolder_->registerClient(shellFCElement.get());
 551         neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(shellFCElement.get()));
 552         energySignallerBuilder->registerSignallerClient(compat::make_not_null(shellFCElement.get()));
 553
 554         // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 555         return std::move(shellFCElement);
 556     }
 557     else
 558     {
 559         auto forceElement = std::make_unique<ForceElement>(
 560                     statePropagatorDataPtr, energyElementPtr, isDynamicBox, fplog,
 561                     cr, inputrec, mdAtoms, nrnb, fr, fcd, wcycle,
 562                     runScheduleWork, vsite, imdSession, pull_work);
 563         topologyHolder_->registerClient(forceElement.get());
 564         neighborSearchSignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
 565         energySignallerBuilder->registerSignallerClient(compat::make_not_null(forceElement.get()));
 566
 567         // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 568         return std::move(forceElement);
 569     }
 570 }
 571
 572 std::unique_ptr<ISimulatorElement> ModularSimulator::buildIntegrator(
 573         SignallerBuilder<NeighborSearchSignaller> *neighborSearchSignallerBuilder,
 574         SignallerBuilder<EnergySignaller>         *energySignallerBuilder,
 575         SignallerBuilder<LoggingSignaller>        *loggingSignallerBuilder,
 576         TrajectoryElementBuilder                  *trajectoryElementBuilder,
 577         CheckBondedInteractionsCallbackPtr        *checkBondedInteractionsCallback,
 578         compat::not_null<StatePropagatorData*>     statePropagatorDataPtr,
 579         compat::not_null<EnergyElement*>           energyElementPtr,
 580         bool                                       hasReadEkinState)
 581 {
 582     auto forceElement = buildForces(
 583                 neighborSearchSignallerBuilder,
 584                 energySignallerBuilder,
 585                 statePropagatorDataPtr,
 586                 energyElementPtr);
 587
 588     // list of elements owned by the simulator composite object
 589     std::vector< std::unique_ptr<ISimulatorElement> >   elementsOwnershipList;
 590     // call list of the simulator composite object
 591     std::vector< compat::not_null<ISimulatorElement*> > elementCallList;
 592
 593     std::function<void()> needToCheckNumberOfBondedInteractions;
 594     if (inputrec->eI == eiMD)
 595     {
 596         auto computeGlobalsElement =
 597             std::make_unique< ComputeGlobalsElement<ComputeGlobalsAlgorithm::LeapFrog> >(
 598                     statePropagatorDataPtr, energyElementPtr, nstglobalcomm_, fplog, mdlog, cr,
 599                     inputrec, mdAtoms, nrnb, wcycle, fr,
 600                     &topologyHolder_->globalTopology(), constr, hasReadEkinState);
 601         topologyHolder_->registerClient(computeGlobalsElement.get());
 602         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
 603         trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElement.get()));
 604
 605         *checkBondedInteractionsCallback = computeGlobalsElement->getCheckNumberOfBondedInteractionsCallback();
 606
 607         auto propagator = std::make_unique< Propagator<IntegrationStep::LeapFrog> >(
 608                     inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
 609
 610         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
 611         addToCallList(statePropagatorDataPtr, elementCallList);  // we have a full microstate at time t here!
 612         addToCallListAndMove(std::move(propagator), elementCallList, elementsOwnershipList);
 613         if (constr)
 614         {
 615             auto constraintElement = std::make_unique< ConstraintsElement<ConstraintVariable::Positions> >(
 616                         constr, statePropagatorDataPtr, energyElementPtr, MASTER(cr),
 617                         fplog, inputrec, mdAtoms->mdatoms());
 618             auto constraintElementPtr = compat::make_not_null(constraintElement.get());
 619             energySignallerBuilder->registerSignallerClient(constraintElementPtr);
 620             trajectoryElementBuilder->registerSignallerClient(constraintElementPtr);
 621             loggingSignallerBuilder->registerSignallerClient(constraintElementPtr);
 622
 623             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 624         }
 625
 626         addToCallListAndMove(std::move(computeGlobalsElement), elementCallList, elementsOwnershipList);
 627         addToCallList(energyElementPtr, elementCallList);  // we have the energies at time t here!
 628     }
 629     else if (inputrec->eI == eiVV)
 630     {
 631         auto computeGlobalsElementAtFullTimeStep =
 632             std::make_unique< ComputeGlobalsElement<ComputeGlobalsAlgorithm::VelocityVerletAtFullTimeStep> >(
 633                     statePropagatorDataPtr, energyElementPtr, nstglobalcomm_, fplog, mdlog, cr,
 634                     inputrec, mdAtoms, nrnb, wcycle, fr,
 635                     &topologyHolder_->globalTopology(), constr, hasReadEkinState);
 636         topologyHolder_->registerClient(computeGlobalsElementAtFullTimeStep.get());
 637         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElementAtFullTimeStep.get()));
 638         trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElementAtFullTimeStep.get()));
 639
 640         auto computeGlobalsElementAfterCoordinateUpdate =
 641             std::make_unique<ComputeGlobalsElement <ComputeGlobalsAlgorithm::VelocityVerletAfterCoordinateUpdate> >(
 642                     statePropagatorDataPtr, energyElementPtr, nstglobalcomm_, fplog, mdlog, cr,
 643                     inputrec, mdAtoms, nrnb, wcycle, fr,
 644                     &topologyHolder_->globalTopology(), constr, hasReadEkinState);
 645         topologyHolder_->registerClient(computeGlobalsElementAfterCoordinateUpdate.get());
 646         energySignallerBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElementAfterCoordinateUpdate.get()));
 647         trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(computeGlobalsElementAfterCoordinateUpdate.get()));
 648
 649         *checkBondedInteractionsCallback = computeGlobalsElementAfterCoordinateUpdate->getCheckNumberOfBondedInteractionsCallback();
 650
 651         auto propagatorVelocities = std::make_unique< Propagator <IntegrationStep::VelocitiesOnly> >(
 652                     inputrec->delta_t * 0.5, statePropagatorDataPtr, mdAtoms, wcycle);
 653         auto propagatorVelocitiesAndPositions = std::make_unique< Propagator <IntegrationStep::VelocityVerletPositionsAndVelocities> >(
 654                     inputrec->delta_t, statePropagatorDataPtr, mdAtoms, wcycle);
 655
 656         addToCallListAndMove(std::move(forceElement), elementCallList, elementsOwnershipList);
 657         addToCallListAndMove(std::move(propagatorVelocities), elementCallList, elementsOwnershipList);
 658         if (constr)
 659         {
 660             auto constraintElement = std::make_unique< ConstraintsElement<ConstraintVariable::Velocities> >(
 661                         constr, statePropagatorDataPtr, energyElementPtr, MASTER(cr),
 662                         fplog, inputrec, mdAtoms->mdatoms());
 663             energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 664             trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 665             loggingSignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 666
 667             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 668         }
 669         addToCallListAndMove(std::move(computeGlobalsElementAtFullTimeStep), elementCallList, elementsOwnershipList);
 670         addToCallList(statePropagatorDataPtr, elementCallList);  // we have a full microstate at time t here!
 671         addToCallListAndMove(std::move(propagatorVelocitiesAndPositions), elementCallList, elementsOwnershipList);
 672         if (constr)
 673         {
 674             auto constraintElement = std::make_unique< ConstraintsElement<ConstraintVariable::Positions> >(
 675                         constr, statePropagatorDataPtr, energyElementPtr, MASTER(cr),
 676                         fplog, inputrec, mdAtoms->mdatoms());
 677             energySignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 678             trajectoryElementBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 679             loggingSignallerBuilder->registerSignallerClient(compat::make_not_null(constraintElement.get()));
 680
 681             addToCallListAndMove(std::move(constraintElement), elementCallList, elementsOwnershipList);
 682         }
 683         addToCallListAndMove(std::move(computeGlobalsElementAfterCoordinateUpdate), elementCallList, elementsOwnershipList);
 684         addToCallList(energyElementPtr, elementCallList);  // we have the energies at time t here!
 685     }
 686     else
 687     {
 688         gmx_fatal(FARGS, "Integrator not implemented for the modular simulator.");
 689     }
 690
 691     auto integrator = std::make_unique<CompositeSimulatorElement>(
 692                 std::move(elementCallList), std::move(elementsOwnershipList));
 693     // std::move *should* not be needed with c++-14, but clang-3.6 still requires it
 694     return std::move(integrator);
 695 }
 696
 697 bool ModularSimulator::isInputCompatible(
 698         bool                             exitOnFailure,
 699         const t_inputrec                *inputrec,
 700         bool                             doRerun,
 701         const gmx_vsite_t               *vsite,
 702         const gmx_multisim_t            *ms,
 703         const ReplicaExchangeParameters &replExParams,
 704         const t_fcdata                  *fcd,
 705         int                              nfile,
 706         const t_filenm                  *fnm,
 707         ObservablesHistory              *observablesHistory,
 708         const gmx_membed_t              *membed)
 709 {
 710     auto conditionalAssert =
 711         [exitOnFailure](bool condition, const char* message)
 712         {
 713             if (exitOnFailure)
 714             {
 715                 GMX_RELEASE_ASSERT(condition, message);
 716             }
 717             return condition;
 718         };
 719
 720     bool isInputCompatible = true;
 721
 722     // GMX_USE_MODULAR_SIMULATOR allows to use modular simulator also for non-standard uses,
 723     // such as the leap-frog integrator
 724     const auto modularSimulatorExplicitlyTurnedOn =
 725         (getenv("GMX_USE_MODULAR_SIMULATOR") != nullptr);
 726
 727     isInputCompatible = isInputCompatible && conditionalAssert(
 728                 inputrec->eI == eiMD || inputrec->eI == eiVV,
 729                 "Only integrators md and md-vv are supported by the modular simulator.");
 730     isInputCompatible = isInputCompatible && conditionalAssert(
 731                 inputrec->eI != eiMD || modularSimulatorExplicitlyTurnedOn,
 732                 "Set GMX_USE_MODULAR_SIMULATOR=ON to use the modular simulator with integrator md.");
 733     isInputCompatible = isInputCompatible && conditionalAssert(
 734                 !doRerun,
 735                 "Rerun is not supported by the modular simulator.");
 736     isInputCompatible = isInputCompatible && conditionalAssert(
 737                 inputrec->etc == etcNO,
 738                 "Temperature coupling is not supported by the modular simulator.");
 739     isInputCompatible = isInputCompatible && conditionalAssert(
 740                 inputrec->epc == epcNO,
 741                 "Pressure coupling is not supported by the modular simulator.");
 742     isInputCompatible = isInputCompatible && conditionalAssert(
 743                 !(inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec) || inputrecNvtTrotter(inputrec)),
 744                 "Legacy Trotter decomposition is not supported by the modular simulator.");
 745     isInputCompatible = isInputCompatible && conditionalAssert(
 746                 inputrec->efep == efepNO,
 747                 "Free energy calculation is not supported by the modular simulator.");
 748     isInputCompatible = isInputCompatible && conditionalAssert(
 749                 vsite == nullptr,
 750                 "Virtual sites are not supported by the modular simulator.");
 751     isInputCompatible = isInputCompatible && conditionalAssert(
 752                 !inputrec->bDoAwh,
 753                 "AWH is not supported by the modular simulator.");
 754     isInputCompatible = isInputCompatible && conditionalAssert(
 755                 ms == nullptr,
 756                 "Multi-sim are not supported by the modular simulator.");
 757     isInputCompatible = isInputCompatible && conditionalAssert(
 758                 replExParams.exchangeInterval == 0,
 759                 "Replica exchange is not supported by the modular simulator.");
 760     isInputCompatible = isInputCompatible && conditionalAssert(
 761                 fcd->disres.nsystems <= 1,
 762                 "Ensemble restraints are not supported by the modular simulator.");
 763     isInputCompatible = isInputCompatible && conditionalAssert(
 764                 !doSimulatedAnnealing(inputrec),
 765                 "Simulated annealing is not supported by the modular simulator.");
 766     isInputCompatible = isInputCompatible && conditionalAssert(
 767                 !inputrec->bSimTemp,
 768                 "Simulated tempering is not supported by the modular simulator.");
 769     isInputCompatible = isInputCompatible && conditionalAssert(
 770                 !inputrec->bExpanded,
 771                 "Expanded ensemble simulations are not supported by the modular simulator.");
 772     isInputCompatible = isInputCompatible && conditionalAssert(
 773                 !(opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr),
 774                 "Essential dynamics is not supported by the modular simulator.");
 775     isInputCompatible = isInputCompatible && conditionalAssert(
 776                 inputrec->eSwapCoords == eswapNO,
 777                 "Ion / water position swapping is not supported by the modular simulator.");
 778     isInputCompatible = isInputCompatible && conditionalAssert(
 779                 !inputrec->bIMD,
 780                 "Interactive MD is not supported by the modular simulator.");
 781     isInputCompatible = isInputCompatible && conditionalAssert(
 782                 membed == nullptr,
 783                 "Membrane embedding is not supported by the modular simulator.");
 784     // TODO: Change this to the boolean passed when we merge the user interface change for the GPU update.
 785     isInputCompatible = isInputCompatible && conditionalAssert(
 786                 getenv("GMX_UPDATE_CONSTRAIN_GPU") == nullptr,
 787                 "Integration on the GPU is not supported by the modular simulator.");
 788     // Modular simulator is centered around NS updates
 789     // TODO: think how to handle nstlist == 0
 790     isInputCompatible = isInputCompatible && conditionalAssert(
 791                 inputrec->nstlist != 0,
 792                 "Simulations without neighbor list update are not supported by the modular simulator.");
 793     isInputCompatible = isInputCompatible && conditionalAssert(
 794                 !GMX_FAHCORE,
 795                 "GMX_FAHCORE not supported by the modular simulator.");
 796
 797     return isInputCompatible;
 798 }
 799
 800 void ModularSimulator::checkInputForDisabledFunctionality()
 801 {
 802     isInputCompatible(
 803             true,
 804             inputrec, doRerun, vsite, ms, replExParams,
 805             fcd, nfile, fnm, observablesHistory, membed);
 806 }
 807
 808 SignallerCallbackPtr ModularSimulator::SignalHelper::registerLastStepCallback()
 809 {
 810     return std::make_unique<SignallerCallback>(
 811             [this](Step step, Time gmx_unused time){this->lastStep_ = step; });
 812 }
 813
 814 SignallerCallbackPtr ModularSimulator::SignalHelper::registerNSCallback()
 815 {
 816     return std::make_unique<SignallerCallback>(
 817             [this](Step step, Time gmx_unused time)
 818             {this->nextNSStep_ = step; });
 819 }
 820 }