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