src/gromacs/mdrun/runner.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2011-2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  *
  39  * \brief Implements the MD runner routine calling all integrators.
  40  *
  41  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  42  * \ingroup module_mdrun
  43  */
  44 #include "gmxpre.h"
  45
  46 #include "runner.h"
  47
  48 #include "config.h"
  49
  50 #include <cassert>
  51 #include <cinttypes>
  52 #include <csignal>
  53 #include <cstdlib>
  54 #include <cstring>
  55
  56 #include <algorithm>
  57 #include <memory>
  58
  59 #include "gromacs/commandline/filenm.h"
  60 #include "gromacs/domdec/builder.h"
  61 #include "gromacs/domdec/domdec.h"
  62 #include "gromacs/domdec/domdec_struct.h"
  63 #include "gromacs/domdec/gpuhaloexchange.h"
  64 #include "gromacs/domdec/localatomsetmanager.h"
  65 #include "gromacs/domdec/partition.h"
  66 #include "gromacs/ewald/ewald_utils.h"
  67 #include "gromacs/ewald/pme_gpu_program.h"
  68 #include "gromacs/ewald/pme_only.h"
  69 #include "gromacs/ewald/pme_pp_comm_gpu.h"
  70 #include "gromacs/fileio/checkpoint.h"
  71 #include "gromacs/fileio/gmxfio.h"
  72 #include "gromacs/fileio/oenv.h"
  73 #include "gromacs/fileio/tpxio.h"
  74 #include "gromacs/gmxlib/network.h"
  75 #include "gromacs/gmxlib/nrnb.h"
  76 #include "gromacs/gpu_utils/device_stream_manager.h"
  77 #include "gromacs/hardware/cpuinfo.h"
  78 #include "gromacs/hardware/detecthardware.h"
  79 #include "gromacs/hardware/device_management.h"
  80 #include "gromacs/hardware/hardwaretopology.h"
  81 #include "gromacs/hardware/printhardware.h"
  82 #include "gromacs/imd/imd.h"
  83 #include "gromacs/listed_forces/disre.h"
  84 #include "gromacs/listed_forces/gpubonded.h"
  85 #include "gromacs/listed_forces/listed_forces.h"
  86 #include "gromacs/listed_forces/orires.h"
  87 #include "gromacs/math/functions.h"
  88 #include "gromacs/math/utilities.h"
  89 #include "gromacs/math/vec.h"
  90 #include "gromacs/mdlib/boxdeformation.h"
  91 #include "gromacs/mdlib/broadcaststructs.h"
  92 #include "gromacs/mdlib/calc_verletbuf.h"
  93 #include "gromacs/mdlib/dispersioncorrection.h"
  94 #include "gromacs/mdlib/enerdata_utils.h"
  95 #include "gromacs/mdlib/force.h"
  96 #include "gromacs/mdlib/forcerec.h"
  97 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  98 #include "gromacs/mdlib/gpuforcereduction.h"
  99 #include "gromacs/mdlib/makeconstraints.h"
 100 #include "gromacs/mdlib/md_support.h"
 101 #include "gromacs/mdlib/mdatoms.h"
 102 #include "gromacs/mdlib/sighandler.h"
 103 #include "gromacs/mdlib/stophandler.h"
 104 #include "gromacs/mdlib/tgroup.h"
 105 #include "gromacs/mdlib/updategroups.h"
 106 #include "gromacs/mdlib/vsite.h"
 107 #include "gromacs/mdrun/mdmodules.h"
 108 #include "gromacs/mdrun/simulationcontext.h"
 109 #include "gromacs/mdrun/simulationinput.h"
 110 #include "gromacs/mdrun/simulationinputhandle.h"
 111 #include "gromacs/mdrunutility/handlerestart.h"
 112 #include "gromacs/mdrunutility/logging.h"
 113 #include "gromacs/mdrunutility/multisim.h"
 114 #include "gromacs/mdrunutility/printtime.h"
 115 #include "gromacs/mdrunutility/threadaffinity.h"
 116 #include "gromacs/mdtypes/checkpointdata.h"
 117 #include "gromacs/mdtypes/commrec.h"
 118 #include "gromacs/mdtypes/enerdata.h"
 119 #include "gromacs/mdtypes/fcdata.h"
 120 #include "gromacs/mdtypes/forcerec.h"
 121 #include "gromacs/mdtypes/group.h"
 122 #include "gromacs/mdtypes/inputrec.h"
 123 #include "gromacs/mdtypes/interaction_const.h"
 124 #include "gromacs/mdtypes/md_enums.h"
 125 #include "gromacs/mdtypes/mdatom.h"
 126 #include "gromacs/mdtypes/mdrunoptions.h"
 127 #include "gromacs/mdtypes/observableshistory.h"
 128 #include "gromacs/mdtypes/simulation_workload.h"
 129 #include "gromacs/mdtypes/state.h"
 130 #include "gromacs/mdtypes/state_propagator_data_gpu.h"
 131 #include "gromacs/modularsimulator/modularsimulator.h"
 132 #include "gromacs/nbnxm/gpu_data_mgmt.h"
 133 #include "gromacs/nbnxm/nbnxm.h"
 134 #include "gromacs/nbnxm/pairlist_tuning.h"
 135 #include "gromacs/pbcutil/pbc.h"
 136 #include "gromacs/pulling/output.h"
 137 #include "gromacs/pulling/pull.h"
 138 #include "gromacs/pulling/pull_rotation.h"
 139 #include "gromacs/restraint/manager.h"
 140 #include "gromacs/restraint/restraintmdmodule.h"
 141 #include "gromacs/restraint/restraintpotential.h"
 142 #include "gromacs/swap/swapcoords.h"
 143 #include "gromacs/taskassignment/decidegpuusage.h"
 144 #include "gromacs/taskassignment/decidesimulationworkload.h"
 145 #include "gromacs/taskassignment/resourcedivision.h"
 146 #include "gromacs/taskassignment/taskassignment.h"
 147 #include "gromacs/taskassignment/usergpuids.h"
 148 #include "gromacs/timing/gpu_timing.h"
 149 #include "gromacs/timing/wallcycle.h"
 150 #include "gromacs/timing/wallcyclereporting.h"
 151 #include "gromacs/topology/mtop_util.h"
 152 #include "gromacs/trajectory/trajectoryframe.h"
 153 #include "gromacs/utility/basenetwork.h"
 154 #include "gromacs/utility/cstringutil.h"
 155 #include "gromacs/utility/exceptions.h"
 156 #include "gromacs/utility/fatalerror.h"
 157 #include "gromacs/utility/filestream.h"
 158 #include "gromacs/utility/gmxassert.h"
 159 #include "gromacs/utility/gmxmpi.h"
 160 #include "gromacs/utility/keyvaluetree.h"
 161 #include "gromacs/utility/logger.h"
 162 #include "gromacs/utility/loggerbuilder.h"
 163 #include "gromacs/utility/mdmodulenotification.h"
 164 #include "gromacs/utility/physicalnodecommunicator.h"
 165 #include "gromacs/utility/pleasecite.h"
 166 #include "gromacs/utility/programcontext.h"
 167 #include "gromacs/utility/smalloc.h"
 168 #include "gromacs/utility/stringutil.h"
 169
 170 #include "isimulator.h"
 171 #include "membedholder.h"
 172 #include "replicaexchange.h"
 173 #include "simulatorbuilder.h"
 174
 175 namespace gmx
 176 {
 177
 178
 179 /*! \brief Manage any development feature flag variables encountered
 180  *
 181  * The use of dev features indicated by environment variables is
 182  * logged in order to ensure that runs with such features enabled can
 183  * be identified from their log and standard output. Any cross
 184  * dependencies are also checked, and if unsatisfied, a fatal error
 185  * issued.
 186  *
 187  * Note that some development features overrides are applied already here:
 188  * the GPU communication flags are set to false in non-tMPI and non-CUDA builds.
 189  *
 190  * \param[in]  mdlog                Logger object.
 191  * \param[in]  useGpuForNonbonded   True if the nonbonded task is offloaded in this run.
 192  * \param[in]  pmeRunMode           The PME run mode for this run
 193  * \returns                         The object populated with development feature flags.
 194  */
 195 static DevelopmentFeatureFlags manageDevelopmentFeatures(const gmx::MDLogger& mdlog,
 196                                                          const bool           useGpuForNonbonded,
 197                                                          const PmeRunMode     pmeRunMode)
 198 {
 199     DevelopmentFeatureFlags devFlags;
 200
 201     // Some builds of GCC 5 give false positive warnings that these
 202     // getenv results are ignored when clearly they are used.
 203 #pragma GCC diagnostic push
 204 #pragma GCC diagnostic ignored "-Wunused-result"
 205
 206     devFlags.enableGpuBufferOps =
 207             GMX_GPU_CUDA && useGpuForNonbonded && (getenv("GMX_USE_GPU_BUFFER_OPS") != nullptr);
 208     devFlags.enableGpuHaloExchange = GMX_GPU_CUDA && GMX_THREAD_MPI && getenv("GMX_GPU_DD_COMMS") != nullptr;
 209     devFlags.forceGpuUpdateDefault = (getenv("GMX_FORCE_UPDATE_DEFAULT_GPU") != nullptr) || GMX_FAHCORE;
 210     devFlags.enableGpuPmePPComm =
 211             GMX_GPU_CUDA && GMX_THREAD_MPI && getenv("GMX_GPU_PME_PP_COMMS") != nullptr;
 212
 213 #pragma GCC diagnostic pop
 214
 215     if (devFlags.enableGpuBufferOps)
 216     {
 217         GMX_LOG(mdlog.warning)
 218                 .asParagraph()
 219                 .appendTextFormatted(
 220                         "This run uses the 'GPU buffer ops' feature, enabled by the "
 221                         "GMX_USE_GPU_BUFFER_OPS environment variable.");
 222     }
 223
 224     if (devFlags.forceGpuUpdateDefault)
 225     {
 226         GMX_LOG(mdlog.warning)
 227                 .asParagraph()
 228                 .appendTextFormatted(
 229                         "This run will default to '-update gpu' as requested by the "
 230                         "GMX_FORCE_UPDATE_DEFAULT_GPU environment variable. GPU update with domain "
 231                         "decomposition lacks substantial testing and should be used with caution.");
 232     }
 233
 234     if (devFlags.enableGpuHaloExchange)
 235     {
 236         if (useGpuForNonbonded)
 237         {
 238             if (!devFlags.enableGpuBufferOps)
 239             {
 240                 GMX_LOG(mdlog.warning)
 241                         .asParagraph()
 242                         .appendTextFormatted(
 243                                 "Enabling GPU buffer operations required by GMX_GPU_DD_COMMS "
 244                                 "(equivalent with GMX_USE_GPU_BUFFER_OPS=1).");
 245                 devFlags.enableGpuBufferOps = true;
 246             }
 247             GMX_LOG(mdlog.warning)
 248                     .asParagraph()
 249                     .appendTextFormatted(
 250                             "This run has requested the 'GPU halo exchange' feature, enabled by "
 251                             "the "
 252                             "GMX_GPU_DD_COMMS environment variable.");
 253         }
 254         else
 255         {
 256             GMX_LOG(mdlog.warning)
 257                     .asParagraph()
 258                     .appendTextFormatted(
 259                             "GMX_GPU_DD_COMMS environment variable detected, but the 'GPU "
 260                             "halo exchange' feature will not be enabled as nonbonded interactions "
 261                             "are not offloaded.");
 262             devFlags.enableGpuHaloExchange = false;
 263         }
 264     }
 265
 266     if (devFlags.enableGpuPmePPComm)
 267     {
 268         if (pmeRunMode == PmeRunMode::GPU)
 269         {
 270             if (!devFlags.enableGpuBufferOps)
 271             {
 272                 GMX_LOG(mdlog.warning)
 273                         .asParagraph()
 274                         .appendTextFormatted(
 275                                 "Enabling GPU buffer operations required by GMX_GPU_PME_PP_COMMS "
 276                                 "(equivalent with GMX_USE_GPU_BUFFER_OPS=1).");
 277                 devFlags.enableGpuBufferOps = true;
 278             }
 279             GMX_LOG(mdlog.warning)
 280                     .asParagraph()
 281                     .appendTextFormatted(
 282                             "This run uses the 'GPU PME-PP communications' feature, enabled "
 283                             "by the GMX_GPU_PME_PP_COMMS environment variable.");
 284         }
 285         else
 286         {
 287             std::string clarification;
 288             if (pmeRunMode == PmeRunMode::Mixed)
 289             {
 290                 clarification =
 291                         "PME FFT and gather are not offloaded to the GPU (PME is running in mixed "
 292                         "mode).";
 293             }
 294             else
 295             {
 296                 clarification = "PME is not offloaded to the GPU.";
 297             }
 298             GMX_LOG(mdlog.warning)
 299                     .asParagraph()
 300                     .appendText(
 301                             "GMX_GPU_PME_PP_COMMS environment variable detected, but the "
 302                             "'GPU PME-PP communications' feature was not enabled as "
 303                             + clarification);
 304             devFlags.enableGpuPmePPComm = false;
 305         }
 306     }
 307
 308     return devFlags;
 309 }
 310
 311 /*! \brief Barrier for safe simultaneous thread access to mdrunner data
 312  *
 313  * Used to ensure that the master thread does not modify mdrunner during copy
 314  * on the spawned threads. */
 315 static void threadMpiMdrunnerAccessBarrier()
 316 {
 317 #if GMX_THREAD_MPI
 318     MPI_Barrier(MPI_COMM_WORLD);
 319 #endif
 320 }
 321
 322 Mdrunner Mdrunner::cloneOnSpawnedThread() const
 323 {
 324     auto newRunner = Mdrunner(std::make_unique<MDModules>());
 325
 326     // All runners in the same process share a restraint manager resource because it is
 327     // part of the interface to the client code, which is associated only with the
 328     // original thread. Handles to the same resources can be obtained by copy.
 329     {
 330         newRunner.restraintManager_ = std::make_unique<RestraintManager>(*restraintManager_);
 331     }
 332
 333     // Copy members of master runner.
 334     // \todo Replace with builder when Simulation context and/or runner phases are better defined.
 335     // Ref https://gitlab.com/gromacs/gromacs/-/issues/2587 and https://gitlab.com/gromacs/gromacs/-/issues/2375
 336     newRunner.hw_opt    = hw_opt;
 337     newRunner.filenames = filenames;
 338
 339     newRunner.hwinfo_         = hwinfo_;
 340     newRunner.oenv            = oenv;
 341     newRunner.mdrunOptions    = mdrunOptions;
 342     newRunner.domdecOptions   = domdecOptions;
 343     newRunner.nbpu_opt        = nbpu_opt;
 344     newRunner.pme_opt         = pme_opt;
 345     newRunner.pme_fft_opt     = pme_fft_opt;
 346     newRunner.bonded_opt      = bonded_opt;
 347     newRunner.update_opt      = update_opt;
 348     newRunner.nstlist_cmdline = nstlist_cmdline;
 349     newRunner.replExParams    = replExParams;
 350     newRunner.pforce          = pforce;
 351     // Give the spawned thread the newly created valid communicator
 352     // for the simulation.
 353     newRunner.libraryWorldCommunicator = MPI_COMM_WORLD;
 354     newRunner.simulationCommunicator   = MPI_COMM_WORLD;
 355     newRunner.ms                       = ms;
 356     newRunner.startingBehavior         = startingBehavior;
 357     newRunner.stopHandlerBuilder_      = std::make_unique<StopHandlerBuilder>(*stopHandlerBuilder_);
 358     newRunner.inputHolder_             = inputHolder_;
 359
 360     threadMpiMdrunnerAccessBarrier();
 361
 362     return newRunner;
 363 }
 364
 365 /*! \brief The callback used for running on spawned threads.
 366  *
 367  * Obtains the pointer to the master mdrunner object from the one
 368  * argument permitted to the thread-launch API call, copies it to make
 369  * a new runner for this thread, reinitializes necessary data, and
 370  * proceeds to the simulation. */
 371 static void mdrunner_start_fn(const void* arg)
 372 {
 373     try
 374     {
 375         auto masterMdrunner = reinterpret_cast<const gmx::Mdrunner*>(arg);
 376         /* copy the arg list to make sure that it's thread-local. This
 377            doesn't copy pointed-to items, of course; fnm, cr and fplog
 378            are reset in the call below, all others should be const. */
 379         gmx::Mdrunner mdrunner = masterMdrunner->cloneOnSpawnedThread();
 380         mdrunner.mdrunner();
 381     }
 382     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 383 }
 384
 385
 386 void Mdrunner::spawnThreads(int numThreadsToLaunch)
 387 {
 388 #if GMX_THREAD_MPI
 389     /* now spawn new threads that start mdrunner_start_fn(), while
 390        the main thread returns. Thread affinity is handled later. */
 391     if (tMPI_Init_fn(TRUE, numThreadsToLaunch, TMPI_AFFINITY_NONE, mdrunner_start_fn,
 392                      static_cast<const void*>(this))
 393         != TMPI_SUCCESS)
 394     {
 395         GMX_THROW(gmx::InternalError("Failed to spawn thread-MPI threads"));
 396     }
 397
 398     // Give the master thread the newly created valid communicator for
 399     // the simulation.
 400     libraryWorldCommunicator = MPI_COMM_WORLD;
 401     simulationCommunicator   = MPI_COMM_WORLD;
 402     threadMpiMdrunnerAccessBarrier();
 403 #else
 404     GMX_UNUSED_VALUE(numThreadsToLaunch);
 405     GMX_UNUSED_VALUE(mdrunner_start_fn);
 406 #endif
 407 }
 408
 409 } // namespace gmx
 410
 411 /*! \brief Initialize variables for Verlet scheme simulation */
 412 static void prepare_verlet_scheme(FILE*               fplog,
 413                                   t_commrec*          cr,
 414                                   t_inputrec*         ir,
 415                                   int                 nstlist_cmdline,
 416                                   const gmx_mtop_t*   mtop,
 417                                   const matrix        box,
 418                                   bool                makeGpuPairList,
 419                                   const gmx::CpuInfo& cpuinfo)
 420 {
 421     // We checked the cut-offs in grompp, but double-check here.
 422     // We have PME+LJcutoff kernels for rcoulomb>rvdw.
 423     if (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == eelCUT)
 424     {
 425         GMX_RELEASE_ASSERT(ir->rcoulomb >= ir->rvdw,
 426                            "With Verlet lists and PME we should have rcoulomb>=rvdw");
 427     }
 428     else
 429     {
 430         GMX_RELEASE_ASSERT(ir->rcoulomb == ir->rvdw,
 431                            "With Verlet lists and no PME rcoulomb and rvdw should be identical");
 432     }
 433     /* For NVE simulations, we will retain the initial list buffer */
 434     if (EI_DYNAMICS(ir->eI) && ir->verletbuf_tol > 0 && !(EI_MD(ir->eI) && ir->etc == etcNO))
 435     {
 436         /* Update the Verlet buffer size for the current run setup */
 437
 438         /* Here we assume SIMD-enabled kernels are being used. But as currently
 439          * calc_verlet_buffer_size gives the same results for 4x8 and 4x4
 440          * and 4x2 gives a larger buffer than 4x4, this is ok.
 441          */
 442         ListSetupType listType =
 443                 (makeGpuPairList ? ListSetupType::Gpu : ListSetupType::CpuSimdWhenSupported);
 444         VerletbufListSetup listSetup = verletbufGetSafeListSetup(listType);
 445
 446         const real rlist_new =
 447                 calcVerletBufferSize(*mtop, det(box), *ir, ir->nstlist, ir->nstlist - 1, -1, listSetup);
 448
 449         if (rlist_new != ir->rlist)
 450         {
 451             if (fplog != nullptr)
 452             {
 453                 fprintf(fplog,
 454                         "\nChanging rlist from %g to %g for non-bonded %dx%d atom kernels\n\n",
 455                         ir->rlist, rlist_new, listSetup.cluster_size_i, listSetup.cluster_size_j);
 456             }
 457             ir->rlist = rlist_new;
 458         }
 459     }
 460
 461     if (nstlist_cmdline > 0 && (!EI_DYNAMICS(ir->eI) || ir->verletbuf_tol <= 0))
 462     {
 463         gmx_fatal(FARGS, "Can not set nstlist without %s",
 464                   !EI_DYNAMICS(ir->eI) ? "dynamics" : "verlet-buffer-tolerance");
 465     }
 466
 467     if (EI_DYNAMICS(ir->eI))
 468     {
 469         /* Set or try nstlist values */
 470         increaseNstlist(fplog, cr, ir, nstlist_cmdline, mtop, box, makeGpuPairList, cpuinfo);
 471     }
 472 }
 473
 474 /*! \brief Override the nslist value in inputrec
 475  *
 476  * with value passed on the command line (if any)
 477  */
 478 static void override_nsteps_cmdline(const gmx::MDLogger& mdlog, int64_t nsteps_cmdline, t_inputrec* ir)
 479 {
 480     assert(ir);
 481
 482     /* override with anything else than the default -2 */
 483     if (nsteps_cmdline > -2)
 484     {
 485         char sbuf_steps[STEPSTRSIZE];
 486         char sbuf_msg[STRLEN];
 487
 488         ir->nsteps = nsteps_cmdline;
 489         if (EI_DYNAMICS(ir->eI) && nsteps_cmdline != -1)
 490         {
 491             sprintf(sbuf_msg,
 492                     "Overriding nsteps with value passed on the command line: %s steps, %.3g ps",
 493                     gmx_step_str(nsteps_cmdline, sbuf_steps), fabs(nsteps_cmdline * ir->delta_t));
 494         }
 495         else
 496         {
 497             sprintf(sbuf_msg, "Overriding nsteps with value passed on the command line: %s steps",
 498                     gmx_step_str(nsteps_cmdline, sbuf_steps));
 499         }
 500
 501         GMX_LOG(mdlog.warning).asParagraph().appendText(sbuf_msg);
 502     }
 503     else if (nsteps_cmdline < -2)
 504     {
 505         gmx_fatal(FARGS, "Invalid nsteps value passed on the command line: %" PRId64, nsteps_cmdline);
 506     }
 507     /* Do nothing if nsteps_cmdline == -2 */
 508 }
 509
 510 namespace gmx
 511 {
 512
 513 /*! \brief Return whether GPU acceleration of nonbondeds is supported with the given settings.
 514  *
 515  * If not, and if a warning may be issued, logs a warning about
 516  * falling back to CPU code. With thread-MPI, only the first
 517  * call to this function should have \c issueWarning true. */
 518 static bool gpuAccelerationOfNonbondedIsUseful(const MDLogger& mdlog, const t_inputrec& ir, bool issueWarning)
 519 {
 520     bool        gpuIsUseful = true;
 521     std::string warning;
 522
 523     if (ir.opts.ngener - ir.nwall > 1)
 524     {
 525         /* The GPU code does not support more than one energy group.
 526          * If the user requested GPUs explicitly, a fatal error is given later.
 527          */
 528         gpuIsUseful = false;
 529         warning =
 530                 "Multiple energy groups is not implemented for GPUs, falling back to the CPU. "
 531                 "For better performance, run on the GPU without energy groups and then do "
 532                 "gmx mdrun -rerun option on the trajectory with an energy group .tpr file.";
 533     }
 534
 535     if (EI_TPI(ir.eI))
 536     {
 537         gpuIsUseful = false;
 538         warning     = "TPI is not implemented for GPUs.";
 539     }
 540
 541     if (!gpuIsUseful && issueWarning)
 542     {
 543         GMX_LOG(mdlog.warning).asParagraph().appendText(warning);
 544     }
 545
 546     return gpuIsUseful;
 547 }
 548
 549 //! Initializes the logger for mdrun.
 550 static gmx::LoggerOwner buildLogger(FILE* fplog, const bool isSimulationMasterRank)
 551 {
 552     gmx::LoggerBuilder builder;
 553     if (fplog != nullptr)
 554     {
 555         builder.addTargetFile(gmx::MDLogger::LogLevel::Info, fplog);
 556     }
 557     if (isSimulationMasterRank)
 558     {
 559         builder.addTargetStream(gmx::MDLogger::LogLevel::Warning, &gmx::TextOutputFile::standardError());
 560     }
 561     return builder.build();
 562 }
 563
 564 //! Make a TaskTarget from an mdrun argument string.
 565 static TaskTarget findTaskTarget(const char* optionString)
 566 {
 567     TaskTarget returnValue = TaskTarget::Auto;
 568
 569     if (strncmp(optionString, "auto", 3) == 0)
 570     {
 571         returnValue = TaskTarget::Auto;
 572     }
 573     else if (strncmp(optionString, "cpu", 3) == 0)
 574     {
 575         returnValue = TaskTarget::Cpu;
 576     }
 577     else if (strncmp(optionString, "gpu", 3) == 0)
 578     {
 579         returnValue = TaskTarget::Gpu;
 580     }
 581     else
 582     {
 583         GMX_ASSERT(false, "Option string should have been checked for sanity already");
 584     }
 585
 586     return returnValue;
 587 }
 588
 589 //! Finish run, aggregate data to print performance info.
 590 static void finish_run(FILE*                     fplog,
 591                        const gmx::MDLogger&      mdlog,
 592                        const t_commrec*          cr,
 593                        const t_inputrec*         inputrec,
 594                        t_nrnb                    nrnb[],
 595                        gmx_wallcycle_t           wcycle,
 596                        gmx_walltime_accounting_t walltime_accounting,
 597                        nonbonded_verlet_t*       nbv,
 598                        const gmx_pme_t*          pme,
 599                        gmx_bool                  bWriteStat)
 600 {
 601     double delta_t = 0;
 602     double nbfs = 0, mflop = 0;
 603     double elapsed_time, elapsed_time_over_all_ranks, elapsed_time_over_all_threads,
 604             elapsed_time_over_all_threads_over_all_ranks;
 605     /* Control whether it is valid to print a report. Only the
 606        simulation master may print, but it should not do so if the run
 607        terminated e.g. before a scheduled reset step. This is
 608        complicated by the fact that PME ranks are unaware of the
 609        reason why they were sent a pmerecvqxFINISH. To avoid
 610        communication deadlocks, we always do the communication for the
 611        report, even if we've decided not to write the report, because
 612        how long it takes to finish the run is not important when we've
 613        decided not to report on the simulation performance.
 614
 615        Further, we only report performance for dynamical integrators,
 616        because those are the only ones for which we plan to
 617        consider doing any optimizations. */
 618     bool printReport = EI_DYNAMICS(inputrec->eI) && SIMMASTER(cr);
 619
 620     if (printReport && !walltime_accounting_get_valid_finish(walltime_accounting))
 621     {
 622         GMX_LOG(mdlog.warning)
 623                 .asParagraph()
 624                 .appendText("Simulation ended prematurely, no performance report will be written.");
 625         printReport = false;
 626     }
 627
 628     t_nrnb*                 nrnb_tot;
 629     std::unique_ptr<t_nrnb> nrnbTotalStorage;
 630     if (cr->nnodes > 1)
 631     {
 632         nrnbTotalStorage = std::make_unique<t_nrnb>();
 633         nrnb_tot         = nrnbTotalStorage.get();
 634 #if GMX_MPI
 635         MPI_Allreduce(nrnb->n, nrnb_tot->n, eNRNB, MPI_DOUBLE, MPI_SUM, cr->mpi_comm_mysim);
 636 #endif
 637     }
 638     else
 639     {
 640         nrnb_tot = nrnb;
 641     }
 642
 643     elapsed_time = walltime_accounting_get_time_since_reset(walltime_accounting);
 644     elapsed_time_over_all_threads =
 645             walltime_accounting_get_time_since_reset_over_all_threads(walltime_accounting);
 646     if (cr->nnodes > 1)
 647     {
 648 #if GMX_MPI
 649         /* reduce elapsed_time over all MPI ranks in the current simulation */
 650         MPI_Allreduce(&elapsed_time, &elapsed_time_over_all_ranks, 1, MPI_DOUBLE, MPI_SUM,
 651                       cr->mpi_comm_mysim);
 652         elapsed_time_over_all_ranks /= cr->nnodes;
 653         /* Reduce elapsed_time_over_all_threads over all MPI ranks in the
 654          * current simulation. */
 655         MPI_Allreduce(&elapsed_time_over_all_threads, &elapsed_time_over_all_threads_over_all_ranks,
 656                       1, MPI_DOUBLE, MPI_SUM, cr->mpi_comm_mysim);
 657 #endif
 658     }
 659     else
 660     {
 661         elapsed_time_over_all_ranks                  = elapsed_time;
 662         elapsed_time_over_all_threads_over_all_ranks = elapsed_time_over_all_threads;
 663     }
 664
 665     if (printReport)
 666     {
 667         print_flop(fplog, nrnb_tot, &nbfs, &mflop);
 668     }
 669
 670     if (thisRankHasDuty(cr, DUTY_PP) && DOMAINDECOMP(cr))
 671     {
 672         print_dd_statistics(cr, inputrec, fplog);
 673     }
 674
 675     /* TODO Move the responsibility for any scaling by thread counts
 676      * to the code that handled the thread region, so that there's a
 677      * mechanism to keep cycle counting working during the transition
 678      * to task parallelism. */
 679     int nthreads_pp  = gmx_omp_nthreads_get(emntNonbonded);
 680     int nthreads_pme = gmx_omp_nthreads_get(emntPME);
 681     wallcycle_scale_by_num_threads(wcycle, thisRankHasDuty(cr, DUTY_PME) && !thisRankHasDuty(cr, DUTY_PP),
 682                                    nthreads_pp, nthreads_pme);
 683     auto cycle_sum(wallcycle_sum(cr, wcycle));
 684
 685     if (printReport)
 686     {
 687         auto nbnxn_gpu_timings =
 688                 (nbv != nullptr && nbv->useGpu()) ? Nbnxm::gpu_get_timings(nbv->gpu_nbv) : nullptr;
 689         gmx_wallclock_gpu_pme_t pme_gpu_timings = {};
 690
 691         if (pme_gpu_task_enabled(pme))
 692         {
 693             pme_gpu_get_timings(pme, &pme_gpu_timings);
 694         }
 695         wallcycle_print(fplog, mdlog, cr->nnodes, cr->npmenodes, nthreads_pp, nthreads_pme,
 696                         elapsed_time_over_all_ranks, wcycle, cycle_sum, nbnxn_gpu_timings,
 697                         &pme_gpu_timings);
 698
 699         if (EI_DYNAMICS(inputrec->eI))
 700         {
 701             delta_t = inputrec->delta_t;
 702         }
 703
 704         if (fplog)
 705         {
 706             print_perf(fplog, elapsed_time_over_all_threads_over_all_ranks, elapsed_time_over_all_ranks,
 707                        walltime_accounting_get_nsteps_done_since_reset(walltime_accounting),
 708                        delta_t, nbfs, mflop);
 709         }
 710         if (bWriteStat)
 711         {
 712             print_perf(stderr, elapsed_time_over_all_threads_over_all_ranks, elapsed_time_over_all_ranks,
 713                        walltime_accounting_get_nsteps_done_since_reset(walltime_accounting),
 714                        delta_t, nbfs, mflop);
 715         }
 716     }
 717 }
 718
 719 int Mdrunner::mdrunner()
 720 {
 721     matrix                    box;
 722     t_forcerec*               fr               = nullptr;
 723     real                      ewaldcoeff_q     = 0;
 724     real                      ewaldcoeff_lj    = 0;
 725     int                       nChargePerturbed = -1, nTypePerturbed = 0;
 726     gmx_wallcycle_t           wcycle;
 727     gmx_walltime_accounting_t walltime_accounting = nullptr;
 728     MembedHolder              membedHolder(filenames.size(), filenames.data());
 729
 730     /* CAUTION: threads may be started later on in this function, so
 731        cr doesn't reflect the final parallel state right now */
 732     gmx_mtop_t mtop;
 733
 734     /* TODO: inputrec should tell us whether we use an algorithm, not a file option */
 735     const bool doEssentialDynamics = opt2bSet("-ei", filenames.size(), filenames.data());
 736     const bool doRerun             = mdrunOptions.rerun;
 737
 738     // Handle task-assignment related user options.
 739     EmulateGpuNonbonded emulateGpuNonbonded =
 740             (getenv("GMX_EMULATE_GPU") != nullptr ? EmulateGpuNonbonded::Yes : EmulateGpuNonbonded::No);
 741
 742     std::vector<int> userGpuTaskAssignment;
 743     try
 744     {
 745         userGpuTaskAssignment = parseUserTaskAssignmentString(hw_opt.userGpuTaskAssignment);
 746     }
 747     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 748     auto nonbondedTarget = findTaskTarget(nbpu_opt);
 749     auto pmeTarget       = findTaskTarget(pme_opt);
 750     auto pmeFftTarget    = findTaskTarget(pme_fft_opt);
 751     auto bondedTarget    = findTaskTarget(bonded_opt);
 752     auto updateTarget    = findTaskTarget(update_opt);
 753
 754     FILE* fplog = nullptr;
 755     // If we are appending, we don't write log output because we need
 756     // to check that the old log file matches what the checkpoint file
 757     // expects. Otherwise, we should start to write log output now if
 758     // there is a file ready for it.
 759     if (logFileHandle != nullptr && startingBehavior != StartingBehavior::RestartWithAppending)
 760     {
 761         fplog = gmx_fio_getfp(logFileHandle);
 762     }
 763     const bool isSimulationMasterRank = findIsSimulationMasterRank(ms, simulationCommunicator);
 764     gmx::LoggerOwner logOwner(buildLogger(fplog, isSimulationMasterRank));
 765     gmx::MDLogger    mdlog(logOwner.logger());
 766
 767     gmx_print_detected_hardware(fplog, isSimulationMasterRank && isMasterSim(ms), mdlog, hwinfo_);
 768
 769     std::vector<int> gpuIdsToUse = makeGpuIdsToUse(hwinfo_->deviceInfoList, hw_opt.gpuIdsAvailable);
 770     const int        numDevicesToUse = gmx::ssize(gpuIdsToUse);
 771
 772     // Print citation requests after all software/hardware printing
 773     pleaseCiteGromacs(fplog);
 774
 775     // Note: legacy program logic relies on checking whether these pointers are assigned.
 776     // Objects may or may not be allocated later.
 777     std::unique_ptr<t_inputrec> inputrec;
 778     std::unique_ptr<t_state>    globalState;
 779
 780     auto partialDeserializedTpr = std::make_unique<PartialDeserializedTprFile>();
 781
 782     if (isSimulationMasterRank)
 783     {
 784         // Allocate objects to be initialized by later function calls.
 785         /* Only the master rank has the global state */
 786         globalState = std::make_unique<t_state>();
 787         inputrec    = std::make_unique<t_inputrec>();
 788
 789         /* Read (nearly) all data required for the simulation
 790          * and keep the partly serialized tpr contents to send to other ranks later
 791          */
 792         applyGlobalSimulationState(*inputHolder_.get(), partialDeserializedTpr.get(),
 793                                    globalState.get(), inputrec.get(), &mtop);
 794     }
 795
 796     /* Check and update the hardware options for internal consistency */
 797     checkAndUpdateHardwareOptions(mdlog, &hw_opt, isSimulationMasterRank, domdecOptions.numPmeRanks,
 798                                   inputrec.get());
 799
 800     if (GMX_THREAD_MPI && isSimulationMasterRank)
 801     {
 802         bool useGpuForNonbonded = false;
 803         bool useGpuForPme       = false;
 804         try
 805         {
 806             GMX_RELEASE_ASSERT(inputrec != nullptr, "Keep the compiler happy");
 807
 808             // If the user specified the number of ranks, then we must
 809             // respect that, but in default mode, we need to allow for
 810             // the number of GPUs to choose the number of ranks.
 811             auto canUseGpuForNonbonded = buildSupportsNonbondedOnGpu(nullptr);
 812             useGpuForNonbonded         = decideWhetherToUseGpusForNonbondedWithThreadMpi(
 813                     nonbondedTarget, numDevicesToUse, userGpuTaskAssignment, emulateGpuNonbonded,
 814                     canUseGpuForNonbonded,
 815                     gpuAccelerationOfNonbondedIsUseful(mdlog, *inputrec, GMX_THREAD_MPI),
 816                     hw_opt.nthreads_tmpi);
 817             useGpuForPme = decideWhetherToUseGpusForPmeWithThreadMpi(
 818                     useGpuForNonbonded, pmeTarget, numDevicesToUse, userGpuTaskAssignment, *hwinfo_,
 819                     *inputrec, hw_opt.nthreads_tmpi, domdecOptions.numPmeRanks);
 820         }
 821         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 822
 823         /* Determine how many thread-MPI ranks to start.
 824          *
 825          * TODO Over-writing the user-supplied value here does
 826          * prevent any possible subsequent checks from working
 827          * correctly. */
 828         hw_opt.nthreads_tmpi =
 829                 get_nthreads_mpi(hwinfo_, &hw_opt, numDevicesToUse, useGpuForNonbonded, useGpuForPme,
 830                                  inputrec.get(), &mtop, mdlog, membedHolder.doMembed());
 831
 832         // Now start the threads for thread MPI.
 833         spawnThreads(hw_opt.nthreads_tmpi);
 834         // The spawned threads enter mdrunner() and execution of
 835         // master and spawned threads joins at the end of this block.
 836     }
 837
 838     GMX_RELEASE_ASSERT(ms || simulationCommunicator != MPI_COMM_NULL,
 839                        "Must have valid communicator unless running a multi-simulation");
 840     CommrecHandle crHandle = init_commrec(simulationCommunicator);
 841     t_commrec*    cr       = crHandle.get();
 842     GMX_RELEASE_ASSERT(cr != nullptr, "Must have valid commrec");
 843
 844     PhysicalNodeCommunicator physicalNodeComm(libraryWorldCommunicator, gmx_physicalnode_id_hash());
 845
 846     // If we detected the topology on this system, double-check that it makes sense
 847     if (hwinfo_->hardwareTopology->isThisSystem())
 848     {
 849         hardwareTopologyDoubleCheckDetection(mdlog, *hwinfo_->hardwareTopology);
 850     }
 851
 852     if (PAR(cr))
 853     {
 854         /* now broadcast everything to the non-master nodes/threads: */
 855         if (!isSimulationMasterRank)
 856         {
 857             // Until now, only the master rank has a non-null pointer.
 858             // On non-master ranks, allocate the object that will receive data in the following call.
 859             inputrec = std::make_unique<t_inputrec>();
 860         }
 861         init_parallel(cr->mpiDefaultCommunicator, MASTER(cr), inputrec.get(), &mtop,
 862                       partialDeserializedTpr.get());
 863     }
 864     GMX_RELEASE_ASSERT(inputrec != nullptr, "All ranks should have a valid inputrec now");
 865     partialDeserializedTpr.reset(nullptr);
 866
 867     // Now the number of ranks is known to all ranks, and each knows
 868     // the inputrec read by the master rank. The ranks can now all run
 869     // the task-deciding functions and will agree on the result
 870     // without needing to communicate.
 871     const bool useDomainDecomposition = (PAR(cr) && !(EI_TPI(inputrec->eI) || inputrec->eI == eiNM));
 872
 873     // Note that these variables describe only their own node.
 874     //
 875     // Note that when bonded interactions run on a GPU they always run
 876     // alongside a nonbonded task, so do not influence task assignment
 877     // even though they affect the force calculation workload.
 878     bool useGpuForNonbonded = false;
 879     bool useGpuForPme       = false;
 880     bool useGpuForBonded    = false;
 881     bool useGpuForUpdate    = false;
 882     bool gpusWereDetected   = hwinfo_->ngpu_compatible_tot > 0;
 883     try
 884     {
 885         // It's possible that there are different numbers of GPUs on
 886         // different nodes, which is the user's responsibility to
 887         // handle. If unsuitable, we will notice that during task
 888         // assignment.
 889         auto canUseGpuForNonbonded = buildSupportsNonbondedOnGpu(nullptr);
 890         useGpuForNonbonded         = decideWhetherToUseGpusForNonbonded(
 891                 nonbondedTarget, userGpuTaskAssignment, emulateGpuNonbonded, canUseGpuForNonbonded,
 892                 gpuAccelerationOfNonbondedIsUseful(mdlog, *inputrec, !GMX_THREAD_MPI), gpusWereDetected);
 893         useGpuForPme = decideWhetherToUseGpusForPme(
 894                 useGpuForNonbonded, pmeTarget, userGpuTaskAssignment, *hwinfo_, *inputrec,
 895                 cr->sizeOfDefaultCommunicator, domdecOptions.numPmeRanks, gpusWereDetected);
 896         useGpuForBonded = decideWhetherToUseGpusForBonded(useGpuForNonbonded, useGpuForPme,
 897                                                           bondedTarget, *inputrec, mtop,
 898                                                           domdecOptions.numPmeRanks, gpusWereDetected);
 899     }
 900     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 901
 902     const PmeRunMode pmeRunMode = determinePmeRunMode(useGpuForPme, pmeFftTarget, *inputrec);
 903
 904     // Initialize development feature flags that enabled by environment variable
 905     // and report those features that are enabled.
 906     const DevelopmentFeatureFlags devFlags =
 907             manageDevelopmentFeatures(mdlog, useGpuForNonbonded, pmeRunMode);
 908
 909     const bool useModularSimulator =
 910             checkUseModularSimulator(false, inputrec.get(), doRerun, mtop, ms, replExParams,
 911                                      nullptr, doEssentialDynamics, membedHolder.doMembed());
 912
 913     // Build restraints.
 914     // TODO: hide restraint implementation details from Mdrunner.
 915     // There is nothing unique about restraints at this point as far as the
 916     // Mdrunner is concerned. The Mdrunner should just be getting a sequence of
 917     // factory functions from the SimulationContext on which to call mdModules_->add().
 918     // TODO: capture all restraints into a single RestraintModule, passed to the runner builder.
 919     for (auto&& restraint : restraintManager_->getRestraints())
 920     {
 921         auto module = RestraintMDModule::create(restraint, restraint->sites());
 922         mdModules_->add(std::move(module));
 923     }
 924
 925     // TODO: Error handling
 926     mdModules_->assignOptionsToModules(*inputrec->params, nullptr);
 927     // now that the MdModules know their options, they know which callbacks to sign up to
 928     mdModules_->subscribeToSimulationSetupNotifications();
 929     const auto& mdModulesNotifier = mdModules_->notifier().simulationSetupNotifications_;
 930
 931     if (inputrec->internalParameters != nullptr)
 932     {
 933         mdModulesNotifier.notify(*inputrec->internalParameters);
 934     }
 935
 936     if (fplog != nullptr)
 937     {
 938         pr_inputrec(fplog, 0, "Input Parameters", inputrec.get(), FALSE);
 939         fprintf(fplog, "\n");
 940     }
 941
 942     if (SIMMASTER(cr))
 943     {
 944         /* In rerun, set velocities to zero if present */
 945         if (doRerun && ((globalState->flags & (1 << estV)) != 0))
 946         {
 947             // rerun does not use velocities
 948             GMX_LOG(mdlog.info)
 949                     .asParagraph()
 950                     .appendText(
 951                             "Rerun trajectory contains velocities. Rerun does only evaluate "
 952                             "potential energy and forces. The velocities will be ignored.");
 953             for (int i = 0; i < globalState->natoms; i++)
 954             {
 955                 clear_rvec(globalState->v[i]);
 956             }
 957             globalState->flags &= ~(1 << estV);
 958         }
 959
 960         /* now make sure the state is initialized and propagated */
 961         set_state_entries(globalState.get(), inputrec.get(), useModularSimulator);
 962     }
 963
 964     /* NM and TPI parallelize over force/energy calculations, not atoms,
 965      * so we need to initialize and broadcast the global state.
 966      */
 967     if (inputrec->eI == eiNM || inputrec->eI == eiTPI)
 968     {
 969         if (!MASTER(cr))
 970         {
 971             globalState = std::make_unique<t_state>();
 972         }
 973         broadcastStateWithoutDynamics(cr->mpiDefaultCommunicator, DOMAINDECOMP(cr), PAR(cr),
 974                                       globalState.get());
 975     }
 976
 977     /* A parallel command line option consistency check that we can
 978        only do after any threads have started. */
 979     if (!PAR(cr)
 980         && (domdecOptions.numCells[XX] > 1 || domdecOptions.numCells[YY] > 1
 981             || domdecOptions.numCells[ZZ] > 1 || domdecOptions.numPmeRanks > 0))
 982     {
 983         gmx_fatal(FARGS,
 984                   "The -dd or -npme option request a parallel simulation, "
 985 #if !GMX_MPI
 986                   "but %s was compiled without threads or MPI enabled",
 987                   output_env_get_program_display_name(oenv));
 988 #elif GMX_THREAD_MPI
 989                   "but the number of MPI-threads (option -ntmpi) is not set or is 1");
 990 #else
 991                   "but %s was not started through mpirun/mpiexec or only one rank was requested "
 992                   "through mpirun/mpiexec",
 993                   output_env_get_program_display_name(oenv));
 994 #endif
 995     }
 996
 997     if (doRerun && (EI_ENERGY_MINIMIZATION(inputrec->eI) || eiNM == inputrec->eI))
 998     {
 999         gmx_fatal(FARGS,
1000                   "The .mdp file specified an energy mininization or normal mode algorithm, and "
1001                   "these are not compatible with mdrun -rerun");
1002     }
1003
1004     if (!(EEL_PME(inputrec->coulombtype) || EVDW_PME(inputrec->vdwtype)))
1005     {
1006         if (domdecOptions.numPmeRanks > 0)
1007         {
1008             gmx_fatal_collective(FARGS, cr->mpiDefaultCommunicator, MASTER(cr),
1009                                  "PME-only ranks are requested, but the system does not use PME "
1010                                  "for electrostatics or LJ");
1011         }
1012
1013         domdecOptions.numPmeRanks = 0;
1014     }
1015
1016     if (useGpuForNonbonded && domdecOptions.numPmeRanks < 0)
1017     {
1018         /* With NB GPUs we don't automatically use PME-only CPU ranks. PME ranks can
1019          * improve performance with many threads per GPU, since our OpenMP
1020          * scaling is bad, but it's difficult to automate the setup.
1021          */
1022         domdecOptions.numPmeRanks = 0;
1023     }
1024     if (useGpuForPme)
1025     {
1026         if (domdecOptions.numPmeRanks < 0)
1027         {
1028             domdecOptions.numPmeRanks = 0;
1029             // TODO possibly print a note that one can opt-in for a separate PME GPU rank?
1030         }
1031         else
1032         {
1033             GMX_RELEASE_ASSERT(domdecOptions.numPmeRanks <= 1,
1034                                "PME GPU decomposition is not supported");
1035         }
1036     }
1037
1038     /* NMR restraints must be initialized before load_checkpoint,
1039      * since with time averaging the history is added to t_state.
1040      * For proper consistency check we therefore need to extend
1041      * t_state here.
1042      * So the PME-only nodes (if present) will also initialize
1043      * the distance restraints.
1044      */
1045
1046     /* This needs to be called before read_checkpoint to extend the state */
1047     t_disresdata* disresdata;
1048     snew(disresdata, 1);
1049     init_disres(fplog, &mtop, inputrec.get(), DisResRunMode::MDRun,
1050                 MASTER(cr) ? DDRole::Master : DDRole::Agent,
1051                 PAR(cr) ? NumRanks::Multiple : NumRanks::Single, cr->mpi_comm_mysim, ms, disresdata,
1052                 globalState.get(), replExParams.exchangeInterval > 0);
1053
1054     t_oriresdata* oriresdata;
1055     snew(oriresdata, 1);
1056     init_orires(fplog, &mtop, inputrec.get(), cr, ms, globalState.get(), oriresdata);
1057
1058     auto deform = prepareBoxDeformation(
1059             globalState != nullptr ? globalState->box : box, MASTER(cr) ? DDRole::Master : DDRole::Agent,
1060             PAR(cr) ? NumRanks::Multiple : NumRanks::Single, cr->mpi_comm_mygroup, *inputrec);
1061
1062 #if GMX_FAHCORE
1063     /* We have to remember the generation's first step before reading checkpoint.
1064        This way, we can report to the F@H core both the generation's first step
1065        and the restored first step, thus making it able to distinguish between
1066        an interruption/resume and start of the n-th generation simulation.
1067        Having this information, the F@H core can correctly calculate and report
1068        the progress.
1069      */
1070     int gen_first_step = 0;
1071     if (MASTER(cr))
1072     {
1073         gen_first_step = inputrec->init_step;
1074     }
1075 #endif
1076
1077     ObservablesHistory observablesHistory = {};
1078
1079     auto modularSimulatorCheckpointData = std::make_unique<ReadCheckpointDataHolder>();
1080     if (startingBehavior != StartingBehavior::NewSimulation)
1081     {
1082         /* Check if checkpoint file exists before doing continuation.
1083          * This way we can use identical input options for the first and subsequent runs...
1084          */
1085         if (mdrunOptions.numStepsCommandline > -2)
1086         {
1087             /* Temporarily set the number of steps to unlimited to avoid
1088              * triggering the nsteps check in load_checkpoint().
1089              * This hack will go away soon when the -nsteps option is removed.
1090              */
1091             inputrec->nsteps = -1;
1092         }
1093
1094         // Finish applying initial simulation state information from external sources on all ranks.
1095         // Reconcile checkpoint file data with Mdrunner state established up to this point.
1096         applyLocalState(*inputHolder_.get(), logFileHandle, cr, domdecOptions.numCells,
1097                         inputrec.get(), globalState.get(), &observablesHistory,
1098                         mdrunOptions.reproducible, mdModules_->notifier(),
1099                         modularSimulatorCheckpointData.get(), useModularSimulator);
1100         // TODO: (#3652) Synchronize filesystem state, SimulationInput contents, and program
1101         //  invariants
1102         //  on all code paths.
1103         // Write checkpoint or provide hook to update SimulationInput.
1104         // If there was a checkpoint file, SimulationInput contains more information
1105         // than if there wasn't. At this point, we have synchronized the in-memory
1106         // state with the filesystem state only for restarted simulations. We should
1107         // be calling applyLocalState unconditionally and expect that the completeness
1108         // of SimulationInput is not dependent on its creation method.
1109
1110         if (startingBehavior == StartingBehavior::RestartWithAppending && logFileHandle)
1111         {
1112             // Now we can start normal logging to the truncated log file.
1113             fplog = gmx_fio_getfp(logFileHandle);
1114             prepareLogAppending(fplog);
1115             logOwner = buildLogger(fplog, MASTER(cr));
1116             mdlog    = logOwner.logger();
1117         }
1118     }
1119
1120 #if GMX_FAHCORE
1121     if (MASTER(cr))
1122     {
1123         fcRegisterSteps(inputrec->nsteps + inputrec->init_step, gen_first_step);
1124     }
1125 #endif
1126
1127     if (mdrunOptions.numStepsCommandline > -2)
1128     {
1129         GMX_LOG(mdlog.info)
1130                 .asParagraph()
1131                 .appendText(
1132                         "The -nsteps functionality is deprecated, and may be removed in a future "
1133                         "version. "
1134                         "Consider using gmx convert-tpr -nsteps or changing the appropriate .mdp "
1135                         "file field.");
1136     }
1137     /* override nsteps with value set on the commandline */
1138     override_nsteps_cmdline(mdlog, mdrunOptions.numStepsCommandline, inputrec.get());
1139
1140     if (isSimulationMasterRank)
1141     {
1142         copy_mat(globalState->box, box);
1143     }
1144
1145     if (PAR(cr))
1146     {
1147         gmx_bcast(sizeof(box), box, cr->mpiDefaultCommunicator);
1148     }
1149
1150     if (inputrec->cutoff_scheme != ecutsVERLET)
1151     {
1152         gmx_fatal(FARGS,
1153                   "This group-scheme .tpr file can no longer be run by mdrun. Please update to the "
1154                   "Verlet scheme, or use an earlier version of GROMACS if necessary.");
1155     }
1156     /* Update rlist and nstlist. */
1157     /* Note: prepare_verlet_scheme is calling increaseNstlist(...), which (while attempting to
1158      * increase rlist) tries to check if the newly chosen value fits with the DD scheme. As this is
1159      * run before any DD scheme is set up, this check is never executed. See #3334 for more details.
1160      */
1161     prepare_verlet_scheme(fplog, cr, inputrec.get(), nstlist_cmdline, &mtop, box,
1162                           useGpuForNonbonded || (emulateGpuNonbonded == EmulateGpuNonbonded::Yes),
1163                           *hwinfo_->cpuInfo);
1164
1165     // This builder is necessary while we have multi-part construction
1166     // of DD. Before DD is constructed, we use the existence of
1167     // the builder object to indicate that further construction of DD
1168     // is needed.
1169     std::unique_ptr<DomainDecompositionBuilder> ddBuilder;
1170     if (useDomainDecomposition)
1171     {
1172         ddBuilder = std::make_unique<DomainDecompositionBuilder>(
1173                 mdlog, cr, domdecOptions, mdrunOptions, mtop, *inputrec, box,
1174                 positionsFromStatePointer(globalState.get()));
1175     }
1176     else
1177     {
1178         /* PME, if used, is done on all nodes with 1D decomposition */
1179         cr->nnodes     = cr->sizeOfDefaultCommunicator;
1180         cr->sim_nodeid = cr->rankInDefaultCommunicator;
1181         cr->nodeid     = cr->rankInDefaultCommunicator;
1182         cr->npmenodes  = 0;
1183         cr->duty       = (DUTY_PP | DUTY_PME);
1184
1185         if (inputrec->pbcType == PbcType::Screw)
1186         {
1187             gmx_fatal(FARGS, "pbc=screw is only implemented with domain decomposition");
1188         }
1189     }
1190
1191     // Produce the task assignment for this rank - done after DD is constructed
1192     GpuTaskAssignments gpuTaskAssignments = GpuTaskAssignmentsBuilder::build(
1193             gpuIdsToUse, userGpuTaskAssignment, *hwinfo_, simulationCommunicator, physicalNodeComm,
1194             nonbondedTarget, pmeTarget, bondedTarget, updateTarget, useGpuForNonbonded,
1195             useGpuForPme, thisRankHasDuty(cr, DUTY_PP),
1196             // TODO cr->duty & DUTY_PME should imply that a PME
1197             // algorithm is active, but currently does not.
1198             EEL_PME(inputrec->coulombtype) && thisRankHasDuty(cr, DUTY_PME));
1199
1200     // Get the device handles for the modules, nullptr when no task is assigned.
1201     int                deviceId   = -1;
1202     DeviceInformation* deviceInfo = gpuTaskAssignments.initDevice(&deviceId);
1203
1204     // timing enabling - TODO put this in gpu_utils (even though generally this is just option handling?)
1205     bool useTiming = true;
1206
1207     if (GMX_GPU_CUDA)
1208     {
1209         /* WARNING: CUDA timings are incorrect with multiple streams.
1210          *          This is the main reason why they are disabled by default.
1211          */
1212         // TODO: Consider turning on by default when we can detect nr of streams.
1213         useTiming = (getenv("GMX_ENABLE_GPU_TIMING") != nullptr);
1214     }
1215     else if (GMX_GPU_OPENCL)
1216     {
1217         useTiming = (getenv("GMX_DISABLE_GPU_TIMING") == nullptr);
1218     }
1219
1220     // TODO Currently this is always built, yet DD partition code
1221     // checks if it is built before using it. Probably it should
1222     // become an MDModule that is made only when another module
1223     // requires it (e.g. pull, CompEl, density fitting), so that we
1224     // don't update the local atom sets unilaterally every step.
1225     LocalAtomSetManager atomSets;
1226     if (ddBuilder)
1227     {
1228         // TODO Pass the GPU streams to ddBuilder to use in buffer
1229         // transfers (e.g. halo exchange)
1230         cr->dd = ddBuilder->build(&atomSets);
1231         // The builder's job is done, so destruct it
1232         ddBuilder.reset(nullptr);
1233         // Note that local state still does not exist yet.
1234     }
1235
1236     // The GPU update is decided here because we need to know whether the constraints or
1237     // SETTLEs can span accross the domain borders (i.e. whether or not update groups are
1238     // defined). This is only known after DD is initialized, hence decision on using GPU
1239     // update is done so late.
1240     try
1241     {
1242         const bool useUpdateGroups = cr->dd ? ddUsesUpdateGroups(*cr->dd) : false;
1243
1244         useGpuForUpdate = decideWhetherToUseGpuForUpdate(
1245                 useDomainDecomposition, useUpdateGroups, pmeRunMode, domdecOptions.numPmeRanks > 0,
1246                 useGpuForNonbonded, updateTarget, gpusWereDetected, *inputrec, mtop,
1247                 doEssentialDynamics, gmx_mtop_ftype_count(mtop, F_ORIRES) > 0,
1248                 replExParams.exchangeInterval > 0, doRerun, devFlags, mdlog);
1249     }
1250     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1251
1252     const bool printHostName = (cr->nnodes > 1);
1253     gpuTaskAssignments.reportGpuUsage(mdlog, printHostName, useGpuForBonded, pmeRunMode, useGpuForUpdate);
1254
1255     const bool disableNonbondedCalculation = (getenv("GMX_NO_NONBONDED") != nullptr);
1256     if (disableNonbondedCalculation)
1257     {
1258         /* turn off non-bonded calculations */
1259         GMX_LOG(mdlog.warning)
1260                 .asParagraph()
1261                 .appendText(
1262                         "Found environment variable GMX_NO_NONBONDED.\n"
1263                         "Disabling nonbonded calculations.");
1264     }
1265
1266     MdrunScheduleWorkload runScheduleWork;
1267
1268     bool useGpuDirectHalo = decideWhetherToUseGpuForHalo(
1269             devFlags, havePPDomainDecomposition(cr), useGpuForNonbonded, useModularSimulator,
1270             doRerun, EI_ENERGY_MINIMIZATION(inputrec->eI));
1271
1272     // Also populates the simulation constant workload description.
1273     runScheduleWork.simulationWork = createSimulationWorkload(
1274             *inputrec, disableNonbondedCalculation, devFlags, useGpuForNonbonded, pmeRunMode,
1275             useGpuForBonded, useGpuForUpdate, useGpuDirectHalo);
1276
1277     std::unique_ptr<DeviceStreamManager> deviceStreamManager = nullptr;
1278
1279     if (deviceInfo != nullptr)
1280     {
1281         if (DOMAINDECOMP(cr) && thisRankHasDuty(cr, DUTY_PP))
1282         {
1283             dd_setup_dlb_resource_sharing(cr, deviceId);
1284         }
1285         deviceStreamManager = std::make_unique<DeviceStreamManager>(
1286                 *deviceInfo, havePPDomainDecomposition(cr), runScheduleWork.simulationWork, useTiming);
1287     }
1288
1289     // If the user chose a task assignment, give them some hints
1290     // where appropriate.
1291     if (!userGpuTaskAssignment.empty())
1292     {
1293         gpuTaskAssignments.logPerformanceHints(mdlog, numDevicesToUse);
1294     }
1295
1296     if (PAR(cr))
1297     {
1298         /* After possible communicator splitting in make_dd_communicators.
1299          * we can set up the intra/inter node communication.
1300          */
1301         gmx_setup_nodecomm(fplog, cr);
1302     }
1303
1304 #if GMX_MPI
1305     if (isMultiSim(ms))
1306     {
1307         GMX_LOG(mdlog.warning)
1308                 .asParagraph()
1309                 .appendTextFormatted(
1310                         "This is simulation %d out of %d running as a composite GROMACS\n"
1311                         "multi-simulation job. Setup for this simulation:\n",
1312                         ms->simulationIndex_, ms->numSimulations_);
1313     }
1314     GMX_LOG(mdlog.warning)
1315             .appendTextFormatted("Using %d MPI %s\n", cr->nnodes,
1316 #    if GMX_THREAD_MPI
1317                                  cr->nnodes == 1 ? "thread" : "threads"
1318 #    else
1319                                  cr->nnodes == 1 ? "process" : "processes"
1320 #    endif
1321             );
1322     fflush(stderr);
1323 #endif
1324
1325     // If mdrun -pin auto honors any affinity setting that already
1326     // exists. If so, it is nice to provide feedback about whether
1327     // that existing affinity setting was from OpenMP or something
1328     // else, so we run this code both before and after we initialize
1329     // the OpenMP support.
1330     gmx_check_thread_affinity_set(mdlog, &hw_opt, hwinfo_->nthreads_hw_avail, FALSE);
1331     /* Check and update the number of OpenMP threads requested */
1332     checkAndUpdateRequestedNumOpenmpThreads(&hw_opt, *hwinfo_, cr, ms, physicalNodeComm.size_,
1333                                             pmeRunMode, mtop, *inputrec);
1334
1335     gmx_omp_nthreads_init(mdlog, cr, hwinfo_->nthreads_hw_avail, physicalNodeComm.size_,
1336                           hw_opt.nthreads_omp, hw_opt.nthreads_omp_pme, !thisRankHasDuty(cr, DUTY_PP));
1337
1338     // Enable FP exception detection, but not in
1339     // Release mode and not for compilers with known buggy FP
1340     // exception support (clang with any optimization) or suspected
1341     // buggy FP exception support (gcc 7.* with optimization).
1342 #if !defined NDEBUG                                                                         \
1343         && !((defined __clang__ || (defined(__GNUC__) && !defined(__ICC) && __GNUC__ == 7)) \
1344              && defined __OPTIMIZE__)
1345     const bool bEnableFPE = true;
1346 #else
1347     const bool bEnableFPE = false;
1348 #endif
1349     // FIXME - reconcile with gmx_feenableexcept() call from CommandLineModuleManager::run()
1350     if (bEnableFPE)
1351     {
1352         gmx_feenableexcept();
1353     }
1354
1355     /* Now that we know the setup is consistent, check for efficiency */
1356     check_resource_division_efficiency(hwinfo_, gpuTaskAssignments.thisRankHasAnyGpuTask(),
1357                                        mdrunOptions.ntompOptionIsSet, cr, mdlog);
1358
1359     /* getting number of PP/PME threads on this MPI / tMPI rank.
1360        PME: env variable should be read only on one node to make sure it is
1361        identical everywhere;
1362      */
1363     const int numThreadsOnThisRank = thisRankHasDuty(cr, DUTY_PP) ? gmx_omp_nthreads_get(emntNonbonded)
1364                                                                   : gmx_omp_nthreads_get(emntPME);
1365     checkHardwareOversubscription(numThreadsOnThisRank, cr->nodeid, *hwinfo_->hardwareTopology,
1366                                   physicalNodeComm, mdlog);
1367
1368     // Enable Peer access between GPUs where available
1369     // Only for DD, only master PP rank needs to perform setup, and only if thread MPI plus
1370     // any of the GPU communication features are active.
1371     if (DOMAINDECOMP(cr) && MASTER(cr) && thisRankHasDuty(cr, DUTY_PP) && GMX_THREAD_MPI
1372         && (runScheduleWork.simulationWork.useGpuHaloExchange
1373             || runScheduleWork.simulationWork.useGpuPmePpCommunication))
1374     {
1375         setupGpuDevicePeerAccess(gpuIdsToUse, mdlog);
1376     }
1377
1378     if (hw_opt.threadAffinity != ThreadAffinity::Off)
1379     {
1380         /* Before setting affinity, check whether the affinity has changed
1381          * - which indicates that probably the OpenMP library has changed it
1382          * since we first checked).
1383          */
1384         gmx_check_thread_affinity_set(mdlog, &hw_opt, hwinfo_->nthreads_hw_avail, TRUE);
1385
1386         int numThreadsOnThisNode, intraNodeThreadOffset;
1387         analyzeThreadsOnThisNode(physicalNodeComm, numThreadsOnThisRank, &numThreadsOnThisNode,
1388                                  &intraNodeThreadOffset);
1389
1390         /* Set the CPU affinity */
1391         gmx_set_thread_affinity(mdlog, cr, &hw_opt, *hwinfo_->hardwareTopology, numThreadsOnThisRank,
1392                                 numThreadsOnThisNode, intraNodeThreadOffset, nullptr);
1393     }
1394
1395     if (mdrunOptions.timingOptions.resetStep > -1)
1396     {
1397         GMX_LOG(mdlog.info)
1398                 .asParagraph()
1399                 .appendText(
1400                         "The -resetstep functionality is deprecated, and may be removed in a "
1401                         "future version.");
1402     }
1403     wcycle = wallcycle_init(fplog, mdrunOptions.timingOptions.resetStep, cr);
1404
1405     if (PAR(cr))
1406     {
1407         /* Master synchronizes its value of reset_counters with all nodes
1408          * including PME only nodes */
1409         int64_t reset_counters = wcycle_get_reset_counters(wcycle);
1410         gmx_bcast(sizeof(reset_counters), &reset_counters, cr->mpi_comm_mysim);
1411         wcycle_set_reset_counters(wcycle, reset_counters);
1412     }
1413
1414     // Membrane embedding must be initialized before we call init_forcerec()
1415     membedHolder.initializeMembed(fplog, filenames.size(), filenames.data(), &mtop, inputrec.get(),
1416                                   globalState.get(), cr, &mdrunOptions.checkpointOptions.period);
1417
1418     const bool               thisRankHasPmeGpuTask = gpuTaskAssignments.thisRankHasPmeGpuTask();
1419     std::unique_ptr<MDAtoms> mdAtoms;
1420     std::unique_ptr<VirtualSitesHandler> vsite;
1421     std::unique_ptr<GpuBonded>           gpuBonded;
1422
1423     t_nrnb nrnb;
1424     if (thisRankHasDuty(cr, DUTY_PP))
1425     {
1426         mdModulesNotifier.notify(*cr);
1427         mdModulesNotifier.notify(&atomSets);
1428         mdModulesNotifier.notify(inputrec->pbcType);
1429         mdModulesNotifier.notify(SimulationTimeStep{ inputrec->delta_t });
1430         /* Initiate forcerecord */
1431         fr                 = new t_forcerec;
1432         fr->forceProviders = mdModules_->initForceProviders();
1433         init_forcerec(fplog, mdlog, fr, inputrec.get(), &mtop, cr, box,
1434                       opt2fn("-table", filenames.size(), filenames.data()),
1435                       opt2fn("-tablep", filenames.size(), filenames.data()),
1436                       opt2fns("-tableb", filenames.size(), filenames.data()), pforce);
1437         // Dirty hack, for fixing disres and orires should be made mdmodules
1438         fr->fcdata->disres = disresdata;
1439         fr->fcdata->orires = oriresdata;
1440
1441         // Save a handle to device stream manager to use elsewhere in the code
1442         // TODO: Forcerec is not a correct place to store it.
1443         fr->deviceStreamManager = deviceStreamManager.get();
1444
1445         if (runScheduleWork.simulationWork.useGpuPmePpCommunication && !thisRankHasDuty(cr, DUTY_PME))
1446         {
1447             GMX_RELEASE_ASSERT(
1448                     deviceStreamManager != nullptr,
1449                     "GPU device stream manager should be valid in order to use PME-PP direct "
1450                     "communications.");
1451             GMX_RELEASE_ASSERT(
1452                     deviceStreamManager->streamIsValid(DeviceStreamType::PmePpTransfer),
1453                     "GPU PP-PME stream should be valid in order to use GPU PME-PP direct "
1454                     "communications.");
1455             fr->pmePpCommGpu = std::make_unique<gmx::PmePpCommGpu>(
1456                     cr->mpi_comm_mysim, cr->dd->pme_nodeid, deviceStreamManager->context(),
1457                     deviceStreamManager->stream(DeviceStreamType::PmePpTransfer));
1458         }
1459
1460         fr->nbv = Nbnxm::init_nb_verlet(mdlog, inputrec.get(), fr, cr, *hwinfo_,
1461                                         runScheduleWork.simulationWork.useGpuNonbonded,
1462                                         deviceStreamManager.get(), &mtop, box, wcycle);
1463         // TODO: Move the logic below to a GPU bonded builder
1464         if (runScheduleWork.simulationWork.useGpuBonded)
1465         {
1466             GMX_RELEASE_ASSERT(deviceStreamManager != nullptr,
1467                                "GPU device stream manager should be valid in order to use GPU "
1468                                "version of bonded forces.");
1469             gpuBonded = std::make_unique<GpuBonded>(
1470                     mtop.ffparams, fr->ic->epsfac * fr->fudgeQQ, deviceStreamManager->context(),
1471                     deviceStreamManager->bondedStream(havePPDomainDecomposition(cr)), wcycle);
1472             fr->gpuBonded = gpuBonded.get();
1473         }
1474
1475         /* Initialize the mdAtoms structure.
1476          * mdAtoms is not filled with atom data,
1477          * as this can not be done now with domain decomposition.
1478          */
1479         mdAtoms = makeMDAtoms(fplog, mtop, *inputrec, thisRankHasPmeGpuTask);
1480         if (globalState && thisRankHasPmeGpuTask)
1481         {
1482             // The pinning of coordinates in the global state object works, because we only use
1483             // PME on GPU without DD or on a separate PME rank, and because the local state pointer
1484             // points to the global state object without DD.
1485             // FIXME: MD and EM separately set up the local state - this should happen in the same
1486             // function, which should also perform the pinning.
1487             changePinningPolicy(&globalState->x, pme_get_pinning_policy());
1488         }
1489
1490         /* Initialize the virtual site communication */
1491         vsite = makeVirtualSitesHandler(mtop, cr, fr->pbcType);
1492
1493         calc_shifts(box, fr->shift_vec);
1494
1495         /* With periodic molecules the charge groups should be whole at start up
1496          * and the virtual sites should not be far from their proper positions.
1497          */
1498         if (!inputrec->bContinuation && MASTER(cr)
1499             && !(inputrec->pbcType != PbcType::No && inputrec->bPeriodicMols))
1500         {
1501             /* Make molecules whole at start of run */
1502             if (fr->pbcType != PbcType::No)
1503             {
1504                 do_pbc_first_mtop(fplog, inputrec->pbcType, box, &mtop, globalState->x.rvec_array());
1505             }
1506             if (vsite)
1507             {
1508                 /* Correct initial vsite positions are required
1509                  * for the initial distribution in the domain decomposition
1510                  * and for the initial shell prediction.
1511                  */
1512                 constructVirtualSitesGlobal(mtop, globalState->x);
1513             }
1514         }
1515
1516         if (EEL_PME(fr->ic->eeltype) || EVDW_PME(fr->ic->vdwtype))
1517         {
1518             ewaldcoeff_q  = fr->ic->ewaldcoeff_q;
1519             ewaldcoeff_lj = fr->ic->ewaldcoeff_lj;
1520         }
1521     }
1522     else
1523     {
1524         /* This is a PME only node */
1525
1526         GMX_ASSERT(globalState == nullptr,
1527                    "We don't need the state on a PME only rank and expect it to be unitialized");
1528
1529         ewaldcoeff_q  = calc_ewaldcoeff_q(inputrec->rcoulomb, inputrec->ewald_rtol);
1530         ewaldcoeff_lj = calc_ewaldcoeff_lj(inputrec->rvdw, inputrec->ewald_rtol_lj);
1531     }
1532
1533     gmx_pme_t* sepPmeData = nullptr;
1534     // This reference hides the fact that PME data is owned by runner on PME-only ranks and by forcerec on other ranks
1535     GMX_ASSERT(thisRankHasDuty(cr, DUTY_PP) == (fr != nullptr),
1536                "Double-checking that only PME-only ranks have no forcerec");
1537     gmx_pme_t*& pmedata = fr ? fr->pmedata : sepPmeData;
1538
1539     // TODO should live in ewald module once its testing is improved
1540     //
1541     // Later, this program could contain kernels that might be later
1542     // re-used as auto-tuning progresses, or subsequent simulations
1543     // are invoked.
1544     PmeGpuProgramStorage pmeGpuProgram;
1545     if (thisRankHasPmeGpuTask)
1546     {
1547         GMX_RELEASE_ASSERT(
1548                 (deviceStreamManager != nullptr),
1549                 "GPU device stream manager should be initialized in order to use GPU for PME.");
1550         GMX_RELEASE_ASSERT((deviceInfo != nullptr),
1551                            "GPU device should be initialized in order to use GPU for PME.");
1552         pmeGpuProgram = buildPmeGpuProgram(deviceStreamManager->context());
1553     }
1554
1555     /* Initiate PME if necessary,
1556      * either on all nodes or on dedicated PME nodes only. */
1557     if (EEL_PME(inputrec->coulombtype) || EVDW_PME(inputrec->vdwtype))
1558     {
1559         if (mdAtoms && mdAtoms->mdatoms())
1560         {
1561             nChargePerturbed = mdAtoms->mdatoms()->nChargePerturbed;
1562             if (EVDW_PME(inputrec->vdwtype))
1563             {
1564                 nTypePerturbed = mdAtoms->mdatoms()->nTypePerturbed;
1565             }
1566         }
1567         if (cr->npmenodes > 0)
1568         {
1569             /* The PME only nodes need to know nChargePerturbed(FEP on Q) and nTypePerturbed(FEP on LJ)*/
1570             gmx_bcast(sizeof(nChargePerturbed), &nChargePerturbed, cr->mpi_comm_mysim);
1571             gmx_bcast(sizeof(nTypePerturbed), &nTypePerturbed, cr->mpi_comm_mysim);
1572         }
1573
1574         if (thisRankHasDuty(cr, DUTY_PME))
1575         {
1576             try
1577             {
1578                 // TODO: This should be in the builder.
1579                 GMX_RELEASE_ASSERT(!runScheduleWork.simulationWork.useGpuPme
1580                                            || (deviceStreamManager != nullptr),
1581                                    "Device stream manager should be valid in order to use GPU "
1582                                    "version of PME.");
1583                 GMX_RELEASE_ASSERT(
1584                         !runScheduleWork.simulationWork.useGpuPme
1585                                 || deviceStreamManager->streamIsValid(DeviceStreamType::Pme),
1586                         "GPU PME stream should be valid in order to use GPU version of PME.");
1587
1588                 const DeviceContext* deviceContext = runScheduleWork.simulationWork.useGpuPme
1589                                                              ? &deviceStreamManager->context()
1590                                                              : nullptr;
1591                 const DeviceStream* pmeStream =
1592                         runScheduleWork.simulationWork.useGpuPme
1593                                 ? &deviceStreamManager->stream(DeviceStreamType::Pme)
1594                                 : nullptr;
1595
1596                 pmedata = gmx_pme_init(cr, getNumPmeDomains(cr->dd), inputrec.get(),
1597                                        nChargePerturbed != 0, nTypePerturbed != 0,
1598                                        mdrunOptions.reproducible, ewaldcoeff_q, ewaldcoeff_lj,
1599                                        gmx_omp_nthreads_get(emntPME), pmeRunMode, nullptr,
1600                                        deviceContext, pmeStream, pmeGpuProgram.get(), mdlog);
1601             }
1602             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1603         }
1604     }
1605
1606
1607     if (EI_DYNAMICS(inputrec->eI))
1608     {
1609         /* Turn on signal handling on all nodes */
1610         /*
1611          * (A user signal from the PME nodes (if any)
1612          * is communicated to the PP nodes.
1613          */
1614         signal_handler_install();
1615     }
1616
1617     pull_t* pull_work = nullptr;
1618     if (thisRankHasDuty(cr, DUTY_PP))
1619     {
1620         /* Assumes uniform use of the number of OpenMP threads */
1621         walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntDefault));
1622
1623         if (inputrec->bPull)
1624         {
1625             /* Initialize pull code */
1626             pull_work = init_pull(fplog, inputrec->pull.get(), inputrec.get(), &mtop, cr, &atomSets,
1627                                   inputrec->fepvals->init_lambda);
1628             if (inputrec->pull->bXOutAverage || inputrec->pull->bFOutAverage)
1629             {
1630                 initPullHistory(pull_work, &observablesHistory);
1631             }
1632             if (EI_DYNAMICS(inputrec->eI) && MASTER(cr))
1633             {
1634                 init_pull_output_files(pull_work, filenames.size(), filenames.data(), oenv, startingBehavior);
1635             }
1636         }
1637
1638         std::unique_ptr<EnforcedRotation> enforcedRotation;
1639         if (inputrec->bRot)
1640         {
1641             /* Initialize enforced rotation code */
1642             enforcedRotation = init_rot(fplog, inputrec.get(), filenames.size(), filenames.data(),
1643                                         cr, &atomSets, globalState.get(), &mtop, oenv, mdrunOptions,
1644                                         startingBehavior);
1645         }
1646
1647         t_swap* swap = nullptr;
1648         if (inputrec->eSwapCoords != eswapNO)
1649         {
1650             /* Initialize ion swapping code */
1651             swap = init_swapcoords(fplog, inputrec.get(),
1652                                    opt2fn_master("-swap", filenames.size(), filenames.data(), cr),
1653                                    &mtop, globalState.get(), &observablesHistory, cr, &atomSets,
1654                                    oenv, mdrunOptions, startingBehavior);
1655         }
1656
1657         /* Let makeConstraints know whether we have essential dynamics constraints. */
1658         auto constr = makeConstraints(mtop, *inputrec, pull_work, doEssentialDynamics, fplog, cr,
1659                                       ms, &nrnb, wcycle, fr->bMolPBC);
1660
1661         /* Energy terms and groups */
1662         gmx_enerdata_t enerd(mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
1663                              inputrec->fepvals->n_lambda);
1664
1665         /* Kinetic energy data */
1666         gmx_ekindata_t ekind;
1667         init_ekindata(fplog, &mtop, &(inputrec->opts), &ekind, inputrec->cos_accel);
1668
1669         /* Set up interactive MD (IMD) */
1670         auto imdSession =
1671                 makeImdSession(inputrec.get(), cr, wcycle, &enerd, ms, &mtop, mdlog,
1672                                MASTER(cr) ? globalState->x.rvec_array() : nullptr, filenames.size(),
1673                                filenames.data(), oenv, mdrunOptions.imdOptions, startingBehavior);
1674
1675         if (DOMAINDECOMP(cr))
1676         {
1677             GMX_RELEASE_ASSERT(fr, "fr was NULL while cr->duty was DUTY_PP");
1678             /* This call is not included in init_domain_decomposition mainly
1679              * because fr->cginfo_mb is set later.
1680              */
1681             dd_init_bondeds(fplog, cr->dd, mtop, vsite.get(), inputrec.get(),
1682                             domdecOptions.checkBondedInteractions, fr->cginfo_mb);
1683         }
1684
1685         if (runScheduleWork.simulationWork.useGpuBufferOps)
1686         {
1687             fr->gpuForceReduction[gmx::AtomLocality::Local] = std::make_unique<gmx::GpuForceReduction>(
1688                     deviceStreamManager->context(),
1689                     deviceStreamManager->stream(gmx::DeviceStreamType::NonBondedLocal), wcycle);
1690             fr->gpuForceReduction[gmx::AtomLocality::NonLocal] = std::make_unique<gmx::GpuForceReduction>(
1691                     deviceStreamManager->context(),
1692                     deviceStreamManager->stream(gmx::DeviceStreamType::NonBondedNonLocal), wcycle);
1693         }
1694
1695         std::unique_ptr<gmx::StatePropagatorDataGpu> stateGpu;
1696         if (gpusWereDetected
1697             && ((runScheduleWork.simulationWork.useGpuPme && thisRankHasDuty(cr, DUTY_PME))
1698                 || runScheduleWork.simulationWork.useGpuBufferOps))
1699         {
1700             GpuApiCallBehavior transferKind = (inputrec->eI == eiMD && !doRerun && !useModularSimulator)
1701                                                       ? GpuApiCallBehavior::Async
1702                                                       : GpuApiCallBehavior::Sync;
1703             GMX_RELEASE_ASSERT(deviceStreamManager != nullptr,
1704                                "GPU device stream manager should be initialized to use GPU.");
1705             stateGpu = std::make_unique<gmx::StatePropagatorDataGpu>(
1706                     *deviceStreamManager, transferKind, pme_gpu_get_block_size(fr->pmedata), wcycle);
1707             fr->stateGpu = stateGpu.get();
1708         }
1709
1710         GMX_ASSERT(stopHandlerBuilder_, "Runner must provide StopHandlerBuilder to simulator.");
1711         SimulatorBuilder simulatorBuilder;
1712
1713         simulatorBuilder.add(SimulatorStateData(globalState.get(), &observablesHistory, &enerd, &ekind));
1714         simulatorBuilder.add(std::move(membedHolder));
1715         simulatorBuilder.add(std::move(stopHandlerBuilder_));
1716         simulatorBuilder.add(SimulatorConfig(mdrunOptions, startingBehavior, &runScheduleWork));
1717
1718
1719         simulatorBuilder.add(SimulatorEnv(fplog, cr, ms, mdlog, oenv));
1720         simulatorBuilder.add(Profiling(&nrnb, walltime_accounting, wcycle));
1721         simulatorBuilder.add(ConstraintsParam(
1722                 constr.get(), enforcedRotation ? enforcedRotation->getLegacyEnfrot() : nullptr,
1723                 vsite.get()));
1724         // TODO: Separate `fr` to a separate add, and make the `build` handle the coupling sensibly.
1725         simulatorBuilder.add(LegacyInput(static_cast<int>(filenames.size()), filenames.data(),
1726                                          inputrec.get(), fr));
1727         simulatorBuilder.add(ReplicaExchangeParameters(replExParams));
1728         simulatorBuilder.add(InteractiveMD(imdSession.get()));
1729         simulatorBuilder.add(SimulatorModules(mdModules_->outputProvider(), mdModules_->notifier()));
1730         simulatorBuilder.add(CenterOfMassPulling(pull_work));
1731         // Todo move to an MDModule
1732         simulatorBuilder.add(IonSwapping(swap));
1733         simulatorBuilder.add(TopologyData(&mtop, mdAtoms.get()));
1734         simulatorBuilder.add(BoxDeformationHandle(deform.get()));
1735         simulatorBuilder.add(std::move(modularSimulatorCheckpointData));
1736
1737         // build and run simulator object based on user-input
1738         auto simulator = simulatorBuilder.build(useModularSimulator);
1739         simulator->run();
1740
1741         if (fr->pmePpCommGpu)
1742         {
1743             // destroy object since it is no longer required. (This needs to be done while the GPU context still exists.)
1744             fr->pmePpCommGpu.reset();
1745         }
1746
1747         if (inputrec->bPull)
1748         {
1749             finish_pull(pull_work);
1750         }
1751         finish_swapcoords(swap);
1752     }
1753     else
1754     {
1755         GMX_RELEASE_ASSERT(pmedata, "pmedata was NULL while cr->duty was not DUTY_PP");
1756         /* do PME only */
1757         walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntPME));
1758         gmx_pmeonly(pmedata, cr, &nrnb, wcycle, walltime_accounting, inputrec.get(), pmeRunMode,
1759                     deviceStreamManager.get());
1760     }
1761
1762     wallcycle_stop(wcycle, ewcRUN);
1763
1764     /* Finish up, write some stuff
1765      * if rerunMD, don't write last frame again
1766      */
1767     finish_run(fplog, mdlog, cr, inputrec.get(), &nrnb, wcycle, walltime_accounting,
1768                fr ? fr->nbv.get() : nullptr, pmedata, EI_DYNAMICS(inputrec->eI) && !isMultiSim(ms));
1769
1770     // clean up cycle counter
1771     wallcycle_destroy(wcycle);
1772
1773     deviceStreamManager.reset(nullptr);
1774     // Free PME data
1775     if (pmedata)
1776     {
1777         gmx_pme_destroy(pmedata);
1778         pmedata = nullptr;
1779     }
1780
1781     // FIXME: this is only here to manually unpin mdAtoms->chargeA_ and state->x,
1782     // before we destroy the GPU context(s)
1783     // Pinned buffers are associated with contexts in CUDA.
1784     // As soon as we destroy GPU contexts after mdrunner() exits, these lines should go.
1785     mdAtoms.reset(nullptr);
1786     globalState.reset(nullptr);
1787     mdModules_.reset(nullptr); // destruct force providers here as they might also use the GPU
1788     gpuBonded.reset(nullptr);
1789     /* Free pinned buffers in *fr */
1790     delete fr;
1791     fr = nullptr;
1792     // TODO convert to C++ so we can get rid of these frees
1793     sfree(disresdata);
1794     sfree(oriresdata);
1795
1796     if (!hwinfo_->deviceInfoList.empty())
1797     {
1798         /* stop the GPU profiler (only CUDA) */
1799         stopGpuProfiler();
1800     }
1801
1802     /* With tMPI we need to wait for all ranks to finish deallocation before
1803      * destroying the CUDA context as some tMPI ranks may be sharing
1804      * GPU and context.
1805      *
1806      * This is not a concern in OpenCL where we use one context per rank.
1807      *
1808      * Note: it is safe to not call the barrier on the ranks which do not use GPU,
1809      * but it is easier and more futureproof to call it on the whole node.
1810      *
1811      * Note that this function needs to be called even if GPUs are not used
1812      * in this run because the PME ranks have no knowledge of whether GPUs
1813      * are used or not, but all ranks need to enter the barrier below.
1814      * \todo Remove this physical node barrier after making sure
1815      * that it's not needed anymore (with a shared GPU run).
1816      */
1817     if (GMX_THREAD_MPI)
1818     {
1819         physicalNodeComm.barrier();
1820     }
1821     releaseDevice(deviceInfo);
1822
1823     /* Does what it says */
1824     print_date_and_time(fplog, cr->nodeid, "Finished mdrun", gmx_gettime());
1825     walltime_accounting_destroy(walltime_accounting);
1826
1827     // Ensure log file content is written
1828     if (logFileHandle)
1829     {
1830         gmx_fio_flush(logFileHandle);
1831     }
1832
1833     /* Reset FPEs (important for unit tests) by disabling them. Assumes no
1834      * exceptions were enabled before function was called. */
1835     if (bEnableFPE)
1836     {
1837         gmx_fedisableexcept();
1838     }
1839
1840     auto rc = static_cast<int>(gmx_get_stop_condition());
1841
1842 #if GMX_THREAD_MPI
1843     /* we need to join all threads. The sub-threads join when they
1844        exit this function, but the master thread needs to be told to
1845        wait for that. */
1846     if (MASTER(cr))
1847     {
1848         tMPI_Finalize();
1849     }
1850 #endif
1851     return rc;
1852 } // namespace gmx
1853
1854 Mdrunner::~Mdrunner()
1855 {
1856     // Clean up of the Manager.
1857     // This will end up getting called on every thread-MPI rank, which is unnecessary,
1858     // but okay as long as threads synchronize some time before adding or accessing
1859     // a new set of restraints.
1860     if (restraintManager_)
1861     {
1862         restraintManager_->clear();
1863         GMX_ASSERT(restraintManager_->countRestraints() == 0,
1864                    "restraints added during runner life time should be cleared at runner "
1865                    "destruction.");
1866     }
1867 };
1868
1869 void Mdrunner::addPotential(std::shared_ptr<gmx::IRestraintPotential> puller, const std::string& name)
1870 {
1871     GMX_ASSERT(restraintManager_, "Mdrunner must have a restraint manager.");
1872     // Not sure if this should be logged through the md logger or something else,
1873     // but it is helpful to have some sort of INFO level message sent somewhere.
1874     //    std::cout << "Registering restraint named " << name << std::endl;
1875
1876     // When multiple restraints are used, it may be wasteful to register them separately.
1877     // Maybe instead register an entire Restraint Manager as a force provider.
1878     restraintManager_->addToSpec(std::move(puller), name);
1879 }
1880
1881 Mdrunner::Mdrunner(std::unique_ptr<MDModules> mdModules) : mdModules_(std::move(mdModules)) {}
1882
1883 Mdrunner::Mdrunner(Mdrunner&&) noexcept = default;
1884
1885 //NOLINTNEXTLINE(performance-noexcept-move-constructor) working around GCC bug 58265 in CentOS 7
1886 Mdrunner& Mdrunner::operator=(Mdrunner&& /*handle*/) noexcept(BUGFREE_NOEXCEPT_STRING) = default;
1887
1888 class Mdrunner::BuilderImplementation
1889 {
1890 public:
1891     BuilderImplementation() = delete;
1892     BuilderImplementation(std::unique_ptr<MDModules> mdModules, compat::not_null<SimulationContext*> context);
1893     ~BuilderImplementation();
1894
1895     BuilderImplementation& setExtraMdrunOptions(const MdrunOptions& options,
1896                                                 real                forceWarningThreshold,
1897                                                 StartingBehavior    startingBehavior);
1898
1899     void addHardwareDetectionResult(const gmx_hw_info_t* hwinfo);
1900
1901     void addDomdec(const DomdecOptions& options);
1902
1903     void addInput(SimulationInputHandle inputHolder);
1904
1905     void addVerletList(int nstlist);
1906
1907     void addReplicaExchange(const ReplicaExchangeParameters& params);
1908
1909     void addNonBonded(const char* nbpu_opt);
1910
1911     void addPME(const char* pme_opt_, const char* pme_fft_opt_);
1912
1913     void addBondedTaskAssignment(const char* bonded_opt);
1914
1915     void addUpdateTaskAssignment(const char* update_opt);
1916
1917     void addHardwareOptions(const gmx_hw_opt_t& hardwareOptions);
1918
1919     void addFilenames(ArrayRef<const t_filenm> filenames);
1920
1921     void addOutputEnvironment(gmx_output_env_t* outputEnvironment);
1922
1923     void addLogFile(t_fileio* logFileHandle);
1924
1925     void addStopHandlerBuilder(std::unique_ptr<StopHandlerBuilder> builder);
1926
1927     Mdrunner build();
1928
1929 private:
1930     // Default parameters copied from runner.h
1931     // \todo Clarify source(s) of default parameters.
1932
1933     const char* nbpu_opt_    = nullptr;
1934     const char* pme_opt_     = nullptr;
1935     const char* pme_fft_opt_ = nullptr;
1936     const char* bonded_opt_  = nullptr;
1937     const char* update_opt_  = nullptr;
1938
1939     MdrunOptions mdrunOptions_;
1940
1941     DomdecOptions domdecOptions_;
1942
1943     ReplicaExchangeParameters replicaExchangeParameters_;
1944
1945     //! Command-line override for the duration of a neighbor list with the Verlet scheme.
1946     int nstlist_ = 0;
1947
1948     //! World communicator, used for hardware detection and task assignment
1949     MPI_Comm libraryWorldCommunicator_ = MPI_COMM_NULL;
1950
1951     //! Multisim communicator handle.
1952     gmx_multisim_t* multiSimulation_;
1953
1954     //! mdrun communicator
1955     MPI_Comm simulationCommunicator_ = MPI_COMM_NULL;
1956
1957     //! Print a warning if any force is larger than this (in kJ/mol nm).
1958     real forceWarningThreshold_ = -1;
1959
1960     //! Whether the simulation will start afresh, or restart with/without appending.
1961     StartingBehavior startingBehavior_ = StartingBehavior::NewSimulation;
1962
1963     //! The modules that comprise the functionality of mdrun.
1964     std::unique_ptr<MDModules> mdModules_;
1965
1966     //! Detected hardware.
1967     const gmx_hw_info_t* hwinfo_ = nullptr;
1968
1969     //! \brief Parallelism information.
1970     gmx_hw_opt_t hardwareOptions_;
1971
1972     //! filename options for simulation.
1973     ArrayRef<const t_filenm> filenames_;
1974
1975     /*! \brief Handle to output environment.
1976      *
1977      * \todo gmx_output_env_t needs lifetime management.
1978      */
1979     gmx_output_env_t* outputEnvironment_ = nullptr;
1980
1981     /*! \brief Non-owning handle to MD log file.
1982      *
1983      * \todo Context should own output facilities for client.
1984      * \todo Improve log file handle management.
1985      * \internal
1986      * Code managing the FILE* relies on the ability to set it to
1987      * nullptr to check whether the filehandle is valid.
1988      */
1989     t_fileio* logFileHandle_ = nullptr;
1990
1991     /*!
1992      * \brief Builder for simulation stop signal handler.
1993      */
1994     std::unique_ptr<StopHandlerBuilder> stopHandlerBuilder_ = nullptr;
1995
1996     /*!
1997      * \brief Sources for initial simulation state.
1998      *
1999      * See issue #3652 for near-term refinements to the SimulationInput interface.
2000      *
2001      * See issue #3379 for broader discussion on API aspects of simulation inputs and outputs.
2002      */
2003     SimulationInputHandle inputHolder_;
2004 };
2005
2006 Mdrunner::BuilderImplementation::BuilderImplementation(std::unique_ptr<MDModules> mdModules,
2007                                                        compat::not_null<SimulationContext*> context) :
2008     mdModules_(std::move(mdModules))
2009 {
2010     libraryWorldCommunicator_ = context->libraryWorldCommunicator_;
2011     simulationCommunicator_   = context->simulationCommunicator_;
2012     multiSimulation_          = context->multiSimulation_.get();
2013 }
2014
2015 Mdrunner::BuilderImplementation::~BuilderImplementation() = default;
2016
2017 Mdrunner::BuilderImplementation&
2018 Mdrunner::BuilderImplementation::setExtraMdrunOptions(const MdrunOptions&    options,
2019                                                       const real             forceWarningThreshold,
2020                                                       const StartingBehavior startingBehavior)
2021 {
2022     mdrunOptions_          = options;
2023     forceWarningThreshold_ = forceWarningThreshold;
2024     startingBehavior_      = startingBehavior;
2025     return *this;
2026 }
2027
2028 void Mdrunner::BuilderImplementation::addDomdec(const DomdecOptions& options)
2029 {
2030     domdecOptions_ = options;
2031 }
2032
2033 void Mdrunner::BuilderImplementation::addVerletList(int nstlist)
2034 {
2035     nstlist_ = nstlist;
2036 }
2037
2038 void Mdrunner::BuilderImplementation::addReplicaExchange(const ReplicaExchangeParameters& params)
2039 {
2040     replicaExchangeParameters_ = params;
2041 }
2042
2043 Mdrunner Mdrunner::BuilderImplementation::build()
2044 {
2045     auto newRunner = Mdrunner(std::move(mdModules_));
2046
2047     newRunner.mdrunOptions     = mdrunOptions_;
2048     newRunner.pforce           = forceWarningThreshold_;
2049     newRunner.startingBehavior = startingBehavior_;
2050     newRunner.domdecOptions    = domdecOptions_;
2051
2052     // \todo determine an invariant to check or confirm that all gmx_hw_opt_t objects are valid
2053     newRunner.hw_opt = hardwareOptions_;
2054
2055     // No invariant to check. This parameter exists to optionally override other behavior.
2056     newRunner.nstlist_cmdline = nstlist_;
2057
2058     newRunner.replExParams = replicaExchangeParameters_;
2059
2060     newRunner.filenames = filenames_;
2061
2062     newRunner.libraryWorldCommunicator = libraryWorldCommunicator_;
2063
2064     newRunner.simulationCommunicator = simulationCommunicator_;
2065
2066     // nullptr is a valid value for the multisim handle
2067     newRunner.ms = multiSimulation_;
2068
2069     if (hwinfo_)
2070     {
2071         newRunner.hwinfo_ = hwinfo_;
2072     }
2073     else
2074     {
2075         GMX_THROW(gmx::APIError(
2076                 "MdrunnerBuilder::addHardwareDetectionResult() is required before build()"));
2077     }
2078
2079     if (inputHolder_)
2080     {
2081         newRunner.inputHolder_ = std::move(inputHolder_);
2082     }
2083     else
2084     {
2085         GMX_THROW(gmx::APIError("MdrunnerBuilder::addInput() is required before build()."));
2086     }
2087
2088     // \todo Clarify ownership and lifetime management for gmx_output_env_t
2089     // \todo Update sanity checking when output environment has clearly specified invariants.
2090     // Initialization and default values for oenv are not well specified in the current version.
2091     if (outputEnvironment_)
2092     {
2093         newRunner.oenv = outputEnvironment_;
2094     }
2095     else
2096     {
2097         GMX_THROW(gmx::APIError(
2098                 "MdrunnerBuilder::addOutputEnvironment() is required before build()"));
2099     }
2100
2101     newRunner.logFileHandle = logFileHandle_;
2102
2103     if (nbpu_opt_)
2104     {
2105         newRunner.nbpu_opt = nbpu_opt_;
2106     }
2107     else
2108     {
2109         GMX_THROW(gmx::APIError("MdrunnerBuilder::addNonBonded() is required before build()"));
2110     }
2111
2112     if (pme_opt_ && pme_fft_opt_)
2113     {
2114         newRunner.pme_opt     = pme_opt_;
2115         newRunner.pme_fft_opt = pme_fft_opt_;
2116     }
2117     else
2118     {
2119         GMX_THROW(gmx::APIError("MdrunnerBuilder::addElectrostatics() is required before build()"));
2120     }
2121
2122     if (bonded_opt_)
2123     {
2124         newRunner.bonded_opt = bonded_opt_;
2125     }
2126     else
2127     {
2128         GMX_THROW(gmx::APIError(
2129                 "MdrunnerBuilder::addBondedTaskAssignment() is required before build()"));
2130     }
2131
2132     if (update_opt_)
2133     {
2134         newRunner.update_opt = update_opt_;
2135     }
2136     else
2137     {
2138         GMX_THROW(gmx::APIError(
2139                 "MdrunnerBuilder::addUpdateTaskAssignment() is required before build()  "));
2140     }
2141
2142
2143     newRunner.restraintManager_ = std::make_unique<gmx::RestraintManager>();
2144
2145     if (stopHandlerBuilder_)
2146     {
2147         newRunner.stopHandlerBuilder_ = std::move(stopHandlerBuilder_);
2148     }
2149     else
2150     {
2151         newRunner.stopHandlerBuilder_ = std::make_unique<StopHandlerBuilder>();
2152     }
2153
2154     return newRunner;
2155 }
2156
2157 void Mdrunner::BuilderImplementation::addHardwareDetectionResult(const gmx_hw_info_t* hwinfo)
2158 {
2159     hwinfo_ = hwinfo;
2160 }
2161
2162 void Mdrunner::BuilderImplementation::addNonBonded(const char* nbpu_opt)
2163 {
2164     nbpu_opt_ = nbpu_opt;
2165 }
2166
2167 void Mdrunner::BuilderImplementation::addPME(const char* pme_opt, const char* pme_fft_opt)
2168 {
2169     pme_opt_     = pme_opt;
2170     pme_fft_opt_ = pme_fft_opt;
2171 }
2172
2173 void Mdrunner::BuilderImplementation::addBondedTaskAssignment(const char* bonded_opt)
2174 {
2175     bonded_opt_ = bonded_opt;
2176 }
2177
2178 void Mdrunner::BuilderImplementation::addUpdateTaskAssignment(const char* update_opt)
2179 {
2180     update_opt_ = update_opt;
2181 }
2182
2183 void Mdrunner::BuilderImplementation::addHardwareOptions(const gmx_hw_opt_t& hardwareOptions)
2184 {
2185     hardwareOptions_ = hardwareOptions;
2186 }
2187
2188 void Mdrunner::BuilderImplementation::addFilenames(ArrayRef<const t_filenm> filenames)
2189 {
2190     filenames_ = filenames;
2191 }
2192
2193 void Mdrunner::BuilderImplementation::addOutputEnvironment(gmx_output_env_t* outputEnvironment)
2194 {
2195     outputEnvironment_ = outputEnvironment;
2196 }
2197
2198 void Mdrunner::BuilderImplementation::addLogFile(t_fileio* logFileHandle)
2199 {
2200     logFileHandle_ = logFileHandle;
2201 }
2202
2203 void Mdrunner::BuilderImplementation::addStopHandlerBuilder(std::unique_ptr<StopHandlerBuilder> builder)
2204 {
2205     stopHandlerBuilder_ = std::move(builder);
2206 }
2207
2208 void Mdrunner::BuilderImplementation::addInput(SimulationInputHandle inputHolder)
2209 {
2210     inputHolder_ = std::move(inputHolder);
2211 }
2212
2213 MdrunnerBuilder::MdrunnerBuilder(std::unique_ptr<MDModules>           mdModules,
2214                                  compat::not_null<SimulationContext*> context) :
2215     impl_{ std::make_unique<Mdrunner::BuilderImplementation>(std::move(mdModules), context) }
2216 {
2217 }
2218
2219 MdrunnerBuilder::~MdrunnerBuilder() = default;
2220
2221 MdrunnerBuilder& MdrunnerBuilder::addHardwareDetectionResult(const gmx_hw_info_t* hwinfo)
2222 {
2223     impl_->addHardwareDetectionResult(hwinfo);
2224     return *this;
2225 }
2226
2227 MdrunnerBuilder& MdrunnerBuilder::addSimulationMethod(const MdrunOptions&    options,
2228                                                       real                   forceWarningThreshold,
2229                                                       const StartingBehavior startingBehavior)
2230 {
2231     impl_->setExtraMdrunOptions(options, forceWarningThreshold, startingBehavior);
2232     return *this;
2233 }
2234
2235 MdrunnerBuilder& MdrunnerBuilder::addDomainDecomposition(const DomdecOptions& options)
2236 {
2237     impl_->addDomdec(options);
2238     return *this;
2239 }
2240
2241 MdrunnerBuilder& MdrunnerBuilder::addNeighborList(int nstlist)
2242 {
2243     impl_->addVerletList(nstlist);
2244     return *this;
2245 }
2246
2247 MdrunnerBuilder& MdrunnerBuilder::addReplicaExchange(const ReplicaExchangeParameters& params)
2248 {
2249     impl_->addReplicaExchange(params);
2250     return *this;
2251 }
2252
2253 MdrunnerBuilder& MdrunnerBuilder::addNonBonded(const char* nbpu_opt)
2254 {
2255     impl_->addNonBonded(nbpu_opt);
2256     return *this;
2257 }
2258
2259 MdrunnerBuilder& MdrunnerBuilder::addElectrostatics(const char* pme_opt, const char* pme_fft_opt)
2260 {
2261     // The builder method may become more general in the future, but in this version,
2262     // parameters for PME electrostatics are both required and the only parameters
2263     // available.
2264     if (pme_opt && pme_fft_opt)
2265     {
2266         impl_->addPME(pme_opt, pme_fft_opt);
2267     }
2268     else
2269     {
2270         GMX_THROW(
2271                 gmx::InvalidInputError("addElectrostatics() arguments must be non-null pointers."));
2272     }
2273     return *this;
2274 }
2275
2276 MdrunnerBuilder& MdrunnerBuilder::addBondedTaskAssignment(const char* bonded_opt)
2277 {
2278     impl_->addBondedTaskAssignment(bonded_opt);
2279     return *this;
2280 }
2281
2282 MdrunnerBuilder& MdrunnerBuilder::addUpdateTaskAssignment(const char* update_opt)
2283 {
2284     impl_->addUpdateTaskAssignment(update_opt);
2285     return *this;
2286 }
2287
2288 Mdrunner MdrunnerBuilder::build()
2289 {
2290     return impl_->build();
2291 }
2292
2293 MdrunnerBuilder& MdrunnerBuilder::addHardwareOptions(const gmx_hw_opt_t& hardwareOptions)
2294 {
2295     impl_->addHardwareOptions(hardwareOptions);
2296     return *this;
2297 }
2298
2299 MdrunnerBuilder& MdrunnerBuilder::addFilenames(ArrayRef<const t_filenm> filenames)
2300 {
2301     impl_->addFilenames(filenames);
2302     return *this;
2303 }
2304
2305 MdrunnerBuilder& MdrunnerBuilder::addOutputEnvironment(gmx_output_env_t* outputEnvironment)
2306 {
2307     impl_->addOutputEnvironment(outputEnvironment);
2308     return *this;
2309 }
2310
2311 MdrunnerBuilder& MdrunnerBuilder::addLogFile(t_fileio* logFileHandle)
2312 {
2313     impl_->addLogFile(logFileHandle);
2314     return *this;
2315 }
2316
2317 MdrunnerBuilder& MdrunnerBuilder::addStopHandlerBuilder(std::unique_ptr<StopHandlerBuilder> builder)
2318 {
2319     impl_->addStopHandlerBuilder(std::move(builder));
2320     return *this;
2321 }
2322
2323 MdrunnerBuilder& MdrunnerBuilder::addInput(SimulationInputHandle input)
2324 {
2325     impl_->addInput(std::move(input));
2326     return *this;
2327 }
2328
2329 MdrunnerBuilder::MdrunnerBuilder(MdrunnerBuilder&&) noexcept = default;
2330
2331 MdrunnerBuilder& MdrunnerBuilder::operator=(MdrunnerBuilder&&) noexcept = default;
2332
2333 } // namespace gmx