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/freeenergy.h"
  63 #include "gromacs/mdrunutility/handlerestart.h"
  64 #include "gromacs/mdrunutility/printtime.h"
  65 #include "gromacs/mdtypes/commrec.h"
  66 #include "gromacs/mdtypes/fcdata.h"
  67 #include "gromacs/mdtypes/forcerec.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/mdatom.h"
  70 #include "gromacs/mdtypes/mdrunoptions.h"
  71 #include "gromacs/mdtypes/observableshistory.h"
  72 #include "gromacs/nbnxm/nbnxm.h"
  73 #include "gromacs/timing/walltime_accounting.h"
  74 #include "gromacs/topology/topology.h"
  75 #include "gromacs/utility/cstringutil.h"
  76 #include "gromacs/utility/fatalerror.h"
  77
  78 #include "checkpointhelper.h"
  79 #include "domdechelper.h"
  80 #include "energydata.h"
  81 #include "firstorderpressurecoupling.h"
  82 #include "freeenergyperturbationdata.h"
  83 #include "modularsimulator.h"
  84 #include "pmeloadbalancehelper.h"
  85 #include "propagator.h"
  86 #include "referencetemperaturemanager.h"
  87 #include "statepropagatordata.h"
  88
  89 namespace gmx
  90 {
  91 ModularSimulatorAlgorithm::ModularSimulatorAlgorithm(std::string              topologyName,
  92                                                      FILE*                    fplog,
  93                                                      t_commrec*               cr,
  94                                                      const MDLogger&          mdlog,
  95                                                      const MdrunOptions&      mdrunOptions,
  96                                                      const t_inputrec*        inputrec,
  97                                                      t_nrnb*                  nrnb,
  98                                                      gmx_wallcycle*           wcycle,
  99                                                      t_forcerec*              fr,
 100                                                      gmx_walltime_accounting* walltime_accounting) :
 101     taskIterator_(taskQueue_.end()),
 102     statePropagatorData_(nullptr),
 103     energyData_(nullptr),
 104     freeEnergyPerturbationData_(nullptr),
 105     step_(-1),
 106     runFinished_(false),
 107     topologyName_(std::move(topologyName)),
 108     fplog(fplog),
 109     cr(cr),
 110     mdlog(mdlog),
 111     mdrunOptions(mdrunOptions),
 112     inputrec(inputrec),
 113     nrnb(nrnb),
 114     wcycle(wcycle),
 115     fr(fr),
 116     walltime_accounting(walltime_accounting)
 117 {
 118     signalHelper_ = std::make_unique<SignalHelper>();
 119 }
 120
 121 void ModularSimulatorAlgorithm::setup()
 122 {
 123     simulatorSetup();
 124     for (auto& signaller : signallerList_)
 125     {
 126         signaller->setup();
 127     }
 128     if (domDecHelper_)
 129     {
 130         domDecHelper_->setup();
 131     }
 132
 133     for (auto& element : elementSetupTeardownList_)
 134     {
 135         element->elementSetup();
 136     }
 137     statePropagatorData_->setup();
 138     if (pmeLoadBalanceHelper_)
 139     {
 140         // State must have been initialized so pmeLoadBalanceHelper_ gets a valid box
 141         pmeLoadBalanceHelper_->setup();
 142     }
 143 }
 144
 145 const SimulatorRunFunction* ModularSimulatorAlgorithm::getNextTask()
 146 {
 147     if (!taskQueue_.empty())
 148     {
 149         taskIterator_++;
 150     }
 151     if (taskIterator_ == taskQueue_.end())
 152     {
 153         if (runFinished_)
 154         {
 155             return nullptr;
 156         }
 157         updateTaskQueue();
 158         taskIterator_ = taskQueue_.begin();
 159     }
 160     return &*taskIterator_;
 161 }
 162
 163 void ModularSimulatorAlgorithm::updateTaskQueue()
 164 {
 165     // For now, we'll just clean the task queue and then re-populate
 166     // TODO: If tasks are periodic around updates of the task queue,
 167     //       we should reuse it instead
 168     taskQueue_.clear();
 169     populateTaskQueue();
 170 }
 171
 172 void ModularSimulatorAlgorithm::teardown()
 173 {
 174     for (auto& element : elementSetupTeardownList_)
 175     {
 176         element->elementTeardown();
 177     }
 178     energyData_->teardown();
 179     if (pmeLoadBalanceHelper_)
 180     {
 181         pmeLoadBalanceHelper_->teardown();
 182     }
 183     simulatorTeardown();
 184 }
 185
 186 void ModularSimulatorAlgorithm::simulatorSetup()
 187 {
 188     if (!mdrunOptions.writeConfout)
 189     {
 190         // This is on by default, and the main known use case for
 191         // turning it off is for convenience in benchmarking, which is
 192         // something that should not show up in the general user
 193         // interface.
 194         GMX_LOG(mdlog.info)
 195                 .asParagraph()
 196                 .appendText(
 197                         "The -noconfout functionality is deprecated, and "
 198                         "may be removed in a future version.");
 199     }
 200
 201     if (MASTER(cr))
 202     {
 203         char        sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
 204         std::string timeString;
 205         fprintf(stderr, "starting mdrun '%s'\n", topologyName_.c_str());
 206         if (inputrec->nsteps >= 0)
 207         {
 208             timeString = formatString(
 209                     "%8.1f", static_cast<double>(inputrec->init_step + inputrec->nsteps) * inputrec->delta_t);
 210         }
 211         else
 212         {
 213             timeString = "infinite";
 214         }
 215         if (inputrec->init_step > 0)
 216         {
 217             fprintf(stderr,
 218                     "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 219                     gmx_step_str(inputrec->init_step + inputrec->nsteps, sbuf),
 220                     timeString.c_str(),
 221                     gmx_step_str(inputrec->init_step, sbuf2),
 222                     inputrec->init_step * inputrec->delta_t);
 223         }
 224         else
 225         {
 226             fprintf(stderr, "%s steps, %s ps.\n", gmx_step_str(inputrec->nsteps, sbuf), timeString.c_str());
 227         }
 228         fprintf(fplog, "\n");
 229     }
 230
 231     walltime_accounting_start_time(walltime_accounting);
 232     wallcycle_start(wcycle, WallCycleCounter::Run);
 233     print_start(fplog, cr, walltime_accounting, "mdrun");
 234
 235     step_ = inputrec->init_step;
 236 }
 237
 238 void ModularSimulatorAlgorithm::simulatorTeardown()
 239 {
 240
 241     // Stop measuring walltime
 242     walltime_accounting_end_time(walltime_accounting);
 243
 244     if (!thisRankHasDuty(cr, DUTY_PME))
 245     {
 246         /* Tell the PME only node to finish */
 247         gmx_pme_send_finish(cr);
 248     }
 249
 250     walltime_accounting_set_nsteps_done(walltime_accounting, step_ - inputrec->init_step);
 251 }
 252
 253 void ModularSimulatorAlgorithm::preStep(Step step, Time gmx_unused time, bool isNeighborSearchingStep)
 254 {
 255     if (stopHandler_->stoppingAfterCurrentStep(isNeighborSearchingStep) && step != signalHelper_->lastStep_)
 256     {
 257         /*
 258          * Stop handler wants to stop after the current step, which was
 259          * not known when building the current task queue. This happens
 260          * e.g. when a stop is signalled by OS. We therefore want to purge
 261          * the task queue now, and re-schedule this step as last step.
 262          */
 263         // clear task queue
 264         taskQueue_.clear();
 265         // rewind step
 266         step_ = step;
 267         return;
 268     }
 269
 270     resetHandler_->setSignal(walltime_accounting);
 271     // This is a hack to avoid having to rewrite StopHandler to be a NeighborSearchSignaller
 272     // and accept the step as input. Eventually, we want to do that, but currently this would
 273     // require introducing NeighborSearchSignaller in the legacy do_md or a lot of code
 274     // duplication.
 275     stophandlerIsNSStep_    = isNeighborSearchingStep;
 276     stophandlerCurrentStep_ = step;
 277     stopHandler_->setSignal();
 278
 279     wallcycle_start(wcycle, WallCycleCounter::Step);
 280 }
 281
 282 void ModularSimulatorAlgorithm::postStep(Step step, Time gmx_unused time)
 283 {
 284     // Output stuff
 285     if (MASTER(cr))
 286     {
 287         if (do_per_step(step, inputrec->nstlog))
 288         {
 289             if (fflush(fplog) != 0)
 290             {
 291                 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
 292             }
 293         }
 294     }
 295     const bool do_verbose = mdrunOptions.verbose
 296                             && (step % mdrunOptions.verboseStepPrintInterval == 0
 297                                 || step == inputrec->init_step || step == signalHelper_->lastStep_);
 298     // Print the remaining wall clock time for the run
 299     if (MASTER(cr) && (do_verbose || gmx_got_usr_signal())
 300         && !(pmeLoadBalanceHelper_ && pmeLoadBalanceHelper_->pmePrinting()))
 301     {
 302         print_time(stderr, walltime_accounting, step, inputrec, cr);
 303     }
 304
 305     double cycles = wallcycle_stop(wcycle, WallCycleCounter::Step);
 306     if (haveDDAtomOrdering(*cr) && wcycle)
 307     {
 308         dd_cycles_add(cr->dd, static_cast<float>(cycles), ddCyclStep);
 309     }
 310
 311     resetHandler_->resetCounters(step,
 312                                  step - inputrec->init_step,
 313                                  mdlog,
 314                                  fplog,
 315                                  cr,
 316                                  fr->nbv.get(),
 317                                  nrnb,
 318                                  fr->pmedata,
 319                                  pmeLoadBalanceHelper_ ? pmeLoadBalanceHelper_->loadBalancingObject() : nullptr,
 320                                  wcycle,
 321                                  walltime_accounting);
 322 }
 323
 324 void ModularSimulatorAlgorithm::populateTaskQueue()
 325 {
 326     /*
 327      * The registerRunFunction emplaces functions to the task queue.
 328      * All elements are owned by the ModularSimulatorAlgorithm, as is the task queue.
 329      * Elements can hence register lambdas capturing their `this` pointers without expecting
 330      * life time issues, as the task queue and the elements are in the same scope.
 331      */
 332     auto registerRunFunction = [this](SimulatorRunFunction function) {
 333         taskQueue_.emplace_back(std::move(function));
 334     };
 335
 336     Time startTime = inputrec->init_t;
 337     Time timeStep  = inputrec->delta_t;
 338     Time time      = startTime + step_ * timeStep;
 339
 340     // Run an initial call to the signallers
 341     for (auto& signaller : signallerList_)
 342     {
 343         signaller->signal(step_, time);
 344     }
 345
 346     if (checkpointHelper_)
 347     {
 348         checkpointHelper_->run(step_, time);
 349     }
 350
 351     if (pmeLoadBalanceHelper_)
 352     {
 353         pmeLoadBalanceHelper_->run(step_, time);
 354     }
 355     if (domDecHelper_)
 356     {
 357         domDecHelper_->run(step_, time);
 358     }
 359
 360     do
 361     {
 362         // local variables for lambda capturing
 363         const int  step     = step_;
 364         const bool isNSStep = step == signalHelper_->nextNSStep_;
 365
 366         // register pre-step (task queue is local, so no problem with `this`)
 367         registerRunFunction([this, step, time, isNSStep]() { preStep(step, time, isNSStep); });
 368         // register pre step functions
 369         for (const auto& schedulingFunction : preStepScheduling_)
 370         {
 371             schedulingFunction(step_, time, registerRunFunction);
 372         }
 373         // register elements for step
 374         for (auto& element : elementCallList_)
 375         {
 376             element->scheduleTask(step_, time, registerRunFunction);
 377         }
 378         // register post step functions
 379         for (const auto& schedulingFunction : postStepScheduling_)
 380         {
 381             schedulingFunction(step_, time, registerRunFunction);
 382         }
 383         // register post-step (task queue is local, so no problem with `this`)
 384         registerRunFunction([this, step, time]() { postStep(step, time); });
 385
 386         // prepare next step
 387         step_++;
 388         time = startTime + step_ * timeStep;
 389         for (auto& signaller : signallerList_)
 390         {
 391             signaller->signal(step_, time);
 392         }
 393     } while (step_ != signalHelper_->nextNSStep_ && step_ <= signalHelper_->lastStep_);
 394
 395     runFinished_ = (step_ > signalHelper_->lastStep_);
 396
 397     if (runFinished_)
 398     {
 399         // task queue is local, so no problem with `this`
 400         registerRunFunction([this]() { teardown(); });
 401     }
 402 }
 403
 404 ModularSimulatorAlgorithmBuilder::ModularSimulatorAlgorithmBuilder(
 405         compat::not_null<LegacySimulatorData*>    legacySimulatorData,
 406         std::unique_ptr<ReadCheckpointDataHolder> checkpointDataHolder) :
 407     legacySimulatorData_(legacySimulatorData),
 408     signals_(std::make_unique<SimulationSignals>()),
 409     elementAdditionHelper_(this),
 410     globalCommunicationHelper_(computeGlobalCommunicationPeriod(legacySimulatorData->mdlog,
 411                                                                 legacySimulatorData->inputrec,
 412                                                                 legacySimulatorData->cr),
 413                                signals_.get()),
 414     observablesReducer_(legacySimulatorData->observablesReducerBuilder->build()),
 415     checkpointHelperBuilder_(std::move(checkpointDataHolder),
 416                              legacySimulatorData->startingBehavior,
 417                              legacySimulatorData->cr)
 418 {
 419     if (legacySimulatorData->inputrec->efep != FreeEnergyPerturbationType::No)
 420     {
 421         freeEnergyPerturbationData_ = std::make_unique<FreeEnergyPerturbationData>(
 422                 legacySimulatorData->fplog, *legacySimulatorData->inputrec, legacySimulatorData->mdAtoms);
 423         registerExistingElement(freeEnergyPerturbationData_->element());
 424     }
 425
 426     statePropagatorData_ = std::make_unique<StatePropagatorData>(
 427             legacySimulatorData->top_global.natoms,
 428             legacySimulatorData->fplog,
 429             legacySimulatorData->cr,
 430             legacySimulatorData->state_global,
 431             legacySimulatorData->state,
 432             legacySimulatorData->fr->nbv->useGpu(),
 433             legacySimulatorData->fr->bMolPBC,
 434             legacySimulatorData->mdrunOptions.writeConfout,
 435             opt2fn("-c", legacySimulatorData->nfile, legacySimulatorData->fnm),
 436             legacySimulatorData->inputrec,
 437             legacySimulatorData->mdAtoms->mdatoms(),
 438             legacySimulatorData->top_global);
 439     registerExistingElement(statePropagatorData_->element());
 440
 441     // Multi sim is turned off
 442     const bool simulationsShareState = false;
 443
 444     energyData_ = std::make_unique<EnergyData>(statePropagatorData_.get(),
 445                                                freeEnergyPerturbationData_.get(),
 446                                                legacySimulatorData->top_global,
 447                                                legacySimulatorData->inputrec,
 448                                                legacySimulatorData->mdAtoms,
 449                                                legacySimulatorData->enerd,
 450                                                legacySimulatorData->ekind,
 451                                                legacySimulatorData->constr,
 452                                                legacySimulatorData->fplog,
 453                                                legacySimulatorData->fr->fcdata.get(),
 454                                                legacySimulatorData->mdModulesNotifiers,
 455                                                MASTER(legacySimulatorData->cr),
 456                                                legacySimulatorData->observablesHistory,
 457                                                legacySimulatorData->startingBehavior,
 458                                                simulationsShareState,
 459                                                legacySimulatorData->pull_work);
 460     registerExistingElement(energyData_->element());
 461
 462     // This is the only modular simulator object which changes the inputrec
 463     // TODO: Avoid changing inputrec (#3854)
 464     storeSimulationData(
 465             "ReferenceTemperatureManager",
 466             ReferenceTemperatureManager(const_cast<t_inputrec*>(legacySimulatorData->inputrec)));
 467     auto* referenceTemperatureManager =
 468             simulationData<ReferenceTemperatureManager>("ReferenceTemperatureManager").value();
 469
 470     // State propagator data is scaling velocities if reference temperature is updated
 471     auto* statePropagatorDataPtr = statePropagatorData_.get();
 472     referenceTemperatureManager->registerUpdateCallback(
 473             [statePropagatorDataPtr](ArrayRef<const real>                temperatures,
 474                                      ReferenceTemperatureChangeAlgorithm algorithm) {
 475                 statePropagatorDataPtr->updateReferenceTemperature(temperatures, algorithm);
 476             });
 477 }
 478
 479 ModularSimulatorAlgorithm ModularSimulatorAlgorithmBuilder::build()
 480 {
 481     if (algorithmHasBeenBuilt_)
 482     {
 483         GMX_THROW(SimulationAlgorithmSetupError(
 484                 "Tried to build ModularSimulationAlgorithm more than once."));
 485     }
 486     algorithmHasBeenBuilt_ = true;
 487
 488     // Connect propagators with thermostat / barostat
 489     for (const auto& registrationFunction : pressureTemperatureControlRegistrationFunctions_)
 490     {
 491         for (const auto& connection : propagatorConnections_)
 492         {
 493             registrationFunction(connection);
 494         }
 495     }
 496
 497     ModularSimulatorAlgorithm algorithm(*(legacySimulatorData_->top_global.name),
 498                                         legacySimulatorData_->fplog,
 499                                         legacySimulatorData_->cr,
 500                                         legacySimulatorData_->mdlog,
 501                                         legacySimulatorData_->mdrunOptions,
 502                                         legacySimulatorData_->inputrec,
 503                                         legacySimulatorData_->nrnb,
 504                                         legacySimulatorData_->wcycle,
 505                                         legacySimulatorData_->fr,
 506                                         legacySimulatorData_->walltime_accounting);
 507     registerWithInfrastructureAndSignallers(algorithm.signalHelper_.get());
 508     algorithm.statePropagatorData_        = std::move(statePropagatorData_);
 509     algorithm.energyData_                 = std::move(energyData_);
 510     algorithm.freeEnergyPerturbationData_ = std::move(freeEnergyPerturbationData_);
 511     algorithm.signals_                    = std::move(signals_);
 512     algorithm.simulationData_             = std::move(simulationData_);
 513
 514     // Multi sim is turned off
 515     const bool simulationsShareState = false;
 516
 517     // Build stop handler
 518     algorithm.stopHandler_ = legacySimulatorData_->stopHandlerBuilder->getStopHandlerMD(
 519             compat::not_null<SimulationSignal*>(
 520                     &(*globalCommunicationHelper_.simulationSignals())[eglsSTOPCOND]),
 521             simulationsShareState,
 522             MASTER(legacySimulatorData_->cr),
 523             legacySimulatorData_->inputrec->nstlist,
 524             legacySimulatorData_->mdrunOptions.reproducible,
 525             globalCommunicationHelper_.nstglobalcomm(),
 526             legacySimulatorData_->mdrunOptions.maximumHoursToRun,
 527             legacySimulatorData_->inputrec->nstlist == 0,
 528             legacySimulatorData_->fplog,
 529             algorithm.stophandlerCurrentStep_,
 530             algorithm.stophandlerIsNSStep_,
 531             legacySimulatorData_->walltime_accounting);
 532
 533     // Build reset handler
 534     const bool simulationsShareResetCounters = false;
 535     algorithm.resetHandler_                  = std::make_unique<ResetHandler>(
 536             compat::make_not_null<SimulationSignal*>(
 537                     &(*globalCommunicationHelper_.simulationSignals())[eglsRESETCOUNTERS]),
 538             simulationsShareResetCounters,
 539             legacySimulatorData_->inputrec->nsteps,
 540             MASTER(legacySimulatorData_->cr),
 541             legacySimulatorData_->mdrunOptions.timingOptions.resetHalfway,
 542             legacySimulatorData_->mdrunOptions.maximumHoursToRun,
 543             legacySimulatorData_->mdlog,
 544             legacySimulatorData_->wcycle,
 545             legacySimulatorData_->walltime_accounting);
 546
 547     // Build topology holder
 548     algorithm.topologyHolder_ = topologyHolderBuilder_.build(legacySimulatorData_->top_global,
 549                                                              legacySimulatorData_->top,
 550                                                              legacySimulatorData_->cr,
 551                                                              legacySimulatorData_->inputrec,
 552                                                              legacySimulatorData_->fr,
 553                                                              legacySimulatorData_->mdAtoms,
 554                                                              legacySimulatorData_->constr,
 555                                                              legacySimulatorData_->vsite);
 556     registerWithInfrastructureAndSignallers(algorithm.topologyHolder_.get());
 557
 558     // Build PME load balance helper
 559     if (PmeLoadBalanceHelper::doPmeLoadBalancing(legacySimulatorData_->mdrunOptions,
 560                                                  legacySimulatorData_->inputrec,
 561                                                  legacySimulatorData_->fr))
 562     {
 563         algorithm.pmeLoadBalanceHelper_ =
 564                 std::make_unique<PmeLoadBalanceHelper>(legacySimulatorData_->mdrunOptions.verbose,
 565                                                        algorithm.statePropagatorData_.get(),
 566                                                        legacySimulatorData_->fplog,
 567                                                        legacySimulatorData_->cr,
 568                                                        legacySimulatorData_->mdlog,
 569                                                        legacySimulatorData_->inputrec,
 570                                                        legacySimulatorData_->wcycle,
 571                                                        legacySimulatorData_->fr);
 572         registerWithInfrastructureAndSignallers(algorithm.pmeLoadBalanceHelper_.get());
 573     }
 574
 575     // Build trajectory element
 576     auto trajectoryElement = trajectoryElementBuilder_.build(legacySimulatorData_->fplog,
 577                                                              legacySimulatorData_->nfile,
 578                                                              legacySimulatorData_->fnm,
 579                                                              legacySimulatorData_->mdrunOptions,
 580                                                              legacySimulatorData_->cr,
 581                                                              legacySimulatorData_->outputProvider,
 582                                                              legacySimulatorData_->mdModulesNotifiers,
 583                                                              legacySimulatorData_->inputrec,
 584                                                              legacySimulatorData_->top_global,
 585                                                              legacySimulatorData_->oenv,
 586                                                              legacySimulatorData_->wcycle,
 587                                                              legacySimulatorData_->startingBehavior,
 588                                                              simulationsShareState);
 589     registerWithInfrastructureAndSignallers(trajectoryElement.get());
 590
 591     // Build domdec helper (free energy element is a client, so keep this after it is built)
 592     if (haveDDAtomOrdering(*legacySimulatorData_->cr))
 593     {
 594         algorithm.domDecHelper_ =
 595                 domDecHelperBuilder_.build(legacySimulatorData_->mdrunOptions.verbose,
 596                                            legacySimulatorData_->mdrunOptions.verboseStepPrintInterval,
 597                                            algorithm.statePropagatorData_.get(),
 598                                            algorithm.topologyHolder_.get(),
 599                                            globalCommunicationHelper_.nstglobalcomm(),
 600                                            legacySimulatorData_->fplog,
 601                                            legacySimulatorData_->cr,
 602                                            legacySimulatorData_->mdlog,
 603                                            legacySimulatorData_->constr,
 604                                            legacySimulatorData_->inputrec,
 605                                            legacySimulatorData_->mdAtoms,
 606                                            legacySimulatorData_->nrnb,
 607                                            legacySimulatorData_->wcycle,
 608                                            legacySimulatorData_->fr,
 609                                            legacySimulatorData_->vsite,
 610                                            legacySimulatorData_->imdSession,
 611                                            legacySimulatorData_->pull_work);
 612         registerWithInfrastructureAndSignallers(algorithm.domDecHelper_.get());
 613     }
 614     // Build checkpoint helper (do this last so everyone else can be a checkpoint client!)
 615     {
 616         checkpointHelperBuilder_.setCheckpointHandler(std::make_unique<CheckpointHandler>(
 617                 compat::make_not_null<SimulationSignal*>(&(*algorithm.signals_)[eglsCHKPT]),
 618                 simulationsShareState,
 619                 legacySimulatorData_->inputrec->nstlist == 0,
 620                 MASTER(legacySimulatorData_->cr),
 621                 legacySimulatorData_->mdrunOptions.writeConfout,
 622                 legacySimulatorData_->mdrunOptions.checkpointOptions.period));
 623         algorithm.checkpointHelper_ =
 624                 checkpointHelperBuilder_.build(legacySimulatorData_->inputrec->init_step,
 625                                                trajectoryElement.get(),
 626                                                legacySimulatorData_->fplog,
 627                                                legacySimulatorData_->cr,
 628                                                legacySimulatorData_->observablesHistory,
 629                                                legacySimulatorData_->walltime_accounting,
 630                                                legacySimulatorData_->state_global,
 631                                                legacySimulatorData_->mdrunOptions.writeConfout);
 632         registerWithInfrastructureAndSignallers(algorithm.checkpointHelper_.get());
 633     }
 634
 635     // Build signallers
 636     {
 637         /* Signallers need to be called in an exact order. Some signallers are clients
 638          * of other signallers, which requires the clients signallers to be called
 639          * _after_ any signaller they are registered to - otherwise, they couldn't
 640          * adapt their behavior to the information they got signalled.
 641          *
 642          * Signallers being clients of other signallers require registration.
 643          * That registration happens during construction, which in turn means that
 644          * we want to construct the signallers in the reverse order of their later
 645          * call order.
 646          *
 647          * For the above reasons, the `addSignaller` lambda defined below emplaces
 648          * added signallers at the beginning of the signaller list, which will yield
 649          * a signaller list which is inverse to the build order (and hence equal to
 650          * the intended call order).
 651          */
 652         auto addSignaller = [this, &algorithm](auto signaller) {
 653             registerWithInfrastructureAndSignallers(signaller.get());
 654             algorithm.signallerList_.emplace(algorithm.signallerList_.begin(), std::move(signaller));
 655         };
 656         const auto* inputrec   = legacySimulatorData_->inputrec;
 657         auto        virialMode = EnergySignallerVirialMode::Off;
 658         if (inputrec->epc != PressureCoupling::No)
 659         {
 660             if (EI_VV(inputrec->eI))
 661             {
 662                 virialMode = EnergySignallerVirialMode::OnStepAndNext;
 663             }
 664             else
 665             {
 666                 virialMode = EnergySignallerVirialMode::OnStep;
 667             }
 668         }
 669         addSignaller(energySignallerBuilder_.build(
 670                 inputrec->nstcalcenergy,
 671                 computeFepPeriod(*inputrec, legacySimulatorData_->replExParams),
 672                 inputrec->nstpcouple,
 673                 virialMode));
 674         addSignaller(trajectorySignallerBuilder_.build(inputrec->nstxout,
 675                                                        inputrec->nstvout,
 676                                                        inputrec->nstfout,
 677                                                        inputrec->nstxout_compressed,
 678                                                        trajectoryElement->tngBoxOut(),
 679                                                        trajectoryElement->tngLambdaOut(),
 680                                                        trajectoryElement->tngBoxOutCompressed(),
 681                                                        trajectoryElement->tngLambdaOutCompressed(),
 682                                                        inputrec->nstenergy));
 683         addSignaller(loggingSignallerBuilder_.build(
 684                 inputrec->nstlog, inputrec->init_step, legacySimulatorData_->startingBehavior));
 685         addSignaller(lastStepSignallerBuilder_.build(
 686                 inputrec->nsteps, inputrec->init_step, algorithm.stopHandler_.get()));
 687         addSignaller(neighborSearchSignallerBuilder_.build(
 688                 inputrec->nstlist, inputrec->init_step, inputrec->init_t));
 689     }
 690
 691     // Move setup / teardown list
 692     algorithm.elementSetupTeardownList_ = std::move(setupAndTeardownList_);
 693     // Move pre- / post-step scheduling lists
 694     algorithm.preStepScheduling_  = std::move(preStepScheduling_);
 695     algorithm.postStepScheduling_ = std::move(postStepScheduling_);
 696
 697     // Create element list
 698     // Checkpoint helper needs to be in the call list (as first element!) to react to last step
 699     algorithm.elementCallList_.emplace_back(algorithm.checkpointHelper_.get());
 700     // Next, update the free energy lambda vector if needed
 701     if (algorithm.freeEnergyPerturbationData_)
 702     {
 703         algorithm.elementCallList_.emplace_back(algorithm.freeEnergyPerturbationData_->element());
 704     }
 705     // Then, move the built algorithm
 706     algorithm.elementsOwnershipList_.insert(algorithm.elementsOwnershipList_.end(),
 707                                             std::make_move_iterator(elements_.begin()),
 708                                             std::make_move_iterator(elements_.end()));
 709     algorithm.elementCallList_.insert(algorithm.elementCallList_.end(),
 710                                       std::make_move_iterator(callList_.begin()),
 711                                       std::make_move_iterator(callList_.end()));
 712     // Finally, all trajectory writing is happening after the step
 713     // (relevant data was stored by elements through energy signaller)
 714     algorithm.elementsOwnershipList_.emplace_back(std::move(trajectoryElement));
 715     algorithm.elementCallList_.emplace_back(algorithm.elementsOwnershipList_.back().get());
 716     algorithm.elementSetupTeardownList_.emplace_back(algorithm.elementsOwnershipList_.back().get());
 717
 718     algorithm.setup();
 719     return algorithm;
 720 }
 721
 722 bool ModularSimulatorAlgorithmBuilder::elementExists(const ISimulatorElement* element) const
 723 {
 724     // Check whether element exists in element list
 725     if (std::any_of(elements_.begin(), elements_.end(), [element](auto& existingElement) {
 726             return element == existingElement.get();
 727         }))
 728     {
 729         return true;
 730     }
 731     // Check whether element exists in other places controlled by *this
 732     return ((statePropagatorData_ && statePropagatorData_->element() == element)
 733             || (energyData_ && energyData_->element() == element)
 734             || (freeEnergyPerturbationData_ && freeEnergyPerturbationData_->element() == element));
 735 }
 736
 737 std::optional<SignallerCallback> ModularSimulatorAlgorithm::SignalHelper::registerLastStepCallback()
 738 {
 739     return [this](Step step, Time gmx_unused time) { this->lastStep_ = step; };
 740 }
 741
 742 std::optional<SignallerCallback> ModularSimulatorAlgorithm::SignalHelper::registerNSCallback()
 743 {
 744     return [this](Step step, Time gmx_unused time) { this->nextNSStep_ = step; };
 745 }
 746
 747 GlobalCommunicationHelper::GlobalCommunicationHelper(int nstglobalcomm, SimulationSignals* simulationSignals) :
 748     nstglobalcomm_(nstglobalcomm), simulationSignals_(simulationSignals)
 749 {
 750 }
 751
 752 int GlobalCommunicationHelper::nstglobalcomm() const
 753 {
 754     return nstglobalcomm_;
 755 }
 756
 757 SimulationSignals* GlobalCommunicationHelper::simulationSignals()
 758 {
 759     return simulationSignals_;
 760 }
 761
 762 ModularSimulatorAlgorithmBuilderHelper::ModularSimulatorAlgorithmBuilderHelper(
 763         ModularSimulatorAlgorithmBuilder* builder) :
 764     builder_(builder)
 765 {
 766 }
 767
 768 bool ModularSimulatorAlgorithmBuilderHelper::elementIsStored(const ISimulatorElement* element) const
 769 {
 770     return builder_->elementExists(element);
 771 }
 772
 773 [[maybe_unused]] void ModularSimulatorAlgorithmBuilderHelper::registerPreStepScheduling(SchedulingFunction schedulingFunction)
 774 {
 775     builder_->preStepScheduling_.emplace_back(std::move(schedulingFunction));
 776 }
 777
 778 [[maybe_unused]] void ModularSimulatorAlgorithmBuilderHelper::registerPostStepScheduling(SchedulingFunction schedulingFunction)
 779 {
 780     builder_->postStepScheduling_.emplace_back(std::move(schedulingFunction));
 781 }
 782
 783 std::optional<std::any> ModularSimulatorAlgorithmBuilderHelper::builderData(const std::string& key) const
 784 {
 785     const auto iter = builder_->builderData_.find(key);
 786     if (iter == builder_->builderData_.end())
 787     {
 788         return std::nullopt;
 789     }
 790     else
 791     {
 792         return iter->second;
 793     }
 794 }
 795
 796 void ModularSimulatorAlgorithmBuilderHelper::registerTemperaturePressureControl(
 797         std::function<void(const PropagatorConnection&)> registrationFunction)
 798 {
 799     builder_->pressureTemperatureControlRegistrationFunctions_.emplace_back(std::move(registrationFunction));
 800 }
 801
 802 void ModularSimulatorAlgorithmBuilderHelper::registerPropagator(PropagatorConnection connectionData)
 803 {
 804     builder_->propagatorConnections_.emplace_back(std::move(connectionData));
 805 }
 806
 807 void ModularSimulatorAlgorithmBuilderHelper::registerReferenceTemperatureUpdate(
 808         ReferenceTemperatureCallback referenceTemperatureCallback)
 809 {
 810     auto* referenceTemperatureManager =
 811             simulationData<ReferenceTemperatureManager>("ReferenceTemperatureManager").value();
 812     referenceTemperatureManager->registerUpdateCallback(std::move(referenceTemperatureCallback));
 813 }
 814
 815 ReferenceTemperatureCallback ModularSimulatorAlgorithmBuilderHelper::changeReferenceTemperatureCallback()
 816 {
 817     // Capture is safe because SimulatorAlgorithm will manage life time of both the
 818     // recipient of the callback and the reference temperature manager
 819     auto* referenceTemperatureManager =
 820             simulationData<ReferenceTemperatureManager>("ReferenceTemperatureManager").value();
 821     return [referenceTemperatureManager](ArrayRef<const real>                temperatures,
 822                                          ReferenceTemperatureChangeAlgorithm algorithm) {
 823         referenceTemperatureManager->setReferenceTemperature(temperatures, algorithm);
 824     };
 825 }
 826
 827 } // namespace gmx