src/gromacs/modularsimulator/simulatoralgorithm.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2020,2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief Defines the modular simulator algorithm
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "simulatoralgorithm.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/ewald/pme_pp.h"
  51 #include "gromacs/gmxlib/nrnb.h"
  52 #include "gromacs/listed_forces/listed_forces.h"
  53 #include "gromacs/mdlib/checkpointhandler.h"
  54 #include "gromacs/mdlib/constr.h"
  55 #include "gromacs/mdlib/energyoutput.h"
  56 #include "gromacs/mdlib/md_support.h"
  57 #include "gromacs/mdlib/mdatoms.h"
  58 #include "gromacs/mdlib/resethandler.h"
  59 #include "gromacs/mdlib/stat.h"
  60 #include "gromacs/mdrun/replicaexchange.h"
  61 #include "gromacs/mdrun/shellfc.h"
  62 #include "gromacs/mdrunutility/handlerestart.h"
  63 #include "gromacs/mdrunutility/printtime.h"
  64 #include "gromacs/mdtypes/commrec.h"
  65 #include "gromacs/mdtypes/fcdata.h"
  66 #include "gromacs/mdtypes/forcerec.h"
  67 #include "gromacs/mdtypes/inputrec.h"
  68 #include "gromacs/mdtypes/mdatom.h"
  69 #include "gromacs/mdtypes/mdrunoptions.h"
  70 #include "gromacs/mdtypes/observableshistory.h"
  71 #include "gromacs/nbnxm/nbnxm.h"
  72 #include "gromacs/timing/walltime_accounting.h"
  73 #include "gromacs/topology/topology.h"
  74 #include "gromacs/utility/cstringutil.h"
  75 #include "gromacs/utility/fatalerror.h"
  76
  77 #include "checkpointhelper.h"
  78 #include "domdechelper.h"
  79 #include "energydata.h"
  80 #include "freeenergyperturbationdata.h"
  81 #include "modularsimulator.h"
  82 #include "parrinellorahmanbarostat.h"
  83 #include "pmeloadbalancehelper.h"
  84 #include "propagator.h"
  85 #include "statepropagatordata.h"
  86 #include "velocityscalingtemperaturecoupling.h"
  87
  88 namespace gmx
  89 {
  90 ModularSimulatorAlgorithm::ModularSimulatorAlgorithm(std::string              topologyName,
  91                                                      FILE*                    fplog,
  92                                                      t_commrec*               cr,
  93                                                      const MDLogger&          mdlog,
  94                                                      const MdrunOptions&      mdrunOptions,
  95                                                      t_inputrec*              inputrec,
  96                                                      t_nrnb*                  nrnb,
  97                                                      gmx_wallcycle*           wcycle,
  98                                                      t_forcerec*              fr,
  99                                                      gmx_walltime_accounting* walltime_accounting) :
 100     taskIterator_(taskQueue_.end()),
 101     statePropagatorData_(nullptr),
 102     energyData_(nullptr),
 103     freeEnergyPerturbationData_(nullptr),
 104     step_(-1),
 105     runFinished_(false),
 106     topologyName_(std::move(topologyName)),
 107     fplog(fplog),
 108     cr(cr),
 109     mdlog(mdlog),
 110     mdrunOptions(mdrunOptions),
 111     inputrec(inputrec),
 112     nrnb(nrnb),
 113     wcycle(wcycle),
 114     fr(fr),
 115     walltime_accounting(walltime_accounting)
 116 {
 117     signalHelper_ = std::make_unique<SignalHelper>();
 118 }
 119
 120 void ModularSimulatorAlgorithm::setup()
 121 {
 122     simulatorSetup();
 123     for (auto& signaller : signallerList_)
 124     {
 125         signaller->setup();
 126     }
 127     if (domDecHelper_)
 128     {
 129         domDecHelper_->setup();
 130     }
 131
 132     for (auto& element : elementsOwnershipList_)
 133     {
 134         element->elementSetup();
 135     }
 136     statePropagatorData_->setup();
 137     if (pmeLoadBalanceHelper_)
 138     {
 139         // State must have been initialized so pmeLoadBalanceHelper_ gets a valid box
 140         pmeLoadBalanceHelper_->setup();
 141     }
 142 }
 143
 144 const SimulatorRunFunction* ModularSimulatorAlgorithm::getNextTask()
 145 {
 146     if (!taskQueue_.empty())
 147     {
 148         taskIterator_++;
 149     }
 150     if (taskIterator_ == taskQueue_.end())
 151     {
 152         if (runFinished_)
 153         {
 154             return nullptr;
 155         }
 156         updateTaskQueue();
 157         taskIterator_ = taskQueue_.begin();
 158     }
 159     return &*taskIterator_;
 160 }
 161
 162 void ModularSimulatorAlgorithm::updateTaskQueue()
 163 {
 164     // For now, we'll just clean the task queue and then re-populate
 165     // TODO: If tasks are periodic around updates of the task queue,
 166     //       we should reuse it instead
 167     taskQueue_.clear();
 168     populateTaskQueue();
 169 }
 170
 171 void ModularSimulatorAlgorithm::teardown()
 172 {
 173     for (auto& element : elementsOwnershipList_)
 174     {
 175         element->elementTeardown();
 176     }
 177     energyData_->teardown();
 178     if (pmeLoadBalanceHelper_)
 179     {
 180         pmeLoadBalanceHelper_->teardown();
 181     }
 182     simulatorTeardown();
 183 }
 184
 185 void ModularSimulatorAlgorithm::simulatorSetup()
 186 {
 187     if (!mdrunOptions.writeConfout)
 188     {
 189         // This is on by default, and the main known use case for
 190         // turning it off is for convenience in benchmarking, which is
 191         // something that should not show up in the general user
 192         // interface.
 193         GMX_LOG(mdlog.info)
 194                 .asParagraph()
 195                 .appendText(
 196                         "The -noconfout functionality is deprecated, and "
 197                         "may be removed in a future version.");
 198     }
 199
 200     if (MASTER(cr))
 201     {
 202         char        sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
 203         std::string timeString;
 204         fprintf(stderr, "starting mdrun '%s'\n", topologyName_.c_str());
 205         if (inputrec->nsteps >= 0)
 206         {
 207             timeString = formatString("%8.1f", static_cast<double>(inputrec->init_step + inputrec->nsteps)
 208                                                        * inputrec->delta_t);
 209         }
 210         else
 211         {
 212             timeString = "infinite";
 213         }
 214         if (inputrec->init_step > 0)
 215         {
 216             fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 217                     gmx_step_str(inputrec->init_step + inputrec->nsteps, sbuf), timeString.c_str(),
 218                     gmx_step_str(inputrec->init_step, sbuf2), inputrec->init_step * inputrec->delta_t);
 219         }
 220         else
 221         {
 222             fprintf(stderr, "%s steps, %s ps.\n", gmx_step_str(inputrec->nsteps, sbuf),
 223                     timeString.c_str());
 224         }
 225         fprintf(fplog, "\n");
 226     }
 227
 228     walltime_accounting_start_time(walltime_accounting);
 229     wallcycle_start(wcycle, ewcRUN);
 230     print_start(fplog, cr, walltime_accounting, "mdrun");
 231
 232     step_ = inputrec->init_step;
 233 }
 234
 235 void ModularSimulatorAlgorithm::simulatorTeardown()
 236 {
 237
 238     // Stop measuring walltime
 239     walltime_accounting_end_time(walltime_accounting);
 240
 241     if (!thisRankHasDuty(cr, DUTY_PME))
 242     {
 243         /* Tell the PME only node to finish */
 244         gmx_pme_send_finish(cr);
 245     }
 246
 247     walltime_accounting_set_nsteps_done(walltime_accounting, step_ - inputrec->init_step);
 248 }
 249
 250 void ModularSimulatorAlgorithm::preStep(Step step, Time gmx_unused time, bool isNeighborSearchingStep)
 251 {
 252     if (stopHandler_->stoppingAfterCurrentStep(isNeighborSearchingStep) && step != signalHelper_->lastStep_)
 253     {
 254         /*
 255          * Stop handler wants to stop after the current step, which was
 256          * not known when building the current task queue. This happens
 257          * e.g. when a stop is signalled by OS. We therefore want to purge
 258          * the task queue now, and re-schedule this step as last step.
 259          */
 260         // clear task queue
 261         taskQueue_.clear();
 262         // rewind step
 263         step_ = step;
 264         return;
 265     }
 266
 267     resetHandler_->setSignal(walltime_accounting);
 268     // This is a hack to avoid having to rewrite StopHandler to be a NeighborSearchSignaller
 269     // and accept the step as input. Eventually, we want to do that, but currently this would
 270     // require introducing NeighborSearchSignaller in the legacy do_md or a lot of code
 271     // duplication.
 272     stophandlerIsNSStep_    = isNeighborSearchingStep;
 273     stophandlerCurrentStep_ = step;
 274     stopHandler_->setSignal();
 275
 276     wallcycle_start(wcycle, ewcSTEP);
 277 }
 278
 279 void ModularSimulatorAlgorithm::postStep(Step step, Time gmx_unused time)
 280 {
 281     // Output stuff
 282     if (MASTER(cr))
 283     {
 284         if (do_per_step(step, inputrec->nstlog))
 285         {
 286             if (fflush(fplog) != 0)
 287             {
 288                 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
 289             }
 290         }
 291     }
 292     const bool do_verbose = mdrunOptions.verbose
 293                             && (step % mdrunOptions.verboseStepPrintInterval == 0
 294                                 || step == inputrec->init_step || step == signalHelper_->lastStep_);
 295     // Print the remaining wall clock time for the run
 296     if (MASTER(cr) && (do_verbose || gmx_got_usr_signal())
 297         && !(pmeLoadBalanceHelper_ && pmeLoadBalanceHelper_->pmePrinting()))
 298     {
 299         print_time(stderr, walltime_accounting, step, inputrec, cr);
 300     }
 301
 302     double cycles = wallcycle_stop(wcycle, ewcSTEP);
 303     if (DOMAINDECOMP(cr) && wcycle)
 304     {
 305         dd_cycles_add(cr->dd, static_cast<float>(cycles), ddCyclStep);
 306     }
 307
 308     resetHandler_->resetCounters(
 309             step, step - inputrec->init_step, mdlog, fplog, cr, fr->nbv.get(), nrnb, fr->pmedata,
 310             pmeLoadBalanceHelper_ ? pmeLoadBalanceHelper_->loadBalancingObject() : nullptr, wcycle,
 311             walltime_accounting);
 312 }
 313
 314 void ModularSimulatorAlgorithm::populateTaskQueue()
 315 {
 316     /*
 317      * The registerRunFunction emplaces functions to the task queue.
 318      * All elements are owned by the ModularSimulatorAlgorithm, as is the task queue.
 319      * Elements can hence register lambdas capturing their `this` pointers without expecting
 320      * life time issues, as the task queue and the elements are in the same scope.
 321      */
 322     auto registerRunFunction = [this](SimulatorRunFunction function) {
 323         taskQueue_.emplace_back(std::move(function));
 324     };
 325
 326     Time startTime = inputrec->init_t;
 327     Time timeStep  = inputrec->delta_t;
 328     Time time      = startTime + step_ * timeStep;
 329
 330     // Run an initial call to the signallers
 331     for (auto& signaller : signallerList_)
 332     {
 333         signaller->signal(step_, time);
 334     }
 335
 336     if (checkpointHelper_)
 337     {
 338         checkpointHelper_->run(step_, time);
 339     }
 340
 341     if (pmeLoadBalanceHelper_)
 342     {
 343         pmeLoadBalanceHelper_->run(step_, time);
 344     }
 345     if (domDecHelper_)
 346     {
 347         domDecHelper_->run(step_, time);
 348     }
 349
 350     do
 351     {
 352         // local variables for lambda capturing
 353         const int  step     = step_;
 354         const bool isNSStep = step == signalHelper_->nextNSStep_;
 355
 356         // register pre-step (task queue is local, so no problem with `this`)
 357         registerRunFunction([this, step, time, isNSStep]() { preStep(step, time, isNSStep); });
 358         // register elements for step
 359         for (auto& element : elementCallList_)
 360         {
 361             element->scheduleTask(step_, time, registerRunFunction);
 362         }
 363         // register post-step (task queue is local, so no problem with `this`)
 364         registerRunFunction([this, step, time]() { postStep(step, time); });
 365
 366         // prepare next step
 367         step_++;
 368         time = startTime + step_ * timeStep;
 369         for (auto& signaller : signallerList_)
 370         {
 371             signaller->signal(step_, time);
 372         }
 373     } while (step_ != signalHelper_->nextNSStep_ && step_ <= signalHelper_->lastStep_);
 374
 375     runFinished_ = (step_ > signalHelper_->lastStep_);
 376
 377     if (runFinished_)
 378     {
 379         // task queue is local, so no problem with `this`
 380         registerRunFunction([this]() { teardown(); });
 381     }
 382 }
 383
 384 ModularSimulatorAlgorithmBuilder::ModularSimulatorAlgorithmBuilder(
 385         compat::not_null<LegacySimulatorData*>    legacySimulatorData,
 386         std::unique_ptr<ReadCheckpointDataHolder> checkpointDataHolder) :
 387     legacySimulatorData_(legacySimulatorData),
 388     signals_(std::make_unique<SimulationSignals>()),
 389     elementAdditionHelper_(this),
 390     globalCommunicationHelper_(computeGlobalCommunicationPeriod(legacySimulatorData->mdlog,
 391                                                                 legacySimulatorData->inputrec,
 392                                                                 legacySimulatorData->cr),
 393                                signals_.get()),
 394     checkpointHelperBuilder_(std::move(checkpointDataHolder),
 395                              legacySimulatorData->startingBehavior,
 396                              legacySimulatorData->cr)
 397 {
 398     if (legacySimulatorData->inputrec->efep != efepNO)
 399     {
 400         freeEnergyPerturbationData_ = std::make_unique<FreeEnergyPerturbationData>(
 401                 legacySimulatorData->fplog, legacySimulatorData->inputrec, legacySimulatorData->mdAtoms);
 402     }
 403
 404     statePropagatorData_ = std::make_unique<StatePropagatorData>(
 405             legacySimulatorData->top_global->natoms, legacySimulatorData->fplog, legacySimulatorData->cr,
 406             legacySimulatorData->state_global, legacySimulatorData->fr->nbv->useGpu(),
 407             legacySimulatorData->fr->bMolPBC, legacySimulatorData->mdrunOptions.writeConfout,
 408             opt2fn("-c", legacySimulatorData->nfile, legacySimulatorData->fnm), legacySimulatorData->inputrec,
 409             legacySimulatorData->mdAtoms->mdatoms(), legacySimulatorData->top_global);
 410
 411     // Multi sim is turned off
 412     const bool simulationsShareState = false;
 413
 414     energyData_ = std::make_unique<EnergyData>(
 415             statePropagatorData_.get(), freeEnergyPerturbationData_.get(), legacySimulatorData->top_global,
 416             legacySimulatorData->inputrec, legacySimulatorData->mdAtoms, legacySimulatorData->enerd,
 417             legacySimulatorData->ekind, legacySimulatorData->constr, legacySimulatorData->fplog,
 418             legacySimulatorData->fr->fcdata.get(), legacySimulatorData->mdModulesNotifier,
 419             MASTER(legacySimulatorData->cr), legacySimulatorData->observablesHistory,
 420             legacySimulatorData->startingBehavior, simulationsShareState);
 421 }
 422
 423 ModularSimulatorAlgorithm ModularSimulatorAlgorithmBuilder::build()
 424 {
 425     if (algorithmHasBeenBuilt_)
 426     {
 427         throw SimulationAlgorithmSetupError(
 428                 "Tried to build ModularSimulationAlgorithm more than once.");
 429     }
 430     algorithmHasBeenBuilt_ = true;
 431
 432     // Connect propagators with thermostat / barostat
 433     for (const auto& thermostatRegistration : thermostatRegistrationFunctions_)
 434     {
 435         for (const auto& connection : propagatorThermostatConnections_)
 436         {
 437             thermostatRegistration(connection);
 438         }
 439     }
 440     for (const auto& barostatRegistration : barostatRegistrationFunctions_)
 441     {
 442         for (const auto& connection : propagatorBarostatConnections_)
 443         {
 444             barostatRegistration(connection);
 445         }
 446     }
 447
 448     ModularSimulatorAlgorithm algorithm(
 449             *(legacySimulatorData_->top_global->name), legacySimulatorData_->fplog,
 450             legacySimulatorData_->cr, legacySimulatorData_->mdlog, legacySimulatorData_->mdrunOptions,
 451             legacySimulatorData_->inputrec, legacySimulatorData_->nrnb, legacySimulatorData_->wcycle,
 452             legacySimulatorData_->fr, legacySimulatorData_->walltime_accounting);
 453     registerWithInfrastructureAndSignallers(algorithm.signalHelper_.get());
 454     algorithm.statePropagatorData_        = std::move(statePropagatorData_);
 455     algorithm.energyData_                 = std::move(energyData_);
 456     algorithm.freeEnergyPerturbationData_ = std::move(freeEnergyPerturbationData_);
 457     algorithm.signals_                    = std::move(signals_);
 458
 459     // Multi sim is turned off
 460     const bool simulationsShareState = false;
 461
 462     // Build stop handler
 463     algorithm.stopHandler_ = legacySimulatorData_->stopHandlerBuilder->getStopHandlerMD(
 464             compat::not_null<SimulationSignal*>(
 465                     &(*globalCommunicationHelper_.simulationSignals())[eglsSTOPCOND]),
 466             simulationsShareState, MASTER(legacySimulatorData_->cr),
 467             legacySimulatorData_->inputrec->nstlist, legacySimulatorData_->mdrunOptions.reproducible,
 468             globalCommunicationHelper_.nstglobalcomm(), legacySimulatorData_->mdrunOptions.maximumHoursToRun,
 469             legacySimulatorData_->inputrec->nstlist == 0, legacySimulatorData_->fplog,
 470             algorithm.stophandlerCurrentStep_, algorithm.stophandlerIsNSStep_,
 471             legacySimulatorData_->walltime_accounting);
 472
 473     // Build reset handler
 474     const bool simulationsShareResetCounters = false;
 475     algorithm.resetHandler_                  = std::make_unique<ResetHandler>(
 476             compat::make_not_null<SimulationSignal*>(
 477                     &(*globalCommunicationHelper_.simulationSignals())[eglsRESETCOUNTERS]),
 478             simulationsShareResetCounters, legacySimulatorData_->inputrec->nsteps,
 479             MASTER(legacySimulatorData_->cr), legacySimulatorData_->mdrunOptions.timingOptions.resetHalfway,
 480             legacySimulatorData_->mdrunOptions.maximumHoursToRun, legacySimulatorData_->mdlog,
 481             legacySimulatorData_->wcycle, legacySimulatorData_->walltime_accounting);
 482
 483     // Build topology holder
 484     algorithm.topologyHolder_ = topologyHolderBuilder_.build(
 485             *legacySimulatorData_->top_global, legacySimulatorData_->cr,
 486             legacySimulatorData_->inputrec, legacySimulatorData_->fr, legacySimulatorData_->mdAtoms,
 487             legacySimulatorData_->constr, legacySimulatorData_->vsite);
 488
 489     // Build PME load balance helper
 490     if (PmeLoadBalanceHelper::doPmeLoadBalancing(legacySimulatorData_->mdrunOptions,
 491                                                  legacySimulatorData_->inputrec,
 492                                                  legacySimulatorData_->fr))
 493     {
 494         algorithm.pmeLoadBalanceHelper_ = std::make_unique<PmeLoadBalanceHelper>(
 495                 legacySimulatorData_->mdrunOptions.verbose, algorithm.statePropagatorData_.get(),
 496                 legacySimulatorData_->fplog, legacySimulatorData_->cr, legacySimulatorData_->mdlog,
 497                 legacySimulatorData_->inputrec, legacySimulatorData_->wcycle, legacySimulatorData_->fr);
 498         registerWithInfrastructureAndSignallers(algorithm.pmeLoadBalanceHelper_.get());
 499     }
 500     // Build domdec helper
 501     if (DOMAINDECOMP(legacySimulatorData_->cr))
 502     {
 503         algorithm.domDecHelper_ = std::make_unique<DomDecHelper>(
 504                 legacySimulatorData_->mdrunOptions.verbose,
 505                 legacySimulatorData_->mdrunOptions.verboseStepPrintInterval,
 506                 algorithm.statePropagatorData_.get(), algorithm.freeEnergyPerturbationData_.get(),
 507                 algorithm.topologyHolder_.get(),
 508                 globalCommunicationHelper_.moveCheckBondedInteractionsCallback(),
 509                 globalCommunicationHelper_.nstglobalcomm(), legacySimulatorData_->fplog,
 510                 legacySimulatorData_->cr, legacySimulatorData_->mdlog, legacySimulatorData_->constr,
 511                 legacySimulatorData_->inputrec, legacySimulatorData_->mdAtoms, legacySimulatorData_->nrnb,
 512                 legacySimulatorData_->wcycle, legacySimulatorData_->fr, legacySimulatorData_->vsite,
 513                 legacySimulatorData_->imdSession, legacySimulatorData_->pull_work);
 514         registerWithInfrastructureAndSignallers(algorithm.domDecHelper_.get());
 515     }
 516
 517     // Build trajectory element
 518     auto trajectoryElement = trajectoryElementBuilder_.build(
 519             legacySimulatorData_->fplog, legacySimulatorData_->nfile, legacySimulatorData_->fnm,
 520             legacySimulatorData_->mdrunOptions, legacySimulatorData_->cr,
 521             legacySimulatorData_->outputProvider, legacySimulatorData_->mdModulesNotifier,
 522             legacySimulatorData_->inputrec, legacySimulatorData_->top_global,
 523             legacySimulatorData_->oenv, legacySimulatorData_->wcycle,
 524             legacySimulatorData_->startingBehavior, simulationsShareState);
 525     registerWithInfrastructureAndSignallers(trajectoryElement.get());
 526
 527     // Build free energy element
 528     std::unique_ptr<FreeEnergyPerturbationData::Element> freeEnergyPerturbationElement = nullptr;
 529     if (algorithm.freeEnergyPerturbationData_)
 530     {
 531         freeEnergyPerturbationElement = std::make_unique<FreeEnergyPerturbationData::Element>(
 532                 algorithm.freeEnergyPerturbationData_.get(),
 533                 legacySimulatorData_->inputrec->fepvals->delta_lambda);
 534         registerWithInfrastructureAndSignallers(freeEnergyPerturbationElement.get());
 535     }
 536
 537     // Build checkpoint helper (do this last so everyone else can be a checkpoint client!)
 538     {
 539         checkpointHelperBuilder_.setCheckpointHandler(std::make_unique<CheckpointHandler>(
 540                 compat::make_not_null<SimulationSignal*>(&(*algorithm.signals_)[eglsCHKPT]),
 541                 simulationsShareState, legacySimulatorData_->inputrec->nstlist == 0,
 542                 MASTER(legacySimulatorData_->cr), legacySimulatorData_->mdrunOptions.writeConfout,
 543                 legacySimulatorData_->mdrunOptions.checkpointOptions.period));
 544         algorithm.checkpointHelper_ = checkpointHelperBuilder_.build(
 545                 legacySimulatorData_->inputrec->init_step, trajectoryElement.get(),
 546                 legacySimulatorData_->fplog, legacySimulatorData_->cr,
 547                 legacySimulatorData_->observablesHistory, legacySimulatorData_->walltime_accounting,
 548                 legacySimulatorData_->state_global, legacySimulatorData_->mdrunOptions.writeConfout);
 549         registerWithInfrastructureAndSignallers(algorithm.checkpointHelper_.get());
 550     }
 551
 552     // Build signallers
 553     {
 554         /* Signallers need to be called in an exact order. Some signallers are clients
 555          * of other signallers, which requires the clients signallers to be called
 556          * _after_ any signaller they are registered to - otherwise, they couldn't
 557          * adapt their behavior to the information they got signalled.
 558          *
 559          * Signallers being clients of other signallers require registration.
 560          * That registration happens during construction, which in turn means that
 561          * we want to construct the signallers in the reverse order of their later
 562          * call order.
 563          *
 564          * For the above reasons, the `addSignaller` lambda defined below emplaces
 565          * added signallers at the beginning of the signaller list, which will yield
 566          * a signaller list which is inverse to the build order (and hence equal to
 567          * the intended call order).
 568          */
 569         auto addSignaller = [this, &algorithm](auto signaller) {
 570             registerWithInfrastructureAndSignallers(signaller.get());
 571             algorithm.signallerList_.emplace(algorithm.signallerList_.begin(), std::move(signaller));
 572         };
 573         const auto* inputrec = legacySimulatorData_->inputrec;
 574         addSignaller(energySignallerBuilder_.build(
 575                 inputrec->nstcalcenergy, inputrec->fepvals->nstdhdl, inputrec->nstpcouple));
 576         addSignaller(trajectorySignallerBuilder_.build(
 577                 inputrec->nstxout, inputrec->nstvout, inputrec->nstfout,
 578                 inputrec->nstxout_compressed, trajectoryElement->tngBoxOut(),
 579                 trajectoryElement->tngLambdaOut(), trajectoryElement->tngBoxOutCompressed(),
 580                 trajectoryElement->tngLambdaOutCompressed(), inputrec->nstenergy));
 581         addSignaller(loggingSignallerBuilder_.build(inputrec->nstlog, inputrec->init_step,
 582                                                     legacySimulatorData_->startingBehavior));
 583         addSignaller(lastStepSignallerBuilder_.build(inputrec->nsteps, inputrec->init_step,
 584                                                      algorithm.stopHandler_.get()));
 585         addSignaller(neighborSearchSignallerBuilder_.build(inputrec->nstlist, inputrec->init_step,
 586                                                            inputrec->init_t));
 587     }
 588
 589     // Create element list
 590     // Checkpoint helper needs to be in the call list (as first element!) to react to last step
 591     algorithm.elementCallList_.emplace_back(algorithm.checkpointHelper_.get());
 592     // Next, update the free energy lambda vector if needed
 593     if (freeEnergyPerturbationElement)
 594     {
 595         algorithm.elementsOwnershipList_.emplace_back(std::move(freeEnergyPerturbationElement));
 596         algorithm.elementCallList_.emplace_back(algorithm.elementsOwnershipList_.back().get());
 597     }
 598     // Then, move the built algorithm
 599     algorithm.elementsOwnershipList_.insert(algorithm.elementsOwnershipList_.end(),
 600                                             std::make_move_iterator(elements_.begin()),
 601                                             std::make_move_iterator(elements_.end()));
 602     algorithm.elementCallList_.insert(algorithm.elementCallList_.end(),
 603                                       std::make_move_iterator(callList_.begin()),
 604                                       std::make_move_iterator(callList_.end()));
 605     // Finally, all trajectory writing is happening after the step
 606     // (relevant data was stored by elements through energy signaller)
 607     algorithm.elementsOwnershipList_.emplace_back(std::move(trajectoryElement));
 608     algorithm.elementCallList_.emplace_back(algorithm.elementsOwnershipList_.back().get());
 609
 610     algorithm.setup();
 611     return algorithm;
 612 }
 613
 614 ISimulatorElement* ModularSimulatorAlgorithmBuilder::addElementToSimulatorAlgorithm(
 615         std::unique_ptr<ISimulatorElement> element)
 616 {
 617     elements_.emplace_back(std::move(element));
 618     return elements_.back().get();
 619 }
 620
 621 bool ModularSimulatorAlgorithmBuilder::elementExists(const ISimulatorElement* element) const
 622 {
 623     // Check whether element exists in element list
 624     if (std::any_of(elements_.begin(), elements_.end(),
 625                     [element](auto& existingElement) { return element == existingElement.get(); }))
 626     {
 627         return true;
 628     }
 629     // Check whether element exists in other places controlled by *this
 630     return (statePropagatorData_->element() == element || energyData_->element() == element
 631             || (freeEnergyPerturbationData_ && freeEnergyPerturbationData_->element() == element));
 632 }
 633
 634 void ModularSimulatorAlgorithmBuilder::addElementToSetupTeardownList(ISimulatorElement* element)
 635 {
 636     // Add element if it's not already in the list
 637     if (std::find(setupAndTeardownList_.begin(), setupAndTeardownList_.end(), element)
 638         == setupAndTeardownList_.end())
 639     {
 640         setupAndTeardownList_.emplace_back(element);
 641     }
 642 }
 643
 644 std::optional<SignallerCallback> ModularSimulatorAlgorithm::SignalHelper::registerLastStepCallback()
 645 {
 646     return [this](Step step, Time gmx_unused time) { this->lastStep_ = step; };
 647 }
 648
 649 std::optional<SignallerCallback> ModularSimulatorAlgorithm::SignalHelper::registerNSCallback()
 650 {
 651     return [this](Step step, Time gmx_unused time) { this->nextNSStep_ = step; };
 652 }
 653
 654 GlobalCommunicationHelper::GlobalCommunicationHelper(int nstglobalcomm, SimulationSignals* simulationSignals) :
 655     nstglobalcomm_(nstglobalcomm),
 656     simulationSignals_(simulationSignals)
 657 {
 658 }
 659
 660 int GlobalCommunicationHelper::nstglobalcomm() const
 661 {
 662     return nstglobalcomm_;
 663 }
 664
 665 SimulationSignals* GlobalCommunicationHelper::simulationSignals()
 666 {
 667     return simulationSignals_;
 668 }
 669
 670 void GlobalCommunicationHelper::setCheckBondedInteractionsCallback(CheckBondedInteractionsCallback callback)
 671 {
 672     checkBondedInteractionsCallback_ = std::move(callback);
 673 }
 674
 675 CheckBondedInteractionsCallback GlobalCommunicationHelper::moveCheckBondedInteractionsCallback()
 676 {
 677     if (!checkBondedInteractionsCallback_)
 678     {
 679         throw SimulationAlgorithmSetupError(
 680                 "Requested CheckBondedInteractionsCallback before it was set.");
 681     }
 682     return *checkBondedInteractionsCallback_;
 683 }
 684
 685 ModularSimulatorAlgorithmBuilderHelper::ModularSimulatorAlgorithmBuilderHelper(
 686         ModularSimulatorAlgorithmBuilder* builder) :
 687     builder_(builder)
 688 {
 689 }
 690
 691 ISimulatorElement* ModularSimulatorAlgorithmBuilderHelper::storeElement(std::unique_ptr<ISimulatorElement> element)
 692 {
 693     return builder_->addElementToSimulatorAlgorithm(std::move(element));
 694 }
 695
 696 bool ModularSimulatorAlgorithmBuilderHelper::elementIsStored(const ISimulatorElement* element) const
 697 {
 698     return builder_->elementExists(element);
 699 }
 700
 701 std::optional<std::any> ModularSimulatorAlgorithmBuilderHelper::getStoredValue(const std::string& key) const
 702 {
 703     const auto iter = values_.find(key);
 704     if (iter == values_.end())
 705     {
 706         return std::nullopt;
 707     }
 708     else
 709     {
 710         return iter->second;
 711     }
 712 }
 713
 714 void ModularSimulatorAlgorithmBuilderHelper::registerThermostat(
 715         std::function<void(const PropagatorThermostatConnection&)> registrationFunction)
 716 {
 717     builder_->thermostatRegistrationFunctions_.emplace_back(std::move(registrationFunction));
 718 }
 719
 720 void ModularSimulatorAlgorithmBuilderHelper::registerBarostat(
 721         std::function<void(const PropagatorBarostatConnection&)> registrationFunction)
 722 {
 723     builder_->barostatRegistrationFunctions_.emplace_back(std::move(registrationFunction));
 724 }
 725
 726 void ModularSimulatorAlgorithmBuilderHelper::registerWithThermostat(PropagatorThermostatConnection connectionData)
 727 {
 728     builder_->propagatorThermostatConnections_.emplace_back(std::move(connectionData));
 729 }
 730
 731 void ModularSimulatorAlgorithmBuilderHelper::registerWithBarostat(PropagatorBarostatConnection connectionData)
 732 {
 733     builder_->propagatorBarostatConnections_.emplace_back(std::move(connectionData));
 734 }
 735
 736
 737 } // namespace gmx