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