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     // Build restraints.
 897     // TODO: hide restraint implementation details from Mdrunner.
 898     // There is nothing unique about restraints at this point as far as the
 899     // Mdrunner is concerned. The Mdrunner should just be getting a sequence of
 900     // factory functions from the SimulationContext on which to call mdModules_->add().
 901     // TODO: capture all restraints into a single RestraintModule, passed to the runner builder.
 902     for (auto&& restraint : restraintManager_->getRestraints())
 903     {
 904         auto module = RestraintMDModule::create(restraint, restraint->sites());
 905         mdModules_->add(std::move(module));
 906     }
 907
 908     // TODO: Error handling
 909     mdModules_->assignOptionsToModules(*inputrec->params, nullptr);
 910     const auto& mdModulesNotifier = mdModules_->notifier().notifier_;
 911
 912     if (inputrec->internalParameters != nullptr)
 913     {
 914         mdModulesNotifier.notify(*inputrec->internalParameters);
 915     }
 916
 917     if (fplog != nullptr)
 918     {
 919         pr_inputrec(fplog, 0, "Input Parameters", inputrec, FALSE);
 920         fprintf(fplog, "\n");
 921     }
 922
 923     if (SIMMASTER(cr))
 924     {
 925         /* In rerun, set velocities to zero if present */
 926         if (doRerun && ((globalState->flags & (1 << estV)) != 0))
 927         {
 928             // rerun does not use velocities
 929             GMX_LOG(mdlog.info)
 930                     .asParagraph()
 931                     .appendText(
 932                             "Rerun trajectory contains velocities. Rerun does only evaluate "
 933                             "potential energy and forces. The velocities will be ignored.");
 934             for (int i = 0; i < globalState->natoms; i++)
 935             {
 936                 clear_rvec(globalState->v[i]);
 937             }
 938             globalState->flags &= ~(1 << estV);
 939         }
 940
 941         /* now make sure the state is initialized and propagated */
 942         set_state_entries(globalState.get(), inputrec);
 943     }
 944
 945     /* NM and TPI parallelize over force/energy calculations, not atoms,
 946      * so we need to initialize and broadcast the global state.
 947      */
 948     if (inputrec->eI == eiNM || inputrec->eI == eiTPI)
 949     {
 950         if (!MASTER(cr))
 951         {
 952             globalState = std::make_unique<t_state>();
 953         }
 954         broadcastStateWithoutDynamics(cr, globalState.get());
 955     }
 956
 957     /* A parallel command line option consistency check that we can
 958        only do after any threads have started. */
 959     if (!PAR(cr)
 960         && (domdecOptions.numCells[XX] > 1 || domdecOptions.numCells[YY] > 1
 961             || domdecOptions.numCells[ZZ] > 1 || domdecOptions.numPmeRanks > 0))
 962     {
 963         gmx_fatal(FARGS,
 964                   "The -dd or -npme option request a parallel simulation, "
 965 #if !GMX_MPI
 966                   "but %s was compiled without threads or MPI enabled",
 967                   output_env_get_program_display_name(oenv));
 968 #elif GMX_THREAD_MPI
 969                   "but the number of MPI-threads (option -ntmpi) is not set or is 1");
 970 #else
 971                   "but %s was not started through mpirun/mpiexec or only one rank was requested "
 972                   "through mpirun/mpiexec",
 973                   output_env_get_program_display_name(oenv));
 974 #endif
 975     }
 976
 977     if (doRerun && (EI_ENERGY_MINIMIZATION(inputrec->eI) || eiNM == inputrec->eI))
 978     {
 979         gmx_fatal(FARGS,
 980                   "The .mdp file specified an energy mininization or normal mode algorithm, and "
 981                   "these are not compatible with mdrun -rerun");
 982     }
 983
 984     if (!(EEL_PME(inputrec->coulombtype) || EVDW_PME(inputrec->vdwtype)))
 985     {
 986         if (domdecOptions.numPmeRanks > 0)
 987         {
 988             gmx_fatal_collective(FARGS, cr->mpi_comm_mysim, MASTER(cr),
 989                                  "PME-only ranks are requested, but the system does not use PME "
 990                                  "for electrostatics or LJ");
 991         }
 992
 993         domdecOptions.numPmeRanks = 0;
 994     }
 995
 996     if (useGpuForNonbonded && domdecOptions.numPmeRanks < 0)
 997     {
 998         /* With NB GPUs we don't automatically use PME-only CPU ranks. PME ranks can
 999          * improve performance with many threads per GPU, since our OpenMP
1000          * scaling is bad, but it's difficult to automate the setup.
1001          */
1002         domdecOptions.numPmeRanks = 0;
1003     }
1004     if (useGpuForPme)
1005     {
1006         if (domdecOptions.numPmeRanks < 0)
1007         {
1008             domdecOptions.numPmeRanks = 0;
1009             // TODO possibly print a note that one can opt-in for a separate PME GPU rank?
1010         }
1011         else
1012         {
1013             GMX_RELEASE_ASSERT(domdecOptions.numPmeRanks <= 1,
1014                                "PME GPU decomposition is not supported");
1015         }
1016     }
1017
1018 #if GMX_FAHCORE
1019     if (MASTER(cr))
1020     {
1021         fcRegisterSteps(inputrec->nsteps, inputrec->init_step);
1022     }
1023 #endif
1024
1025     /* NMR restraints must be initialized before load_checkpoint,
1026      * since with time averaging the history is added to t_state.
1027      * For proper consistency check we therefore need to extend
1028      * t_state here.
1029      * So the PME-only nodes (if present) will also initialize
1030      * the distance restraints.
1031      */
1032     snew(fcd, 1);
1033
1034     /* This needs to be called before read_checkpoint to extend the state */
1035     init_disres(fplog, &mtop, inputrec, cr, ms, fcd, globalState.get(), replExParams.exchangeInterval > 0);
1036
1037     init_orires(fplog, &mtop, inputrec, cr, ms, globalState.get(), &(fcd->orires));
1038
1039     auto deform = prepareBoxDeformation(globalState->box, cr, *inputrec);
1040
1041     ObservablesHistory observablesHistory = {};
1042
1043     if (startingBehavior != StartingBehavior::NewSimulation)
1044     {
1045         /* Check if checkpoint file exists before doing continuation.
1046          * This way we can use identical input options for the first and subsequent runs...
1047          */
1048         if (mdrunOptions.numStepsCommandline > -2)
1049         {
1050             /* Temporarily set the number of steps to unlimited to avoid
1051              * triggering the nsteps check in load_checkpoint().
1052              * This hack will go away soon when the -nsteps option is removed.
1053              */
1054             inputrec->nsteps = -1;
1055         }
1056
1057         load_checkpoint(opt2fn_master("-cpi", filenames.size(), filenames.data(), cr),
1058                         logFileHandle, cr, domdecOptions.numCells, inputrec, globalState.get(),
1059                         &observablesHistory, mdrunOptions.reproducible, mdModules_->notifier());
1060
1061         if (startingBehavior == StartingBehavior::RestartWithAppending && logFileHandle)
1062         {
1063             // Now we can start normal logging to the truncated log file.
1064             fplog = gmx_fio_getfp(logFileHandle);
1065             prepareLogAppending(fplog);
1066             logOwner = buildLogger(fplog, MASTER(cr));
1067             mdlog    = logOwner.logger();
1068         }
1069     }
1070
1071     if (mdrunOptions.numStepsCommandline > -2)
1072     {
1073         GMX_LOG(mdlog.info)
1074                 .asParagraph()
1075                 .appendText(
1076                         "The -nsteps functionality is deprecated, and may be removed in a future "
1077                         "version. "
1078                         "Consider using gmx convert-tpr -nsteps or changing the appropriate .mdp "
1079                         "file field.");
1080     }
1081     /* override nsteps with value set on the commandline */
1082     override_nsteps_cmdline(mdlog, mdrunOptions.numStepsCommandline, inputrec);
1083
1084     if (SIMMASTER(cr))
1085     {
1086         copy_mat(globalState->box, box);
1087     }
1088
1089     if (PAR(cr))
1090     {
1091         gmx_bcast(sizeof(box), box, cr);
1092     }
1093
1094     if (inputrec->cutoff_scheme != ecutsVERLET)
1095     {
1096         gmx_fatal(FARGS,
1097                   "This group-scheme .tpr file can no longer be run by mdrun. Please update to the "
1098                   "Verlet scheme, or use an earlier version of GROMACS if necessary.");
1099     }
1100     /* Update rlist and nstlist. */
1101     prepare_verlet_scheme(fplog, cr, inputrec, nstlist_cmdline, &mtop, box,
1102                           useGpuForNonbonded || (emulateGpuNonbonded == EmulateGpuNonbonded::Yes),
1103                           *hwinfo->cpuInfo);
1104
1105     const bool prefer1DAnd1PulseDD = (devFlags.enableGpuHaloExchange && useGpuForNonbonded);
1106     // This builder is necessary while we have multi-part construction
1107     // of DD. Before DD is constructed, we use the existence of
1108     // the builder object to indicate that further construction of DD
1109     // is needed.
1110     std::unique_ptr<DomainDecompositionBuilder> ddBuilder;
1111     if (useDomainDecomposition)
1112     {
1113         ddBuilder = std::make_unique<DomainDecompositionBuilder>(
1114                 mdlog, cr, domdecOptions, mdrunOptions, prefer1DAnd1PulseDD, mtop, *inputrec, box,
1115                 positionsFromStatePointer(globalState.get()));
1116     }
1117     else
1118     {
1119         /* PME, if used, is done on all nodes with 1D decomposition */
1120         cr->npmenodes = 0;
1121         cr->duty      = (DUTY_PP | DUTY_PME);
1122
1123         if (inputrec->ePBC == epbcSCREW)
1124         {
1125             gmx_fatal(FARGS, "pbc=screw is only implemented with domain decomposition");
1126         }
1127     }
1128
1129     // Produce the task assignment for this rank.
1130     GpuTaskAssignmentsBuilder gpuTaskAssignmentsBuilder;
1131     GpuTaskAssignments        gpuTaskAssignments = gpuTaskAssignmentsBuilder.build(
1132             gpuIdsToUse, userGpuTaskAssignment, *hwinfo, communicator, physicalNodeComm,
1133             nonbondedTarget, pmeTarget, bondedTarget, updateTarget, useGpuForNonbonded,
1134             useGpuForPme, thisRankHasDuty(cr, DUTY_PP),
1135             // TODO cr->duty & DUTY_PME should imply that a PME
1136             // algorithm is active, but currently does not.
1137             EEL_PME(inputrec->coulombtype) && thisRankHasDuty(cr, DUTY_PME));
1138
1139     // Get the device handles for the modules, nullptr when no task is assigned.
1140     gmx_device_info_t* nonbondedDeviceInfo = gpuTaskAssignments.initNonbondedDevice(cr);
1141     gmx_device_info_t* pmeDeviceInfo       = gpuTaskAssignments.initPmeDevice();
1142
1143     // TODO Initialize GPU streams here.
1144
1145     // TODO Currently this is always built, yet DD partition code
1146     // checks if it is built before using it. Probably it should
1147     // become an MDModule that is made only when another module
1148     // requires it (e.g. pull, CompEl, density fitting), so that we
1149     // don't update the local atom sets unilaterally every step.
1150     LocalAtomSetManager atomSets;
1151     if (ddBuilder)
1152     {
1153         // TODO Pass the GPU streams to ddBuilder to use in buffer
1154         // transfers (e.g. halo exchange)
1155         cr->dd = ddBuilder->build(&atomSets);
1156         // The builder's job is done, so destruct it
1157         ddBuilder.reset(nullptr);
1158         // Note that local state still does not exist yet.
1159     }
1160
1161     // The GPU update is decided here because we need to know whether the constraints or
1162     // SETTLEs can span accross the domain borders (i.e. whether or not update groups are
1163     // defined). This is only known after DD is initialized, hence decision on using GPU
1164     // update is done so late.
1165     try
1166     {
1167         const bool useUpdateGroups = cr->dd ? ddUsesUpdateGroups(*cr->dd) : false;
1168
1169         useGpuForUpdate = decideWhetherToUseGpuForUpdate(
1170                 devFlags.forceGpuUpdateDefault, useDomainDecomposition, useUpdateGroups, pmeRunMode,
1171                 domdecOptions.numPmeRanks > 0, useGpuForNonbonded, updateTarget, gpusWereDetected,
1172                 *inputrec, mtop, doEssentialDynamics, gmx_mtop_ftype_count(mtop, F_ORIRES) > 0,
1173                 replExParams.exchangeInterval > 0, doRerun, mdlog);
1174     }
1175     GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1176
1177     const bool printHostName = (cr->nnodes > 1);
1178     gpuTaskAssignments.reportGpuUsage(mdlog, printHostName, useGpuForBonded, pmeRunMode, useGpuForUpdate);
1179
1180     // If the user chose a task assignment, give them some hints
1181     // where appropriate.
1182     if (!userGpuTaskAssignment.empty())
1183     {
1184         gpuTaskAssignments.logPerformanceHints(mdlog, ssize(gpuIdsToUse));
1185     }
1186
1187     if (PAR(cr))
1188     {
1189         /* After possible communicator splitting in make_dd_communicators.
1190          * we can set up the intra/inter node communication.
1191          */
1192         gmx_setup_nodecomm(fplog, cr);
1193     }
1194
1195 #if GMX_MPI
1196     if (isMultiSim(ms))
1197     {
1198         GMX_LOG(mdlog.warning)
1199                 .asParagraph()
1200                 .appendTextFormatted(
1201                         "This is simulation %d out of %d running as a composite GROMACS\n"
1202                         "multi-simulation job. Setup for this simulation:\n",
1203                         ms->sim, ms->nsim);
1204     }
1205     GMX_LOG(mdlog.warning)
1206             .appendTextFormatted("Using %d MPI %s\n", cr->nnodes,
1207 #    if GMX_THREAD_MPI
1208                                  cr->nnodes == 1 ? "thread" : "threads"
1209 #    else
1210                                  cr->nnodes == 1 ? "process" : "processes"
1211 #    endif
1212             );
1213     fflush(stderr);
1214 #endif
1215
1216     // If mdrun -pin auto honors any affinity setting that already
1217     // exists. If so, it is nice to provide feedback about whether
1218     // that existing affinity setting was from OpenMP or something
1219     // else, so we run this code both before and after we initialize
1220     // the OpenMP support.
1221     gmx_check_thread_affinity_set(mdlog, &hw_opt, hwinfo->nthreads_hw_avail, FALSE);
1222     /* Check and update the number of OpenMP threads requested */
1223     checkAndUpdateRequestedNumOpenmpThreads(&hw_opt, *hwinfo, cr, ms, physicalNodeComm.size_,
1224                                             pmeRunMode, mtop, *inputrec);
1225
1226     gmx_omp_nthreads_init(mdlog, cr, hwinfo->nthreads_hw_avail, physicalNodeComm.size_,
1227                           hw_opt.nthreads_omp, hw_opt.nthreads_omp_pme, !thisRankHasDuty(cr, DUTY_PP));
1228
1229     // Enable FP exception detection, but not in
1230     // Release mode and not for compilers with known buggy FP
1231     // exception support (clang with any optimization) or suspected
1232     // buggy FP exception support (gcc 7.* with optimization).
1233 #if !defined NDEBUG                                                                         \
1234         && !((defined __clang__ || (defined(__GNUC__) && !defined(__ICC) && __GNUC__ == 7)) \
1235              && defined __OPTIMIZE__)
1236     const bool bEnableFPE = true;
1237 #else
1238     const bool bEnableFPE = false;
1239 #endif
1240     // FIXME - reconcile with gmx_feenableexcept() call from CommandLineModuleManager::run()
1241     if (bEnableFPE)
1242     {
1243         gmx_feenableexcept();
1244     }
1245
1246     /* Now that we know the setup is consistent, check for efficiency */
1247     check_resource_division_efficiency(hwinfo, gpuTaskAssignments.thisRankHasAnyGpuTask(),
1248                                        mdrunOptions.ntompOptionIsSet, cr, mdlog);
1249
1250     /* getting number of PP/PME threads on this MPI / tMPI rank.
1251        PME: env variable should be read only on one node to make sure it is
1252        identical everywhere;
1253      */
1254     const int numThreadsOnThisRank = thisRankHasDuty(cr, DUTY_PP) ? gmx_omp_nthreads_get(emntNonbonded)
1255                                                                   : gmx_omp_nthreads_get(emntPME);
1256     checkHardwareOversubscription(numThreadsOnThisRank, cr->nodeid, *hwinfo->hardwareTopology,
1257                                   physicalNodeComm, mdlog);
1258
1259     // Enable Peer access between GPUs where available
1260     // Only for DD, only master PP rank needs to perform setup, and only if thread MPI plus
1261     // any of the GPU communication features are active.
1262     if (DOMAINDECOMP(cr) && MASTER(cr) && thisRankHasDuty(cr, DUTY_PP) && GMX_THREAD_MPI
1263         && (devFlags.enableGpuHaloExchange || devFlags.enableGpuPmePPComm))
1264     {
1265         setupGpuDevicePeerAccess(gpuIdsToUse, mdlog);
1266     }
1267
1268     if (hw_opt.threadAffinity != ThreadAffinity::Off)
1269     {
1270         /* Before setting affinity, check whether the affinity has changed
1271          * - which indicates that probably the OpenMP library has changed it
1272          * since we first checked).
1273          */
1274         gmx_check_thread_affinity_set(mdlog, &hw_opt, hwinfo->nthreads_hw_avail, TRUE);
1275
1276         int numThreadsOnThisNode, intraNodeThreadOffset;
1277         analyzeThreadsOnThisNode(physicalNodeComm, numThreadsOnThisRank, &numThreadsOnThisNode,
1278                                  &intraNodeThreadOffset);
1279
1280         /* Set the CPU affinity */
1281         gmx_set_thread_affinity(mdlog, cr, &hw_opt, *hwinfo->hardwareTopology, numThreadsOnThisRank,
1282                                 numThreadsOnThisNode, intraNodeThreadOffset, nullptr);
1283     }
1284
1285     if (mdrunOptions.timingOptions.resetStep > -1)
1286     {
1287         GMX_LOG(mdlog.info)
1288                 .asParagraph()
1289                 .appendText(
1290                         "The -resetstep functionality is deprecated, and may be removed in a "
1291                         "future version.");
1292     }
1293     wcycle = wallcycle_init(fplog, mdrunOptions.timingOptions.resetStep, cr);
1294
1295     if (PAR(cr))
1296     {
1297         /* Master synchronizes its value of reset_counters with all nodes
1298          * including PME only nodes */
1299         int64_t reset_counters = wcycle_get_reset_counters(wcycle);
1300         gmx_bcast_sim(sizeof(reset_counters), &reset_counters, cr);
1301         wcycle_set_reset_counters(wcycle, reset_counters);
1302     }
1303
1304     // Membrane embedding must be initialized before we call init_forcerec()
1305     if (doMembed)
1306     {
1307         if (MASTER(cr))
1308         {
1309             fprintf(stderr, "Initializing membed");
1310         }
1311         /* Note that membed cannot work in parallel because mtop is
1312          * changed here. Fix this if we ever want to make it run with
1313          * multiple ranks. */
1314         membed = init_membed(fplog, filenames.size(), filenames.data(), &mtop, inputrec,
1315                              globalState.get(), cr, &mdrunOptions.checkpointOptions.period);
1316     }
1317
1318     const bool                   thisRankHasPmeGpuTask = gpuTaskAssignments.thisRankHasPmeGpuTask();
1319     std::unique_ptr<MDAtoms>     mdAtoms;
1320     std::unique_ptr<gmx_vsite_t> vsite;
1321
1322     t_nrnb nrnb;
1323     if (thisRankHasDuty(cr, DUTY_PP))
1324     {
1325         mdModulesNotifier.notify(*cr);
1326         mdModulesNotifier.notify(&atomSets);
1327         mdModulesNotifier.notify(PeriodicBoundaryConditionType{ inputrec->ePBC });
1328         mdModulesNotifier.notify(SimulationTimeStep{ inputrec->delta_t });
1329         /* Initiate forcerecord */
1330         fr                 = new t_forcerec;
1331         fr->forceProviders = mdModules_->initForceProviders();
1332         init_forcerec(fplog, mdlog, fr, fcd, inputrec, &mtop, cr, box,
1333                       opt2fn("-table", filenames.size(), filenames.data()),
1334                       opt2fn("-tablep", filenames.size(), filenames.data()),
1335                       opt2fns("-tableb", filenames.size(), filenames.data()), *hwinfo,
1336                       nonbondedDeviceInfo, useGpuForBonded,
1337                       pmeRunMode == PmeRunMode::GPU && !thisRankHasDuty(cr, DUTY_PME), pforce, wcycle);
1338
1339         // TODO Move this to happen during domain decomposition setup,
1340         // once stream and event handling works well with that.
1341         // TODO remove need to pass local stream into GPU halo exchange - Redmine #3093
1342         if (havePPDomainDecomposition(cr) && prefer1DAnd1PulseDD && is1DAnd1PulseDD(*cr->dd))
1343         {
1344             GMX_RELEASE_ASSERT(devFlags.enableGpuBufferOps,
1345                                "Must use GMX_USE_GPU_BUFFER_OPS=1 to use GMX_GPU_DD_COMMS=1");
1346             void* streamLocal =
1347                     Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, InteractionLocality::Local);
1348             void* streamNonLocal =
1349                     Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, InteractionLocality::NonLocal);
1350             GMX_LOG(mdlog.warning)
1351                     .asParagraph()
1352                     .appendTextFormatted(
1353                             "NOTE: This run uses the 'GPU halo exchange' feature, enabled by the "
1354                             "GMX_GPU_DD_COMMS environment variable.");
1355             cr->dd->gpuHaloExchange = std::make_unique<GpuHaloExchange>(
1356                     cr->dd, cr->mpi_comm_mysim, streamLocal, streamNonLocal);
1357         }
1358
1359         /* Initialize the mdAtoms structure.
1360          * mdAtoms is not filled with atom data,
1361          * as this can not be done now with domain decomposition.
1362          */
1363         mdAtoms = makeMDAtoms(fplog, mtop, *inputrec, thisRankHasPmeGpuTask);
1364         if (globalState && thisRankHasPmeGpuTask)
1365         {
1366             // The pinning of coordinates in the global state object works, because we only use
1367             // PME on GPU without DD or on a separate PME rank, and because the local state pointer
1368             // points to the global state object without DD.
1369             // FIXME: MD and EM separately set up the local state - this should happen in the same
1370             // function, which should also perform the pinning.
1371             changePinningPolicy(&globalState->x, pme_get_pinning_policy());
1372         }
1373
1374         /* Initialize the virtual site communication */
1375         vsite = initVsite(mtop, cr);
1376
1377         calc_shifts(box, fr->shift_vec);
1378
1379         /* With periodic molecules the charge groups should be whole at start up
1380          * and the virtual sites should not be far from their proper positions.
1381          */
1382         if (!inputrec->bContinuation && MASTER(cr) && !(inputrec->ePBC != epbcNONE && inputrec->bPeriodicMols))
1383         {
1384             /* Make molecules whole at start of run */
1385             if (fr->ePBC != epbcNONE)
1386             {
1387                 do_pbc_first_mtop(fplog, inputrec->ePBC, box, &mtop, globalState->x.rvec_array());
1388             }
1389             if (vsite)
1390             {
1391                 /* Correct initial vsite positions are required
1392                  * for the initial distribution in the domain decomposition
1393                  * and for the initial shell prediction.
1394                  */
1395                 constructVsitesGlobal(mtop, globalState->x);
1396             }
1397         }
1398
1399         if (EEL_PME(fr->ic->eeltype) || EVDW_PME(fr->ic->vdwtype))
1400         {
1401             ewaldcoeff_q  = fr->ic->ewaldcoeff_q;
1402             ewaldcoeff_lj = fr->ic->ewaldcoeff_lj;
1403         }
1404     }
1405     else
1406     {
1407         /* This is a PME only node */
1408
1409         GMX_ASSERT(globalState == nullptr,
1410                    "We don't need the state on a PME only rank and expect it to be unitialized");
1411
1412         ewaldcoeff_q  = calc_ewaldcoeff_q(inputrec->rcoulomb, inputrec->ewald_rtol);
1413         ewaldcoeff_lj = calc_ewaldcoeff_lj(inputrec->rvdw, inputrec->ewald_rtol_lj);
1414     }
1415
1416     gmx_pme_t* sepPmeData = nullptr;
1417     // This reference hides the fact that PME data is owned by runner on PME-only ranks and by forcerec on other ranks
1418     GMX_ASSERT(thisRankHasDuty(cr, DUTY_PP) == (fr != nullptr),
1419                "Double-checking that only PME-only ranks have no forcerec");
1420     gmx_pme_t*& pmedata = fr ? fr->pmedata : sepPmeData;
1421
1422     // TODO should live in ewald module once its testing is improved
1423     //
1424     // Later, this program could contain kernels that might be later
1425     // re-used as auto-tuning progresses, or subsequent simulations
1426     // are invoked.
1427     PmeGpuProgramStorage pmeGpuProgram;
1428     if (thisRankHasPmeGpuTask)
1429     {
1430         pmeGpuProgram = buildPmeGpuProgram(pmeDeviceInfo);
1431     }
1432
1433     /* Initiate PME if necessary,
1434      * either on all nodes or on dedicated PME nodes only. */
1435     if (EEL_PME(inputrec->coulombtype) || EVDW_PME(inputrec->vdwtype))
1436     {
1437         if (mdAtoms && mdAtoms->mdatoms())
1438         {
1439             nChargePerturbed = mdAtoms->mdatoms()->nChargePerturbed;
1440             if (EVDW_PME(inputrec->vdwtype))
1441             {
1442                 nTypePerturbed = mdAtoms->mdatoms()->nTypePerturbed;
1443             }
1444         }
1445         if (cr->npmenodes > 0)
1446         {
1447             /* The PME only nodes need to know nChargePerturbed(FEP on Q) and nTypePerturbed(FEP on LJ)*/
1448             gmx_bcast_sim(sizeof(nChargePerturbed), &nChargePerturbed, cr);
1449             gmx_bcast_sim(sizeof(nTypePerturbed), &nTypePerturbed, cr);
1450         }
1451
1452         if (thisRankHasDuty(cr, DUTY_PME))
1453         {
1454             try
1455             {
1456                 pmedata = gmx_pme_init(cr, getNumPmeDomains(cr->dd), inputrec, nChargePerturbed != 0,
1457                                        nTypePerturbed != 0, mdrunOptions.reproducible, ewaldcoeff_q,
1458                                        ewaldcoeff_lj, gmx_omp_nthreads_get(emntPME), pmeRunMode,
1459                                        nullptr, pmeDeviceInfo, pmeGpuProgram.get(), mdlog);
1460             }
1461             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1462         }
1463     }
1464
1465
1466     if (EI_DYNAMICS(inputrec->eI))
1467     {
1468         /* Turn on signal handling on all nodes */
1469         /*
1470          * (A user signal from the PME nodes (if any)
1471          * is communicated to the PP nodes.
1472          */
1473         signal_handler_install();
1474     }
1475
1476     pull_t* pull_work = nullptr;
1477     if (thisRankHasDuty(cr, DUTY_PP))
1478     {
1479         /* Assumes uniform use of the number of OpenMP threads */
1480         walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntDefault));
1481
1482         if (inputrec->bPull)
1483         {
1484             /* Initialize pull code */
1485             pull_work = init_pull(fplog, inputrec->pull, inputrec, &mtop, cr, &atomSets,
1486                                   inputrec->fepvals->init_lambda);
1487             if (inputrec->pull->bXOutAverage || inputrec->pull->bFOutAverage)
1488             {
1489                 initPullHistory(pull_work, &observablesHistory);
1490             }
1491             if (EI_DYNAMICS(inputrec->eI) && MASTER(cr))
1492             {
1493                 init_pull_output_files(pull_work, filenames.size(), filenames.data(), oenv, startingBehavior);
1494             }
1495         }
1496
1497         std::unique_ptr<EnforcedRotation> enforcedRotation;
1498         if (inputrec->bRot)
1499         {
1500             /* Initialize enforced rotation code */
1501             enforcedRotation =
1502                     init_rot(fplog, inputrec, filenames.size(), filenames.data(), cr, &atomSets,
1503                              globalState.get(), &mtop, oenv, mdrunOptions, startingBehavior);
1504         }
1505
1506         t_swap* swap = nullptr;
1507         if (inputrec->eSwapCoords != eswapNO)
1508         {
1509             /* Initialize ion swapping code */
1510             swap = init_swapcoords(fplog, inputrec,
1511                                    opt2fn_master("-swap", filenames.size(), filenames.data(), cr),
1512                                    &mtop, globalState.get(), &observablesHistory, cr, &atomSets,
1513                                    oenv, mdrunOptions, startingBehavior);
1514         }
1515
1516         /* Let makeConstraints know whether we have essential dynamics constraints. */
1517         auto constr = makeConstraints(mtop, *inputrec, pull_work, doEssentialDynamics, fplog,
1518                                       *mdAtoms->mdatoms(), cr, ms, &nrnb, wcycle, fr->bMolPBC);
1519
1520         /* Energy terms and groups */
1521         gmx_enerdata_t enerd(mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
1522                              inputrec->fepvals->n_lambda);
1523
1524         /* Kinetic energy data */
1525         gmx_ekindata_t ekind;
1526         init_ekindata(fplog, &mtop, &(inputrec->opts), &ekind, inputrec->cos_accel);
1527
1528         /* Set up interactive MD (IMD) */
1529         auto imdSession =
1530                 makeImdSession(inputrec, cr, wcycle, &enerd, ms, &mtop, mdlog,
1531                                MASTER(cr) ? globalState->x.rvec_array() : nullptr, filenames.size(),
1532                                filenames.data(), oenv, mdrunOptions.imdOptions, startingBehavior);
1533
1534         if (DOMAINDECOMP(cr))
1535         {
1536             GMX_RELEASE_ASSERT(fr, "fr was NULL while cr->duty was DUTY_PP");
1537             /* This call is not included in init_domain_decomposition mainly
1538              * because fr->cginfo_mb is set later.
1539              */
1540             dd_init_bondeds(fplog, cr->dd, &mtop, vsite.get(), inputrec,
1541                             domdecOptions.checkBondedInteractions, fr->cginfo_mb);
1542         }
1543
1544         const bool inputIsCompatibleWithModularSimulator = ModularSimulator::isInputCompatible(
1545                 false, inputrec, doRerun, vsite.get(), ms, replExParams, fcd,
1546                 static_cast<int>(filenames.size()), filenames.data(), &observablesHistory, membed);
1547
1548         const bool useModularSimulator = inputIsCompatibleWithModularSimulator
1549                                          && !(getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
1550
1551         // TODO This is not the right place to manage the lifetime of
1552         // this data structure, but currently it's the easiest way to
1553         // make it work.
1554         MdrunScheduleWorkload runScheduleWork;
1555         // Also populates the simulation constant workload description.
1556         runScheduleWork.simulationWork = createSimulationWorkload(
1557                 useGpuForNonbonded, pmeRunMode, useGpuForBonded, useGpuForUpdate,
1558                 devFlags.enableGpuBufferOps, devFlags.enableGpuHaloExchange,
1559                 devFlags.enableGpuPmePPComm, haveEwaldSurfaceContribution(*inputrec));
1560
1561         std::unique_ptr<gmx::StatePropagatorDataGpu> stateGpu;
1562         if (gpusWereDetected
1563             && ((useGpuForPme && thisRankHasDuty(cr, DUTY_PME))
1564                 || runScheduleWork.simulationWork.useGpuBufferOps))
1565         {
1566             const void* pmeStream = pme_gpu_get_device_stream(fr->pmedata);
1567             const void* localStream =
1568                     fr->nbv->gpu_nbv != nullptr
1569                             ? Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, InteractionLocality::Local)
1570                             : nullptr;
1571             const void* nonLocalStream =
1572                     fr->nbv->gpu_nbv != nullptr
1573                             ? Nbnxm::gpu_get_command_stream(fr->nbv->gpu_nbv, InteractionLocality::NonLocal)
1574                             : nullptr;
1575             const void*        deviceContext = pme_gpu_get_device_context(fr->pmedata);
1576             const int          paddingSize   = pme_gpu_get_padding_size(fr->pmedata);
1577             GpuApiCallBehavior transferKind = (inputrec->eI == eiMD && !doRerun && !useModularSimulator)
1578                                                       ? GpuApiCallBehavior::Async
1579                                                       : GpuApiCallBehavior::Sync;
1580
1581             stateGpu = std::make_unique<gmx::StatePropagatorDataGpu>(
1582                     pmeStream, localStream, nonLocalStream, deviceContext, transferKind, paddingSize, wcycle);
1583             fr->stateGpu = stateGpu.get();
1584         }
1585
1586         GMX_ASSERT(stopHandlerBuilder_, "Runner must provide StopHandlerBuilder to simulator.");
1587         SimulatorBuilder simulatorBuilder;
1588
1589         // build and run simulator object based on user-input
1590         auto simulator = simulatorBuilder.build(
1591                 inputIsCompatibleWithModularSimulator, fplog, cr, ms, mdlog,
1592                 static_cast<int>(filenames.size()), filenames.data(), oenv, mdrunOptions,
1593                 startingBehavior, vsite.get(), constr.get(),
1594                 enforcedRotation ? enforcedRotation->getLegacyEnfrot() : nullptr, deform.get(),
1595                 mdModules_->outputProvider(), mdModules_->notifier(), inputrec, imdSession.get(),
1596                 pull_work, swap, &mtop, fcd, globalState.get(), &observablesHistory, mdAtoms.get(),
1597                 &nrnb, wcycle, fr, &enerd, &ekind, &runScheduleWork, replExParams, membed,
1598                 walltime_accounting, std::move(stopHandlerBuilder_), doRerun);
1599         simulator->run();
1600
1601         if (fr->pmePpCommGpu)
1602         {
1603             // destroy object since it is no longer required. (This needs to be done while the GPU context still exists.)
1604             fr->pmePpCommGpu.reset();
1605         }
1606
1607         if (inputrec->bPull)
1608         {
1609             finish_pull(pull_work);
1610         }
1611         finish_swapcoords(swap);
1612     }
1613     else
1614     {
1615         GMX_RELEASE_ASSERT(pmedata, "pmedata was NULL while cr->duty was not DUTY_PP");
1616         /* do PME only */
1617         walltime_accounting = walltime_accounting_init(gmx_omp_nthreads_get(emntPME));
1618         gmx_pmeonly(pmedata, cr, &nrnb, wcycle, walltime_accounting, inputrec, pmeRunMode);
1619     }
1620
1621     wallcycle_stop(wcycle, ewcRUN);
1622
1623     /* Finish up, write some stuff
1624      * if rerunMD, don't write last frame again
1625      */
1626     finish_run(fplog, mdlog, cr, inputrec, &nrnb, wcycle, walltime_accounting,
1627                fr ? fr->nbv.get() : nullptr, pmedata, EI_DYNAMICS(inputrec->eI) && !isMultiSim(ms));
1628
1629     // clean up cycle counter
1630     wallcycle_destroy(wcycle);
1631
1632     // Free PME data
1633     if (pmedata)
1634     {
1635         gmx_pme_destroy(pmedata);
1636         pmedata = nullptr;
1637     }
1638
1639     // FIXME: this is only here to manually unpin mdAtoms->chargeA_ and state->x,
1640     // before we destroy the GPU context(s) in free_gpu_resources().
1641     // Pinned buffers are associated with contexts in CUDA.
1642     // As soon as we destroy GPU contexts after mdrunner() exits, these lines should go.
1643     mdAtoms.reset(nullptr);
1644     globalState.reset(nullptr);
1645     mdModules_.reset(nullptr); // destruct force providers here as they might also use the GPU
1646
1647     /* Free GPU memory and set a physical node tMPI barrier (which should eventually go away) */
1648     free_gpu_resources(fr, physicalNodeComm, hwinfo->gpu_info);
1649     free_gpu(nonbondedDeviceInfo);
1650     free_gpu(pmeDeviceInfo);
1651     done_forcerec(fr, mtop.molblock.size());
1652     sfree(fcd);
1653
1654     if (doMembed)
1655     {
1656         free_membed(membed);
1657     }
1658
1659     /* Does what it says */
1660     print_date_and_time(fplog, cr->nodeid, "Finished mdrun", gmx_gettime());
1661     walltime_accounting_destroy(walltime_accounting);
1662
1663     // Ensure log file content is written
1664     if (logFileHandle)
1665     {
1666         gmx_fio_flush(logFileHandle);
1667     }
1668
1669     /* Reset FPEs (important for unit tests) by disabling them. Assumes no
1670      * exceptions were enabled before function was called. */
1671     if (bEnableFPE)
1672     {
1673         gmx_fedisableexcept();
1674     }
1675
1676     auto rc = static_cast<int>(gmx_get_stop_condition());
1677
1678 #if GMX_THREAD_MPI
1679     /* we need to join all threads. The sub-threads join when they
1680        exit this function, but the master thread needs to be told to
1681        wait for that. */
1682     if (PAR(cr) && MASTER(cr))
1683     {
1684         tMPI_Finalize();
1685     }
1686 #endif
1687     return rc;
1688 }
1689
1690 Mdrunner::~Mdrunner()
1691 {
1692     // Clean up of the Manager.
1693     // This will end up getting called on every thread-MPI rank, which is unnecessary,
1694     // but okay as long as threads synchronize some time before adding or accessing
1695     // a new set of restraints.
1696     if (restraintManager_)
1697     {
1698         restraintManager_->clear();
1699         GMX_ASSERT(restraintManager_->countRestraints() == 0,
1700                    "restraints added during runner life time should be cleared at runner "
1701                    "destruction.");
1702     }
1703 };
1704
1705 void Mdrunner::addPotential(std::shared_ptr<gmx::IRestraintPotential> puller, const std::string& name)
1706 {
1707     GMX_ASSERT(restraintManager_, "Mdrunner must have a restraint manager.");
1708     // Not sure if this should be logged through the md logger or something else,
1709     // but it is helpful to have some sort of INFO level message sent somewhere.
1710     //    std::cout << "Registering restraint named " << name << std::endl;
1711
1712     // When multiple restraints are used, it may be wasteful to register them separately.
1713     // Maybe instead register an entire Restraint Manager as a force provider.
1714     restraintManager_->addToSpec(std::move(puller), name);
1715 }
1716
1717 Mdrunner::Mdrunner(std::unique_ptr<MDModules> mdModules) : mdModules_(std::move(mdModules)) {}
1718
1719 Mdrunner::Mdrunner(Mdrunner&&) noexcept = default;
1720
1721 //NOLINTNEXTLINE(performance-noexcept-move-constructor) working around GCC bug 58265
1722 Mdrunner& Mdrunner::operator=(Mdrunner&& /*handle*/) noexcept(BUGFREE_NOEXCEPT_STRING) = default;
1723
1724 class Mdrunner::BuilderImplementation
1725 {
1726 public:
1727     BuilderImplementation() = delete;
1728     BuilderImplementation(std::unique_ptr<MDModules> mdModules, compat::not_null<SimulationContext*> context);
1729     ~BuilderImplementation();
1730
1731     BuilderImplementation& setExtraMdrunOptions(const MdrunOptions& options,
1732                                                 real                forceWarningThreshold,
1733                                                 StartingBehavior    startingBehavior);
1734
1735     void addDomdec(const DomdecOptions& options);
1736
1737     void addVerletList(int nstlist);
1738
1739     void addReplicaExchange(const ReplicaExchangeParameters& params);
1740
1741     void addNonBonded(const char* nbpu_opt);
1742
1743     void addPME(const char* pme_opt_, const char* pme_fft_opt_);
1744
1745     void addBondedTaskAssignment(const char* bonded_opt);
1746
1747     void addUpdateTaskAssignment(const char* update_opt);
1748
1749     void addHardwareOptions(const gmx_hw_opt_t& hardwareOptions);
1750
1751     void addFilenames(ArrayRef<const t_filenm> filenames);
1752
1753     void addOutputEnvironment(gmx_output_env_t* outputEnvironment);
1754
1755     void addLogFile(t_fileio* logFileHandle);
1756
1757     void addStopHandlerBuilder(std::unique_ptr<StopHandlerBuilder> builder);
1758
1759     Mdrunner build();
1760
1761 private:
1762     // Default parameters copied from runner.h
1763     // \todo Clarify source(s) of default parameters.
1764
1765     const char* nbpu_opt_    = nullptr;
1766     const char* pme_opt_     = nullptr;
1767     const char* pme_fft_opt_ = nullptr;
1768     const char* bonded_opt_  = nullptr;
1769     const char* update_opt_  = nullptr;
1770
1771     MdrunOptions mdrunOptions_;
1772
1773     DomdecOptions domdecOptions_;
1774
1775     ReplicaExchangeParameters replicaExchangeParameters_;
1776
1777     //! Command-line override for the duration of a neighbor list with the Verlet scheme.
1778     int nstlist_ = 0;
1779
1780     //! Multisim communicator handle.
1781     gmx_multisim_t* multiSimulation_;
1782
1783     //! mdrun communicator
1784     MPI_Comm communicator_ = MPI_COMM_NULL;
1785
1786     //! Print a warning if any force is larger than this (in kJ/mol nm).
1787     real forceWarningThreshold_ = -1;
1788
1789     //! Whether the simulation will start afresh, or restart with/without appending.
1790     StartingBehavior startingBehavior_ = StartingBehavior::NewSimulation;
1791
1792     //! The modules that comprise the functionality of mdrun.
1793     std::unique_ptr<MDModules> mdModules_;
1794
1795     //! \brief Parallelism information.
1796     gmx_hw_opt_t hardwareOptions_;
1797
1798     //! filename options for simulation.
1799     ArrayRef<const t_filenm> filenames_;
1800
1801     /*! \brief Handle to output environment.
1802      *
1803      * \todo gmx_output_env_t needs lifetime management.
1804      */
1805     gmx_output_env_t* outputEnvironment_ = nullptr;
1806
1807     /*! \brief Non-owning handle to MD log file.
1808      *
1809      * \todo Context should own output facilities for client.
1810      * \todo Improve log file handle management.
1811      * \internal
1812      * Code managing the FILE* relies on the ability to set it to
1813      * nullptr to check whether the filehandle is valid.
1814      */
1815     t_fileio* logFileHandle_ = nullptr;
1816
1817     /*!
1818      * \brief Builder for simulation stop signal handler.
1819      */
1820     std::unique_ptr<StopHandlerBuilder> stopHandlerBuilder_ = nullptr;
1821 };
1822
1823 Mdrunner::BuilderImplementation::BuilderImplementation(std::unique_ptr<MDModules> mdModules,
1824                                                        compat::not_null<SimulationContext*> context) :
1825     mdModules_(std::move(mdModules))
1826 {
1827     communicator_    = context->communicator_;
1828     multiSimulation_ = context->multiSimulation_.get();
1829 }
1830
1831 Mdrunner::BuilderImplementation::~BuilderImplementation() = default;
1832
1833 Mdrunner::BuilderImplementation&
1834 Mdrunner::BuilderImplementation::setExtraMdrunOptions(const MdrunOptions&    options,
1835                                                       const real             forceWarningThreshold,
1836                                                       const StartingBehavior startingBehavior)
1837 {
1838     mdrunOptions_          = options;
1839     forceWarningThreshold_ = forceWarningThreshold;
1840     startingBehavior_      = startingBehavior;
1841     return *this;
1842 }
1843
1844 void Mdrunner::BuilderImplementation::addDomdec(const DomdecOptions& options)
1845 {
1846     domdecOptions_ = options;
1847 }
1848
1849 void Mdrunner::BuilderImplementation::addVerletList(int nstlist)
1850 {
1851     nstlist_ = nstlist;
1852 }
1853
1854 void Mdrunner::BuilderImplementation::addReplicaExchange(const ReplicaExchangeParameters& params)
1855 {
1856     replicaExchangeParameters_ = params;
1857 }
1858
1859 Mdrunner Mdrunner::BuilderImplementation::build()
1860 {
1861     auto newRunner = Mdrunner(std::move(mdModules_));
1862
1863     newRunner.mdrunOptions     = mdrunOptions_;
1864     newRunner.pforce           = forceWarningThreshold_;
1865     newRunner.startingBehavior = startingBehavior_;
1866     newRunner.domdecOptions    = domdecOptions_;
1867
1868     // \todo determine an invariant to check or confirm that all gmx_hw_opt_t objects are valid
1869     newRunner.hw_opt = hardwareOptions_;
1870
1871     // No invariant to check. This parameter exists to optionally override other behavior.
1872     newRunner.nstlist_cmdline = nstlist_;
1873
1874     newRunner.replExParams = replicaExchangeParameters_;
1875
1876     newRunner.filenames = filenames_;
1877
1878     newRunner.communicator = communicator_;
1879
1880     // nullptr is a valid value for the multisim handle
1881     newRunner.ms = multiSimulation_;
1882
1883     // \todo Clarify ownership and lifetime management for gmx_output_env_t
1884     // \todo Update sanity checking when output environment has clearly specified invariants.
1885     // Initialization and default values for oenv are not well specified in the current version.
1886     if (outputEnvironment_)
1887     {
1888         newRunner.oenv = outputEnvironment_;
1889     }
1890     else
1891     {
1892         GMX_THROW(gmx::APIError(
1893                 "MdrunnerBuilder::addOutputEnvironment() is required before build()"));
1894     }
1895
1896     newRunner.logFileHandle = logFileHandle_;
1897
1898     if (nbpu_opt_)
1899     {
1900         newRunner.nbpu_opt = nbpu_opt_;
1901     }
1902     else
1903     {
1904         GMX_THROW(gmx::APIError("MdrunnerBuilder::addNonBonded() is required before build()"));
1905     }
1906
1907     if (pme_opt_ && pme_fft_opt_)
1908     {
1909         newRunner.pme_opt     = pme_opt_;
1910         newRunner.pme_fft_opt = pme_fft_opt_;
1911     }
1912     else
1913     {
1914         GMX_THROW(gmx::APIError("MdrunnerBuilder::addElectrostatics() is required before build()"));
1915     }
1916
1917     if (bonded_opt_)
1918     {
1919         newRunner.bonded_opt = bonded_opt_;
1920     }
1921     else
1922     {
1923         GMX_THROW(gmx::APIError(
1924                 "MdrunnerBuilder::addBondedTaskAssignment() is required before build()"));
1925     }
1926
1927     if (update_opt_)
1928     {
1929         newRunner.update_opt = update_opt_;
1930     }
1931     else
1932     {
1933         GMX_THROW(gmx::APIError(
1934                 "MdrunnerBuilder::addUpdateTaskAssignment() is required before build()  "));
1935     }
1936
1937
1938     newRunner.restraintManager_ = std::make_unique<gmx::RestraintManager>();
1939
1940     if (stopHandlerBuilder_)
1941     {
1942         newRunner.stopHandlerBuilder_ = std::move(stopHandlerBuilder_);
1943     }
1944     else
1945     {
1946         newRunner.stopHandlerBuilder_ = std::make_unique<StopHandlerBuilder>();
1947     }
1948
1949     return newRunner;
1950 }
1951
1952 void Mdrunner::BuilderImplementation::addNonBonded(const char* nbpu_opt)
1953 {
1954     nbpu_opt_ = nbpu_opt;
1955 }
1956
1957 void Mdrunner::BuilderImplementation::addPME(const char* pme_opt, const char* pme_fft_opt)
1958 {
1959     pme_opt_     = pme_opt;
1960     pme_fft_opt_ = pme_fft_opt;
1961 }
1962
1963 void Mdrunner::BuilderImplementation::addBondedTaskAssignment(const char* bonded_opt)
1964 {
1965     bonded_opt_ = bonded_opt;
1966 }
1967
1968 void Mdrunner::BuilderImplementation::addUpdateTaskAssignment(const char* update_opt)
1969 {
1970     update_opt_ = update_opt;
1971 }
1972
1973 void Mdrunner::BuilderImplementation::addHardwareOptions(const gmx_hw_opt_t& hardwareOptions)
1974 {
1975     hardwareOptions_ = hardwareOptions;
1976 }
1977
1978 void Mdrunner::BuilderImplementation::addFilenames(ArrayRef<const t_filenm> filenames)
1979 {
1980     filenames_ = filenames;
1981 }
1982
1983 void Mdrunner::BuilderImplementation::addOutputEnvironment(gmx_output_env_t* outputEnvironment)
1984 {
1985     outputEnvironment_ = outputEnvironment;
1986 }
1987
1988 void Mdrunner::BuilderImplementation::addLogFile(t_fileio* logFileHandle)
1989 {
1990     logFileHandle_ = logFileHandle;
1991 }
1992
1993 void Mdrunner::BuilderImplementation::addStopHandlerBuilder(std::unique_ptr<StopHandlerBuilder> builder)
1994 {
1995     stopHandlerBuilder_ = std::move(builder);
1996 }
1997
1998 MdrunnerBuilder::MdrunnerBuilder(std::unique_ptr<MDModules>           mdModules,
1999                                  compat::not_null<SimulationContext*> context) :
2000     impl_{ std::make_unique<Mdrunner::BuilderImplementation>(std::move(mdModules), context) }
2001 {
2002 }
2003
2004 MdrunnerBuilder::~MdrunnerBuilder() = default;
2005
2006 MdrunnerBuilder& MdrunnerBuilder::addSimulationMethod(const MdrunOptions&    options,
2007                                                       real                   forceWarningThreshold,
2008                                                       const StartingBehavior startingBehavior)
2009 {
2010     impl_->setExtraMdrunOptions(options, forceWarningThreshold, startingBehavior);
2011     return *this;
2012 }
2013
2014 MdrunnerBuilder& MdrunnerBuilder::addDomainDecomposition(const DomdecOptions& options)
2015 {
2016     impl_->addDomdec(options);
2017     return *this;
2018 }
2019
2020 MdrunnerBuilder& MdrunnerBuilder::addNeighborList(int nstlist)
2021 {
2022     impl_->addVerletList(nstlist);
2023     return *this;
2024 }
2025
2026 MdrunnerBuilder& MdrunnerBuilder::addReplicaExchange(const ReplicaExchangeParameters& params)
2027 {
2028     impl_->addReplicaExchange(params);
2029     return *this;
2030 }
2031
2032 MdrunnerBuilder& MdrunnerBuilder::addNonBonded(const char* nbpu_opt)
2033 {
2034     impl_->addNonBonded(nbpu_opt);
2035     return *this;
2036 }
2037
2038 MdrunnerBuilder& MdrunnerBuilder::addElectrostatics(const char* pme_opt, const char* pme_fft_opt)
2039 {
2040     // The builder method may become more general in the future, but in this version,
2041     // parameters for PME electrostatics are both required and the only parameters
2042     // available.
2043     if (pme_opt && pme_fft_opt)
2044     {
2045         impl_->addPME(pme_opt, pme_fft_opt);
2046     }
2047     else
2048     {
2049         GMX_THROW(
2050                 gmx::InvalidInputError("addElectrostatics() arguments must be non-null pointers."));
2051     }
2052     return *this;
2053 }
2054
2055 MdrunnerBuilder& MdrunnerBuilder::addBondedTaskAssignment(const char* bonded_opt)
2056 {
2057     impl_->addBondedTaskAssignment(bonded_opt);
2058     return *this;
2059 }
2060
2061 MdrunnerBuilder& MdrunnerBuilder::addUpdateTaskAssignment(const char* update_opt)
2062 {
2063     impl_->addUpdateTaskAssignment(update_opt);
2064     return *this;
2065 }
2066
2067 Mdrunner MdrunnerBuilder::build()
2068 {
2069     return impl_->build();
2070 }
2071
2072 MdrunnerBuilder& MdrunnerBuilder::addHardwareOptions(const gmx_hw_opt_t& hardwareOptions)
2073 {
2074     impl_->addHardwareOptions(hardwareOptions);
2075     return *this;
2076 }
2077
2078 MdrunnerBuilder& MdrunnerBuilder::addFilenames(ArrayRef<const t_filenm> filenames)
2079 {
2080     impl_->addFilenames(filenames);
2081     return *this;
2082 }
2083
2084 MdrunnerBuilder& MdrunnerBuilder::addOutputEnvironment(gmx_output_env_t* outputEnvironment)
2085 {
2086     impl_->addOutputEnvironment(outputEnvironment);
2087     return *this;
2088 }
2089
2090 MdrunnerBuilder& MdrunnerBuilder::addLogFile(t_fileio* logFileHandle)
2091 {
2092     impl_->addLogFile(logFileHandle);
2093     return *this;
2094 }
2095
2096 MdrunnerBuilder& MdrunnerBuilder::addStopHandlerBuilder(std::unique_ptr<StopHandlerBuilder> builder)
2097 {
2098     impl_->addStopHandlerBuilder(std::move(builder));
2099     return *this;
2100 }
2101
2102 MdrunnerBuilder::MdrunnerBuilder(MdrunnerBuilder&&) noexcept = default;
2103
2104 MdrunnerBuilder& MdrunnerBuilder::operator=(MdrunnerBuilder&&) noexcept = default;
2105
2106 } // namespace gmx