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