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