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