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