src/gromacs/mdrun/md.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2011-2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  *
  39  * \brief Implements the integrator for normal molecular dynamics simulations
  40  *
  41  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  42  * \ingroup module_mdrun
  43  */
  44 #include "gmxpre.h"
  45
  46 #include <cinttypes>
  47 #include <cmath>
  48 #include <cstdio>
  49 #include <cstdlib>
  50
  51 #include <algorithm>
  52 #include <memory>
  53 #include <numeric>
  54
  55 #include "gromacs/applied_forces/awh/awh.h"
  56 #include "gromacs/commandline/filenm.h"
  57 #include "gromacs/domdec/collect.h"
  58 #include "gromacs/domdec/dlbtiming.h"
  59 #include "gromacs/domdec/domdec.h"
  60 #include "gromacs/domdec/domdec_network.h"
  61 #include "gromacs/domdec/domdec_struct.h"
  62 #include "gromacs/domdec/gpuhaloexchange.h"
  63 #include "gromacs/domdec/mdsetup.h"
  64 #include "gromacs/domdec/partition.h"
  65 #include "gromacs/essentialdynamics/edsam.h"
  66 #include "gromacs/ewald/pme_load_balancing.h"
  67 #include "gromacs/ewald/pme_pp.h"
  68 #include "gromacs/fileio/trxio.h"
  69 #include "gromacs/gmxlib/network.h"
  70 #include "gromacs/gmxlib/nrnb.h"
  71 #include "gromacs/gpu_utils/device_stream_manager.h"
  72 #include "gromacs/gpu_utils/gpu_utils.h"
  73 #include "gromacs/imd/imd.h"
  74 #include "gromacs/listed_forces/listed_forces.h"
  75 #include "gromacs/math/functions.h"
  76 #include "gromacs/math/invertmatrix.h"
  77 #include "gromacs/math/vec.h"
  78 #include "gromacs/math/vectypes.h"
  79 #include "gromacs/mdlib/checkpointhandler.h"
  80 #include "gromacs/mdlib/compute_io.h"
  81 #include "gromacs/mdlib/constr.h"
  82 #include "gromacs/mdlib/coupling.h"
  83 #include "gromacs/mdlib/ebin.h"
  84 #include "gromacs/mdlib/enerdata_utils.h"
  85 #include "gromacs/mdlib/energyoutput.h"
  86 #include "gromacs/mdlib/expanded.h"
  87 #include "gromacs/mdlib/force.h"
  88 #include "gromacs/mdlib/force_flags.h"
  89 #include "gromacs/mdlib/forcerec.h"
  90 #include "gromacs/mdlib/freeenergyparameters.h"
  91 #include "gromacs/mdlib/md_support.h"
  92 #include "gromacs/mdlib/mdatoms.h"
  93 #include "gromacs/mdlib/mdoutf.h"
  94 #include "gromacs/mdlib/membed.h"
  95 #include "gromacs/mdlib/resethandler.h"
  96 #include "gromacs/mdlib/sighandler.h"
  97 #include "gromacs/mdlib/simulationsignal.h"
  98 #include "gromacs/mdlib/stat.h"
  99 #include "gromacs/mdlib/stophandler.h"
 100 #include "gromacs/mdlib/tgroup.h"
 101 #include "gromacs/mdlib/trajectory_writing.h"
 102 #include "gromacs/mdlib/update.h"
 103 #include "gromacs/mdlib/update_constrain_gpu.h"
 104 #include "gromacs/mdlib/vcm.h"
 105 #include "gromacs/mdlib/vsite.h"
 106 #include "gromacs/mdrunutility/handlerestart.h"
 107 #include "gromacs/mdrunutility/multisim.h"
 108 #include "gromacs/mdrunutility/printtime.h"
 109 #include "gromacs/mdtypes/awh_history.h"
 110 #include "gromacs/mdtypes/awh_params.h"
 111 #include "gromacs/mdtypes/commrec.h"
 112 #include "gromacs/mdtypes/df_history.h"
 113 #include "gromacs/mdtypes/energyhistory.h"
 114 #include "gromacs/mdtypes/fcdata.h"
 115 #include "gromacs/mdtypes/forcebuffers.h"
 116 #include "gromacs/mdtypes/forcerec.h"
 117 #include "gromacs/mdtypes/group.h"
 118 #include "gromacs/mdtypes/inputrec.h"
 119 #include "gromacs/mdtypes/interaction_const.h"
 120 #include "gromacs/mdtypes/md_enums.h"
 121 #include "gromacs/mdtypes/mdatom.h"
 122 #include "gromacs/mdtypes/mdrunoptions.h"
 123 #include "gromacs/mdtypes/multipletimestepping.h"
 124 #include "gromacs/mdtypes/observableshistory.h"
 125 #include "gromacs/mdtypes/pullhistory.h"
 126 #include "gromacs/mdtypes/simulation_workload.h"
 127 #include "gromacs/mdtypes/state.h"
 128 #include "gromacs/mdtypes/state_propagator_data_gpu.h"
 129 #include "gromacs/modularsimulator/energydata.h"
 130 #include "gromacs/nbnxm/gpu_data_mgmt.h"
 131 #include "gromacs/nbnxm/nbnxm.h"
 132 #include "gromacs/pbcutil/pbc.h"
 133 #include "gromacs/pulling/output.h"
 134 #include "gromacs/pulling/pull.h"
 135 #include "gromacs/swap/swapcoords.h"
 136 #include "gromacs/timing/wallcycle.h"
 137 #include "gromacs/timing/walltime_accounting.h"
 138 #include "gromacs/topology/atoms.h"
 139 #include "gromacs/topology/idef.h"
 140 #include "gromacs/topology/mtop_util.h"
 141 #include "gromacs/topology/topology.h"
 142 #include "gromacs/trajectory/trajectoryframe.h"
 143 #include "gromacs/utility/basedefinitions.h"
 144 #include "gromacs/utility/cstringutil.h"
 145 #include "gromacs/utility/fatalerror.h"
 146 #include "gromacs/utility/logger.h"
 147 #include "gromacs/utility/real.h"
 148 #include "gromacs/utility/smalloc.h"
 149
 150 #include "legacysimulator.h"
 151 #include "replicaexchange.h"
 152 #include "shellfc.h"
 153
 154 using gmx::SimulationSignaller;
 155
 156 void gmx::LegacySimulator::do_md()
 157 {
 158     // TODO Historically, the EM and MD "integrators" used different
 159     // names for the t_inputrec *parameter, but these must have the
 160     // same name, now that it's a member of a struct. We use this ir
 161     // alias to avoid a large ripple of nearly useless changes.
 162     // t_inputrec is being replaced by IMdpOptionsProvider, so this
 163     // will go away eventually.
 164     t_inputrec*  ir = inputrec;
 165     int64_t      step, step_rel;
 166     double       t, t0 = ir->init_t;
 167     gmx_bool     bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
 168     gmx_bool     bNS = FALSE, bNStList, bStopCM, bFirstStep, bInitStep, bLastStep = FALSE;
 169     gmx_bool     bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
 170     gmx_bool     do_ene, do_log, do_verbose;
 171     gmx_bool     bMasterState;
 172     unsigned int force_flags;
 173     tensor force_vir = { { 0 } }, shake_vir = { { 0 } }, total_vir = { { 0 } }, pres = { { 0 } };
 174     int    i, m;
 175     rvec   mu_tot;
 176     matrix pressureCouplingMu, M;
 177     gmx_repl_ex_t     repl_ex = nullptr;
 178     gmx_global_stat_t gstat;
 179     gmx_shellfc_t*    shellfc;
 180     gmx_bool          bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
 181     gmx_bool          bTemp, bPres, bTrotter;
 182     real              dvdl_constr;
 183     std::vector<RVec> cbuf;
 184     matrix            lastbox;
 185     int               lamnew = 0;
 186     /* for FEP */
 187     int       nstfep = 0;
 188     double    cycles;
 189     real      saved_conserved_quantity = 0;
 190     real      last_ekin                = 0;
 191     t_extmass MassQ;
 192     char      sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
 193
 194     /* PME load balancing data for GPU kernels */
 195     gmx_bool bPMETune         = FALSE;
 196     gmx_bool bPMETunePrinting = FALSE;
 197
 198     bool bInteractiveMDstep = false;
 199
 200     /* Domain decomposition could incorrectly miss a bonded
 201        interaction, but checking for that requires a global
 202        communication stage, which does not otherwise happen in DD
 203        code. So we do that alongside the first global energy reduction
 204        after a new DD is made. These variables handle whether the
 205        check happens, and the result it returns. */
 206     bool shouldCheckNumberOfBondedInteractions = false;
 207     int  totalNumberOfBondedInteractions       = -1;
 208
 209     SimulationSignals signals;
 210     // Most global communnication stages don't propagate mdrun
 211     // signals, and will use this object to achieve that.
 212     SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
 213
 214     if (!mdrunOptions.writeConfout)
 215     {
 216         // This is on by default, and the main known use case for
 217         // turning it off is for convenience in benchmarking, which is
 218         // something that should not show up in the general user
 219         // interface.
 220         GMX_LOG(mdlog.info)
 221                 .asParagraph()
 222                 .appendText(
 223                         "The -noconfout functionality is deprecated, and may be removed in a "
 224                         "future version.");
 225     }
 226
 227     /* md-vv uses averaged full step velocities for T-control
 228        md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
 229        md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
 230     bTrotter = (EI_VV(ir->eI)
 231                 && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
 232
 233     const bool bRerunMD = false;
 234
 235     int nstglobalcomm = computeGlobalCommunicationPeriod(mdlog, ir, cr);
 236     bGStatEveryStep   = (nstglobalcomm == 1);
 237
 238     const SimulationGroups* groups = &top_global->groups;
 239
 240     std::unique_ptr<EssentialDynamics> ed = nullptr;
 241     if (opt2bSet("-ei", nfile, fnm))
 242     {
 243         /* Initialize essential dynamics sampling */
 244         ed = init_edsam(mdlog, opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm), top_global,
 245                         ir, cr, constr, state_global, observablesHistory, oenv, startingBehavior);
 246     }
 247     else if (observablesHistory->edsamHistory)
 248     {
 249         gmx_fatal(FARGS,
 250                   "The checkpoint is from a run with essential dynamics sampling, "
 251                   "but the current run did not specify the -ei option. "
 252                   "Either specify the -ei option to mdrun, or do not use this checkpoint file.");
 253     }
 254
 255     int*                fep_state = MASTER(cr) ? &state_global->fep_state : nullptr;
 256     gmx::ArrayRef<real> lambda    = MASTER(cr) ? state_global->lambda : gmx::ArrayRef<real>();
 257     initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda);
 258     Update     upd(*ir, deform);
 259     const bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
 260     const bool useReplicaExchange   = (replExParams.exchangeInterval > 0);
 261
 262     const t_fcdata& fcdata = *fr->fcdata;
 263
 264     bool simulationsShareState = false;
 265     int  nstSignalComm         = nstglobalcomm;
 266     {
 267         // TODO This implementation of ensemble orientation restraints is nasty because
 268         // a user can't just do multi-sim with single-sim orientation restraints.
 269         bool usingEnsembleRestraints =
 270                 (fcdata.disres->nsystems > 1) || ((ms != nullptr) && (fcdata.orires->nr != 0));
 271         bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
 272
 273         // Replica exchange, ensemble restraints and AWH need all
 274         // simulations to remain synchronized, so they need
 275         // checkpoints and stop conditions to act on the same step, so
 276         // the propagation of such signals must take place between
 277         // simulations, not just within simulations.
 278         // TODO: Make algorithm initializers set these flags.
 279         simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
 280
 281         if (simulationsShareState)
 282         {
 283             // Inter-simulation signal communication does not need to happen
 284             // often, so we use a minimum of 200 steps to reduce overhead.
 285             const int c_minimumInterSimulationSignallingInterval = 200;
 286             nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1) / nstglobalcomm)
 287                             * nstglobalcomm;
 288         }
 289     }
 290
 291     if (startingBehavior != StartingBehavior::RestartWithAppending)
 292     {
 293         pleaseCiteCouplingAlgorithms(fplog, *ir);
 294     }
 295     gmx_mdoutf* outf =
 296             init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, ir,
 297                         top_global, oenv, wcycle, startingBehavior, simulationsShareState, ms);
 298     gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work,
 299                                    mdoutf_get_fp_dhdl(outf), false, startingBehavior,
 300                                    simulationsShareState, mdModulesNotifier);
 301
 302     gstat = global_stat_init(ir);
 303
 304     const auto& simulationWork     = runScheduleWork->simulationWork;
 305     const bool  useGpuForPme       = simulationWork.useGpuPme;
 306     const bool  useGpuForNonbonded = simulationWork.useGpuNonbonded;
 307     const bool  useGpuForBufferOps = simulationWork.useGpuBufferOps;
 308     const bool  useGpuForUpdate    = simulationWork.useGpuUpdate;
 309
 310     /* Check for polarizable models and flexible constraints */
 311     shellfc = init_shell_flexcon(fplog, top_global, constr ? constr->numFlexibleConstraints() : 0,
 312                                  ir->nstcalcenergy, DOMAINDECOMP(cr), useGpuForPme);
 313
 314     {
 315         double io = compute_io(ir, top_global->natoms, *groups, energyOutput.numEnergyTerms(), 1);
 316         if ((io > 2000) && MASTER(cr))
 317         {
 318             fprintf(stderr, "\nWARNING: This run will generate roughly %.0f Mb of data\n\n", io);
 319         }
 320     }
 321
 322     // Local state only becomes valid now.
 323     std::unique_ptr<t_state> stateInstance;
 324     t_state*                 state;
 325
 326     gmx_localtop_t top(top_global->ffparams);
 327
 328     auto mdatoms = mdAtoms->mdatoms();
 329
 330     ForceBuffers f(fr->useMts, ((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
 331                                        ? PinningPolicy::PinnedIfSupported
 332                                        : PinningPolicy::CannotBePinned);
 333     if (DOMAINDECOMP(cr))
 334     {
 335         stateInstance = std::make_unique<t_state>();
 336         state         = stateInstance.get();
 337         dd_init_local_state(cr->dd, state_global, state);
 338
 339         /* Distribute the charge groups over the nodes from the master node */
 340         dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1, state_global, *top_global, ir,
 341                             imdSession, pull_work, state, &f, mdAtoms, &top, fr, vsite, constr,
 342                             nrnb, nullptr, FALSE);
 343         shouldCheckNumberOfBondedInteractions = true;
 344         upd.setNumAtoms(state->natoms);
 345     }
 346     else
 347     {
 348         state_change_natoms(state_global, state_global->natoms);
 349         /* Copy the pointer to the global state */
 350         state = state_global;
 351
 352         /* Generate and initialize new topology */
 353         mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &f, mdAtoms, constr, vsite, shellfc);
 354
 355         upd.setNumAtoms(state->natoms);
 356     }
 357
 358     std::unique_ptr<UpdateConstrainGpu> integrator;
 359
 360     StatePropagatorDataGpu* stateGpu = fr->stateGpu;
 361
 362     // TODO: the assertions below should be handled by UpdateConstraintsBuilder.
 363     if (useGpuForUpdate)
 364     {
 365         GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr) || ddUsesUpdateGroups(*cr->dd) || constr == nullptr
 366                                    || constr->numConstraintsTotal() == 0,
 367                            "Constraints in domain decomposition are only supported with update "
 368                            "groups if using GPU update.\n");
 369         GMX_RELEASE_ASSERT(ir->eConstrAlg != econtSHAKE || constr == nullptr
 370                                    || constr->numConstraintsTotal() == 0,
 371                            "SHAKE is not supported with GPU update.");
 372         GMX_RELEASE_ASSERT(useGpuForPme || (useGpuForNonbonded && simulationWork.useGpuBufferOps),
 373                            "Either PME or short-ranged non-bonded interaction tasks must run on "
 374                            "the GPU to use GPU update.\n");
 375         GMX_RELEASE_ASSERT(ir->eI == eiMD,
 376                            "Only the md integrator is supported with the GPU update.\n");
 377         GMX_RELEASE_ASSERT(
 378                 ir->etc != etcNOSEHOOVER,
 379                 "Nose-Hoover temperature coupling is not supported with the GPU update.\n");
 380         GMX_RELEASE_ASSERT(
 381                 ir->epc == epcNO || ir->epc == epcPARRINELLORAHMAN || ir->epc == epcBERENDSEN
 382                         || ir->epc == epcCRESCALE,
 383                 "Only Parrinello-Rahman, Berendsen, and C-rescale pressure coupling are supported "
 384                 "with the GPU update.\n");
 385         GMX_RELEASE_ASSERT(!mdatoms->haveVsites,
 386                            "Virtual sites are not supported with the GPU update.\n");
 387         GMX_RELEASE_ASSERT(ed == nullptr,
 388                            "Essential dynamics is not supported with the GPU update.\n");
 389         GMX_RELEASE_ASSERT(!ir->bPull || !pull_have_constraint(*ir->pull),
 390                            "Constraints pulling is not supported with the GPU update.\n");
 391         GMX_RELEASE_ASSERT(fcdata.orires->nr == 0,
 392                            "Orientation restraints are not supported with the GPU update.\n");
 393         GMX_RELEASE_ASSERT(
 394                 ir->efep == efepNO
 395                         || (!haveFepPerturbedMasses(*top_global) && !havePerturbedConstraints(*top_global)),
 396                 "Free energy perturbation of masses and constraints are not supported with the GPU "
 397                 "update.");
 398
 399         if (constr != nullptr && constr->numConstraintsTotal() > 0)
 400         {
 401             GMX_LOG(mdlog.info)
 402                     .asParagraph()
 403                     .appendText("Updating coordinates and applying constraints on the GPU.");
 404         }
 405         else
 406         {
 407             GMX_LOG(mdlog.info).asParagraph().appendText("Updating coordinates on the GPU.");
 408         }
 409         GMX_RELEASE_ASSERT(fr->deviceStreamManager != nullptr,
 410                            "Device stream manager should be initialized in order to use GPU "
 411                            "update-constraints.");
 412         GMX_RELEASE_ASSERT(
 413                 fr->deviceStreamManager->streamIsValid(gmx::DeviceStreamType::UpdateAndConstraints),
 414                 "Update stream should be initialized in order to use GPU "
 415                 "update-constraints.");
 416         integrator = std::make_unique<UpdateConstrainGpu>(
 417                 *ir, *top_global, fr->deviceStreamManager->context(),
 418                 fr->deviceStreamManager->stream(gmx::DeviceStreamType::UpdateAndConstraints),
 419                 stateGpu->xUpdatedOnDevice(), wcycle);
 420
 421         integrator->setPbc(PbcType::Xyz, state->box);
 422     }
 423
 424     if (useGpuForPme || (useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
 425     {
 426         changePinningPolicy(&state->x, PinningPolicy::PinnedIfSupported);
 427     }
 428     if (useGpuForUpdate)
 429     {
 430         changePinningPolicy(&state->v, PinningPolicy::PinnedIfSupported);
 431     }
 432
 433     // NOTE: The global state is no longer used at this point.
 434     // But state_global is still used as temporary storage space for writing
 435     // the global state to file and potentially for replica exchange.
 436     // (Global topology should persist.)
 437
 438     update_mdatoms(mdatoms, state->lambda[efptMASS]);
 439
 440     if (ir->bExpanded)
 441     {
 442         /* Check nstexpanded here, because the grompp check was broken */
 443         if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
 444         {
 445             gmx_fatal(FARGS,
 446                       "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
 447         }
 448         init_expanded_ensemble(startingBehavior != StartingBehavior::NewSimulation, ir, state->dfhist);
 449     }
 450
 451     if (MASTER(cr))
 452     {
 453         EnergyData::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput);
 454     }
 455
 456     preparePrevStepPullCom(ir, pull_work, mdatoms->massT, state, state_global, cr,
 457                            startingBehavior != StartingBehavior::NewSimulation);
 458
 459     // TODO: Remove this by converting AWH into a ForceProvider
 460     auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms,
 461                                 startingBehavior != StartingBehavior::NewSimulation,
 462                                 shellfc != nullptr, opt2fn("-awh", nfile, fnm), pull_work);
 463
 464     if (useReplicaExchange && MASTER(cr))
 465     {
 466         repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir, replExParams);
 467     }
 468     /* PME tuning is only supported in the Verlet scheme, with PME for
 469      * Coulomb. It is not supported with only LJ PME. */
 470     bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) && !mdrunOptions.reproducible
 471                 && ir->cutoff_scheme != ecutsGROUP);
 472
 473     pme_load_balancing_t* pme_loadbal = nullptr;
 474     if (bPMETune)
 475     {
 476         pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box, *fr->ic, *fr->nbv, fr->pmedata,
 477                          fr->nbv->useGpu());
 478     }
 479
 480     if (!ir->bContinuation)
 481     {
 482         if (state->flags & (1U << estV))
 483         {
 484             auto v = makeArrayRef(state->v);
 485             /* Set the velocities of vsites, shells and frozen atoms to zero */
 486             for (i = 0; i < mdatoms->homenr; i++)
 487             {
 488                 if (mdatoms->ptype[i] == eptVSite || mdatoms->ptype[i] == eptShell)
 489                 {
 490                     clear_rvec(v[i]);
 491                 }
 492                 else if (mdatoms->cFREEZE)
 493                 {
 494                     for (m = 0; m < DIM; m++)
 495                     {
 496                         if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 497                         {
 498                             v[i][m] = 0;
 499                         }
 500                     }
 501                 }
 502             }
 503         }
 504
 505         if (constr)
 506         {
 507             /* Constrain the initial coordinates and velocities */
 508             do_constrain_first(fplog, constr, ir, mdatoms->nr, mdatoms->homenr,
 509                                state->x.arrayRefWithPadding(), state->v.arrayRefWithPadding(),
 510                                state->box, state->lambda[efptBONDED]);
 511         }
 512         if (vsite)
 513         {
 514             /* Construct the virtual sites for the initial configuration */
 515             vsite->construct(state->x, ir->delta_t, {}, state->box);
 516         }
 517     }
 518
 519     if (ir->efep != efepNO)
 520     {
 521         /* Set free energy calculation frequency as the greatest common
 522          * denominator of nstdhdl and repl_ex_nst. */
 523         nstfep = ir->fepvals->nstdhdl;
 524         if (ir->bExpanded)
 525         {
 526             nstfep = std::gcd(ir->expandedvals->nstexpanded, nstfep);
 527         }
 528         if (useReplicaExchange)
 529         {
 530             nstfep = std::gcd(replExParams.exchangeInterval, nstfep);
 531         }
 532         if (ir->bDoAwh)
 533         {
 534             nstfep = std::gcd(ir->awhParams->nstSampleCoord, nstfep);
 535         }
 536     }
 537
 538     /* Be REALLY careful about what flags you set here. You CANNOT assume
 539      * this is the first step, since we might be restarting from a checkpoint,
 540      * and in that case we should not do any modifications to the state.
 541      */
 542     bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
 543
 544     // When restarting from a checkpoint, it can be appropriate to
 545     // initialize ekind from quantities in the checkpoint. Otherwise,
 546     // compute_globals must initialize ekind before the simulation
 547     // starts/restarts. However, only the master rank knows what was
 548     // found in the checkpoint file, so we have to communicate in
 549     // order to coordinate the restart.
 550     //
 551     // TODO Consider removing this communication if/when checkpoint
 552     // reading directly follows .tpr reading, because all ranks can
 553     // agree on hasReadEkinState at that time.
 554     bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
 555     if (PAR(cr))
 556     {
 557         gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr->mpi_comm_mygroup);
 558     }
 559     if (hasReadEkinState)
 560     {
 561         restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 562     }
 563
 564     unsigned int cglo_flags =
 565             (CGLO_TEMPERATURE | CGLO_GSTAT | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
 566              | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0) | (hasReadEkinState ? CGLO_READEKIN : 0));
 567
 568     bSumEkinhOld = FALSE;
 569
 570     t_vcm vcm(top_global->groups, *ir);
 571     reportComRemovalInfo(fplog, vcm);
 572
 573     /* To minimize communication, compute_globals computes the COM velocity
 574      * and the kinetic energy for the velocities without COM motion removed.
 575      * Thus to get the kinetic energy without the COM contribution, we need
 576      * to call compute_globals twice.
 577      */
 578     for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
 579     {
 580         unsigned int cglo_flags_iteration = cglo_flags;
 581         if (bStopCM && cgloIteration == 0)
 582         {
 583             cglo_flags_iteration |= CGLO_STOPCM;
 584             cglo_flags_iteration &= ~CGLO_TEMPERATURE;
 585         }
 586         compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
 587                         makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, nullptr,
 588                         enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
 589                         state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
 590                         cglo_flags_iteration
 591                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
 592                                                                          : 0));
 593         if (cglo_flags_iteration & CGLO_STOPCM)
 594         {
 595             /* At initialization, do not pass x with acceleration-correction mode
 596              * to avoid (incorrect) correction of the initial coordinates.
 597              */
 598             auto x = (vcm.mode == ecmLINEAR_ACCELERATION_CORRECTION) ? ArrayRef<RVec>()
 599                                                                      : makeArrayRef(state->x);
 600             process_and_stopcm_grp(fplog, &vcm, *mdatoms, x, makeArrayRef(state->v));
 601             inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 602         }
 603     }
 604     checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global, &top,
 605                                     makeConstArrayRef(state->x), state->box,
 606                                     &shouldCheckNumberOfBondedInteractions);
 607     if (ir->eI == eiVVAK)
 608     {
 609         /* a second call to get the half step temperature initialized as well */
 610         /* we do the same call as above, but turn the pressure off -- internally to
 611            compute_globals, this is recognized as a velocity verlet half-step
 612            kinetic energy calculation.  This minimized excess variables, but
 613            perhaps loses some logic?*/
 614
 615         compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
 616                         makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, nullptr,
 617                         enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
 618                         state->box, nullptr, &bSumEkinhOld, cglo_flags & ~CGLO_PRESSURE);
 619     }
 620
 621     /* Calculate the initial half step temperature, and save the ekinh_old */
 622     if (startingBehavior == StartingBehavior::NewSimulation)
 623     {
 624         for (i = 0; (i < ir->opts.ngtc); i++)
 625         {
 626             copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
 627         }
 628     }
 629
 630     /* need to make an initiation call to get the Trotter variables set, as well as other constants
 631        for non-trotter temperature control */
 632     auto trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
 633
 634     if (MASTER(cr))
 635     {
 636         if (!ir->bContinuation)
 637         {
 638             if (constr && ir->eConstrAlg == econtLINCS)
 639             {
 640                 fprintf(fplog, "RMS relative constraint deviation after constraining: %.2e\n",
 641                         constr->rmsd());
 642             }
 643             if (EI_STATE_VELOCITY(ir->eI))
 644             {
 645                 real temp = enerd->term[F_TEMP];
 646                 if (ir->eI != eiVV)
 647                 {
 648                     /* Result of Ekin averaged over velocities of -half
 649                      * and +half step, while we only have -half step here.
 650                      */
 651                     temp *= 2;
 652                 }
 653                 fprintf(fplog, "Initial temperature: %g K\n", temp);
 654             }
 655         }
 656
 657         char tbuf[20];
 658         fprintf(stderr, "starting mdrun '%s'\n", *(top_global->name));
 659         if (ir->nsteps >= 0)
 660         {
 661             sprintf(tbuf, "%8.1f", (ir->init_step + ir->nsteps) * ir->delta_t);
 662         }
 663         else
 664         {
 665             sprintf(tbuf, "%s", "infinite");
 666         }
 667         if (ir->init_step > 0)
 668         {
 669             fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 670                     gmx_step_str(ir->init_step + ir->nsteps, sbuf), tbuf,
 671                     gmx_step_str(ir->init_step, sbuf2), ir->init_step * ir->delta_t);
 672         }
 673         else
 674         {
 675             fprintf(stderr, "%s steps, %s ps.\n", gmx_step_str(ir->nsteps, sbuf), tbuf);
 676         }
 677         fprintf(fplog, "\n");
 678     }
 679
 680     walltime_accounting_start_time(walltime_accounting);
 681     wallcycle_start(wcycle, ewcRUN);
 682     print_start(fplog, cr, walltime_accounting, "mdrun");
 683
 684     /***********************************************************
 685      *
 686      *             Loop over MD steps
 687      *
 688      ************************************************************/
 689
 690     bFirstStep = TRUE;
 691     /* Skip the first Nose-Hoover integration when we get the state from tpx */
 692     bInitStep        = startingBehavior == StartingBehavior::NewSimulation || EI_VV(ir->eI);
 693     bSumEkinhOld     = FALSE;
 694     bExchanged       = FALSE;
 695     bNeedRepartition = FALSE;
 696
 697     step     = ir->init_step;
 698     step_rel = 0;
 699
 700     auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
 701             compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
 702             MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
 703             mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
 704
 705     auto checkpointHandler = std::make_unique<CheckpointHandler>(
 706             compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]), simulationsShareState,
 707             ir->nstlist == 0, MASTER(cr), mdrunOptions.writeConfout,
 708             mdrunOptions.checkpointOptions.period);
 709
 710     const bool resetCountersIsLocal = true;
 711     auto       resetHandler         = std::make_unique<ResetHandler>(
 712             compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]),
 713             !resetCountersIsLocal, ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
 714             mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
 715
 716     const DDBalanceRegionHandler ddBalanceRegionHandler(cr);
 717
 718     if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
 719     {
 720         logInitialMultisimStatus(ms, cr, mdlog, simulationsShareState, ir->nsteps, ir->init_step);
 721     }
 722
 723     /* and stop now if we should */
 724     bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
 725     while (!bLastStep)
 726     {
 727
 728         /* Determine if this is a neighbor search step */
 729         bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
 730
 731         if (bPMETune && bNStList)
 732         {
 733             // This has to be here because PME load balancing is called so early.
 734             // TODO: Move to after all booleans are defined.
 735             if (useGpuForUpdate && !bFirstStep)
 736             {
 737                 stateGpu->copyCoordinatesFromGpu(state->x, AtomLocality::Local);
 738                 stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
 739             }
 740             /* PME grid + cut-off optimization with GPUs or PME nodes */
 741             pme_loadbal_do(pme_loadbal, cr, (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
 742                            fplog, mdlog, *ir, fr, state->box, state->x, wcycle, step, step_rel,
 743                            &bPMETunePrinting, simulationWork.useGpuPmePpCommunication);
 744         }
 745
 746         wallcycle_start(wcycle, ewcSTEP);
 747
 748         bLastStep = (step_rel == ir->nsteps);
 749         t         = t0 + step * ir->delta_t;
 750
 751         // TODO Refactor this, so that nstfep does not need a default value of zero
 752         if (ir->efep != efepNO || ir->bSimTemp)
 753         {
 754             /* find and set the current lambdas */
 755             state->lambda = currentLambdas(step, *(ir->fepvals), state->fep_state);
 756
 757             bDoDHDL     = do_per_step(step, ir->fepvals->nstdhdl);
 758             bDoFEP      = ((ir->efep != efepNO) && do_per_step(step, nstfep));
 759             bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded) && (ir->bExpanded)
 760                            && (!bFirstStep));
 761         }
 762
 763         bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep
 764                      && do_per_step(step, replExParams.exchangeInterval));
 765
 766         if (doSimulatedAnnealing)
 767         {
 768             update_annealing_target_temp(ir, t, &upd);
 769         }
 770
 771         /* Stop Center of Mass motion */
 772         bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
 773
 774         /* Determine whether or not to do Neighbour Searching */
 775         bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
 776
 777         /* Note that the stopHandler will cause termination at nstglobalcomm
 778          * steps. Since this concides with nstcalcenergy, nsttcouple and/or
 779          * nstpcouple steps, we have computed the half-step kinetic energy
 780          * of the previous step and can always output energies at the last step.
 781          */
 782         bLastStep = bLastStep || stopHandler->stoppingAfterCurrentStep(bNS);
 783
 784         /* do_log triggers energy and virial calculation. Because this leads
 785          * to different code paths, forces can be different. Thus for exact
 786          * continuation we should avoid extra log output.
 787          * Note that the || bLastStep can result in non-exact continuation
 788          * beyond the last step. But we don't consider that to be an issue.
 789          */
 790         do_log     = (do_per_step(step, ir->nstlog)
 791                   || (bFirstStep && startingBehavior == StartingBehavior::NewSimulation) || bLastStep);
 792         do_verbose = mdrunOptions.verbose
 793                      && (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
 794
 795         if (useGpuForUpdate && !bFirstStep && bNS)
 796         {
 797             // Copy velocities from the GPU on search steps to keep a copy on host (device buffers are reinitialized).
 798             stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
 799             stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
 800             // Copy coordinate from the GPU when needed at the search step.
 801             // NOTE: The cases when coordinates needed on CPU for force evaluation are handled in sim_utils.
 802             // NOTE: If the coordinates are to be written into output file they are also copied separately before the output.
 803             stateGpu->copyCoordinatesFromGpu(state->x, AtomLocality::Local);
 804             stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
 805         }
 806
 807         if (bNS && !(bFirstStep && ir->bContinuation))
 808         {
 809             bMasterState = FALSE;
 810             /* Correct the new box if it is too skewed */
 811             if (inputrecDynamicBox(ir))
 812             {
 813                 if (correct_box(fplog, step, state->box))
 814                 {
 815                     bMasterState = TRUE;
 816                     // If update is offloaded, it should be informed about the box size change
 817                     if (useGpuForUpdate)
 818                     {
 819                         integrator->setPbc(PbcType::Xyz, state->box);
 820                     }
 821                 }
 822             }
 823             if (DOMAINDECOMP(cr) && bMasterState)
 824             {
 825                 dd_collect_state(cr->dd, state, state_global);
 826             }
 827
 828             if (DOMAINDECOMP(cr))
 829             {
 830                 /* Repartition the domain decomposition */
 831                 dd_partition_system(fplog, mdlog, step, cr, bMasterState, nstglobalcomm, state_global,
 832                                     *top_global, ir, imdSession, pull_work, state, &f, mdAtoms, &top,
 833                                     fr, vsite, constr, nrnb, wcycle, do_verbose && !bPMETunePrinting);
 834                 shouldCheckNumberOfBondedInteractions = true;
 835                 upd.setNumAtoms(state->natoms);
 836             }
 837         }
 838
 839         // Allocate or re-size GPU halo exchange object, if necessary
 840         if (bNS && havePPDomainDecomposition(cr) && simulationWork.useGpuHaloExchange)
 841         {
 842             GMX_RELEASE_ASSERT(fr->deviceStreamManager != nullptr,
 843                                "GPU device manager has to be initialized to use GPU "
 844                                "version of halo exchange.");
 845             constructGpuHaloExchange(mdlog, *cr, *fr->deviceStreamManager, wcycle);
 846         }
 847
 848         if (MASTER(cr) && do_log)
 849         {
 850             gmx::EnergyOutput::printHeader(fplog, step,
 851                                            t); /* can we improve the information printed here? */
 852         }
 853
 854         if (ir->efep != efepNO)
 855         {
 856             update_mdatoms(mdatoms, state->lambda[efptMASS]);
 857         }
 858
 859         if (bExchanged)
 860         {
 861
 862             /* We need the kinetic energy at minus the half step for determining
 863              * the full step kinetic energy and possibly for T-coupling.*/
 864             /* This may not be quite working correctly yet . . . . */
 865             compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
 866                             makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
 867                             enerd, nullptr, nullptr, nullptr, nullptr, constr, &nullSignaller,
 868                             state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
 869                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
 870             checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global,
 871                                             &top, makeConstArrayRef(state->x), state->box,
 872                                             &shouldCheckNumberOfBondedInteractions);
 873         }
 874         clear_mat(force_vir);
 875
 876         checkpointHandler->decideIfCheckpointingThisStep(bNS, bFirstStep, bLastStep);
 877
 878         /* Determine the energy and pressure:
 879          * at nstcalcenergy steps and at energy output steps (set below).
 880          */
 881         if (EI_VV(ir->eI) && (!bInitStep))
 882         {
 883             bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
 884             bCalcVir      = bCalcEnerStep
 885                        || (ir->epc != epcNO
 886                            && (do_per_step(step, ir->nstpcouple) || do_per_step(step - 1, ir->nstpcouple)));
 887         }
 888         else
 889         {
 890             bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
 891             bCalcVir = bCalcEnerStep || (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
 892         }
 893         bCalcEner = bCalcEnerStep;
 894
 895         do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
 896
 897         if (do_ene || do_log || bDoReplEx)
 898         {
 899             bCalcVir  = TRUE;
 900             bCalcEner = TRUE;
 901         }
 902
 903         /* Do we need global communication ? */
 904         bGStat = (bCalcVir || bCalcEner || bStopCM || do_per_step(step, nstglobalcomm)
 905                   || (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step - 1, nstglobalcomm)));
 906
 907         force_flags = (GMX_FORCE_STATECHANGED | ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0)
 908                        | GMX_FORCE_ALLFORCES | (bCalcVir ? GMX_FORCE_VIRIAL : 0)
 909                        | (bCalcEner ? GMX_FORCE_ENERGY : 0) | (bDoFEP ? GMX_FORCE_DHDL : 0));
 910         if (fr->useMts && !do_per_step(step, ir->nstfout))
 911         {
 912             force_flags |= GMX_FORCE_DO_NOT_NEED_NORMAL_FORCE;
 913         }
 914
 915         if (shellfc)
 916         {
 917             /* Now is the time to relax the shells */
 918             relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir,
 919                                 imdSession, pull_work, bNS, force_flags, &top, constr, enerd,
 920                                 state->natoms, state->x.arrayRefWithPadding(),
 921                                 state->v.arrayRefWithPadding(), state->box, state->lambda,
 922                                 &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc,
 923                                 fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler);
 924         }
 925         else
 926         {
 927             /* The AWH history need to be saved _before_ doing force calculations where the AWH bias
 928                is updated (or the AWH update will be performed twice for one step when continuing).
 929                It would be best to call this update function from do_md_trajectory_writing but that
 930                would occur after do_force. One would have to divide the update_awh function into one
 931                function applying the AWH force and one doing the AWH bias update. The update AWH
 932                bias function could then be called after do_md_trajectory_writing (then containing
 933                update_awh_history). The checkpointing will in the future probably moved to the start
 934                of the md loop which will rid of this issue. */
 935             if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
 936             {
 937                 awh->updateHistory(state_global->awhHistory.get());
 938             }
 939
 940             /* The coordinates (x) are shifted (to get whole molecules)
 941              * in do_force.
 942              * This is parallellized as well, and does communication too.
 943              * Check comments in sim_util.c
 944              */
 945             do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession, pull_work, step,
 946                      nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist,
 947                      &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork,
 948                      vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
 949                      (bNS ? GMX_FORCE_NS : 0) | force_flags, ddBalanceRegionHandler);
 950         }
 951
 952         // VV integrators do not need the following velocity half step
 953         // if it is the first step after starting from a checkpoint.
 954         // That is, the half step is needed on all other steps, and
 955         // also the first step when starting from a .tpr file.
 956         if (EI_VV(ir->eI) && (!bFirstStep || startingBehavior == StartingBehavior::NewSimulation))
 957         /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
 958         {
 959             rvec* vbuf = nullptr;
 960
 961             wallcycle_start(wcycle, ewcUPDATE);
 962             if (ir->eI == eiVV && bInitStep)
 963             {
 964                 /* if using velocity verlet with full time step Ekin,
 965                  * take the first half step only to compute the
 966                  * virial for the first step. From there,
 967                  * revert back to the initial coordinates
 968                  * so that the input is actually the initial step.
 969                  */
 970                 snew(vbuf, state->natoms);
 971                 copy_rvecn(state->v.rvec_array(), vbuf, 0,
 972                            state->natoms); /* should make this better for parallelizing? */
 973             }
 974             else
 975             {
 976                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
 977                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ,
 978                                trotter_seq, ettTSEQ1);
 979             }
 980
 981             upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
 982                               M, etrtVELOCITY1, cr, constr != nullptr);
 983
 984             wallcycle_stop(wcycle, ewcUPDATE);
 985             constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir);
 986             wallcycle_start(wcycle, ewcUPDATE);
 987             /* if VV, compute the pressure and constraints */
 988             /* For VV2, we strictly only need this if using pressure
 989              * control, but we really would like to have accurate pressures
 990              * printed out.
 991              * Think about ways around this in the future?
 992              * For now, keep this choice in comments.
 993              */
 994             /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
 995             /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
 996             bPres = TRUE;
 997             bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
 998             if (bCalcEner && ir->eI == eiVVAK)
 999             {
1000                 bSumEkinhOld = TRUE;
1001             }
1002             /* for vv, the first half of the integration actually corresponds to the previous step.
1003                So we need information from the last step in the first half of the integration */
1004             if (bGStat || do_per_step(step - 1, nstglobalcomm))
1005             {
1006                 wallcycle_stop(wcycle, ewcUPDATE);
1007                 compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
1008                                 makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
1009                                 enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
1010                                 state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
1011                                 (bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
1012                                         | (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
1013                                         | (bPres ? CGLO_CONSTRAINT : 0) | (bStopCM ? CGLO_STOPCM : 0)
1014                                         | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
1015                                                                                  : 0)
1016                                         | CGLO_SCALEEKIN);
1017                 /* explanation of above:
1018                    a) We compute Ekin at the full time step
1019                    if 1) we are using the AveVel Ekin, and it's not the
1020                    initial step, or 2) if we are using AveEkin, but need the full
1021                    time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1022                    b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1023                    EkinAveVel because it's needed for the pressure */
1024                 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1025                                                 top_global, &top, makeConstArrayRef(state->x),
1026                                                 state->box, &shouldCheckNumberOfBondedInteractions);
1027                 if (bStopCM)
1028                 {
1029                     process_and_stopcm_grp(fplog, &vcm, *mdatoms, makeArrayRef(state->x),
1030                                            makeArrayRef(state->v));
1031                     inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
1032                 }
1033                 wallcycle_start(wcycle, ewcUPDATE);
1034             }
1035             /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1036             if (!bInitStep)
1037             {
1038                 if (bTrotter)
1039                 {
1040                     m_add(force_vir, shake_vir,
1041                           total_vir); /* we need the un-dispersion corrected total vir here */
1042                     trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ,
1043                                    trotter_seq, ettTSEQ2);
1044
1045                     /* TODO This is only needed when we're about to write
1046                      * a checkpoint, because we use it after the restart
1047                      * (in a kludge?). But what should we be doing if
1048                      * the startingBehavior is NewSimulation or bInitStep are true? */
1049                     if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1050                     {
1051                         copy_mat(shake_vir, state->svir_prev);
1052                         copy_mat(force_vir, state->fvir_prev);
1053                     }
1054                     if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
1055                     {
1056                         /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1057                         enerd->term[F_TEMP] =
1058                                 sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
1059                         enerd->term[F_EKIN] = trace(ekind->ekin);
1060                     }
1061                 }
1062                 else if (bExchanged)
1063                 {
1064                     wallcycle_stop(wcycle, ewcUPDATE);
1065                     /* We need the kinetic energy at minus the half step for determining
1066                      * the full step kinetic energy and possibly for T-coupling.*/
1067                     /* This may not be quite working correctly yet . . . . */
1068                     compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
1069                                     makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
1070                                     enerd, nullptr, nullptr, nullptr, nullptr, constr, &nullSignaller,
1071                                     state->box, nullptr, &bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
1072                     wallcycle_start(wcycle, ewcUPDATE);
1073                 }
1074             }
1075             /* if it's the initial step, we performed this first step just to get the constraint virial */
1076             if (ir->eI == eiVV && bInitStep)
1077             {
1078                 copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
1079                 sfree(vbuf);
1080             }
1081             wallcycle_stop(wcycle, ewcUPDATE);
1082         }
1083
1084         /* compute the conserved quantity */
1085         if (EI_VV(ir->eI))
1086         {
1087             saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1088             if (ir->eI == eiVV)
1089             {
1090                 last_ekin = enerd->term[F_EKIN];
1091             }
1092             if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1093             {
1094                 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1095             }
1096             /* sum up the foreign kinetic energy and dK/dl terms for vv.  currently done every step so that dhdl is correct in the .edr */
1097             if (ir->efep != efepNO)
1098             {
1099                 accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
1100             }
1101         }
1102
1103         /* ########  END FIRST UPDATE STEP  ############## */
1104         /* ########  If doing VV, we now have v(dt) ###### */
1105         if (bDoExpanded)
1106         {
1107             /* perform extended ensemble sampling in lambda - we don't
1108                actually move to the new state before outputting
1109                statistics, but if performing simulated tempering, we
1110                do update the velocities and the tau_t. */
1111
1112             lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state,
1113                                               state->dfhist, step, state->v.rvec_array(), mdatoms);
1114             /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1115             if (MASTER(cr))
1116             {
1117                 copy_df_history(state_global->dfhist, state->dfhist);
1118             }
1119         }
1120
1121         // Copy coordinate from the GPU for the output/checkpointing if the update is offloaded and
1122         // coordinates have not already been copied for i) search or ii) CPU force tasks.
1123         if (useGpuForUpdate && !bNS && !runScheduleWork->domainWork.haveCpuLocalForceWork
1124             && (do_per_step(step, ir->nstxout) || do_per_step(step, ir->nstxout_compressed)
1125                 || checkpointHandler->isCheckpointingStep()))
1126         {
1127             stateGpu->copyCoordinatesFromGpu(state->x, AtomLocality::Local);
1128             stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
1129         }
1130         // Copy velocities if needed for the output/checkpointing.
1131         // NOTE: Copy on the search steps is done at the beginning of the step.
1132         if (useGpuForUpdate && !bNS
1133             && (do_per_step(step, ir->nstvout) || checkpointHandler->isCheckpointingStep()))
1134         {
1135             stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
1136             stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
1137         }
1138         // Copy forces for the output if the forces were reduced on the GPU (not the case on virial steps)
1139         // and update is offloaded hence forces are kept on the GPU for update and have not been
1140         // already transferred in do_force().
1141         // TODO: There should be an improved, explicit mechanism that ensures this copy is only executed
1142         //       when the forces are ready on the GPU -- the same synchronizer should be used as the one
1143         //       prior to GPU update.
1144         // TODO: When the output flags will be included in step workload, this copy can be combined with the
1145         //       copy call in do_force(...).
1146         // NOTE: The forces should not be copied here if the vsites are present, since they were modified
1147         //       on host after the D2H copy in do_force(...).
1148         if (runScheduleWork->stepWork.useGpuFBufferOps && (simulationWork.useGpuUpdate && !vsite)
1149             && do_per_step(step, ir->nstfout))
1150         {
1151             stateGpu->copyForcesFromGpu(f.view().force(), AtomLocality::Local);
1152             stateGpu->waitForcesReadyOnHost(AtomLocality::Local);
1153         }
1154         /* Now we have the energies and forces corresponding to the
1155          * coordinates at time t. We must output all of this before
1156          * the update.
1157          */
1158         do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, state_global,
1159                                  observablesHistory, top_global, fr, outf, energyOutput, ekind,
1160                                  f.view().force(), checkpointHandler->isCheckpointingStep(),
1161                                  bRerunMD, bLastStep, mdrunOptions.writeConfout, bSumEkinhOld);
1162         /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1163         bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x.rvec_array(), t);
1164
1165         /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1166         if (startingBehavior != StartingBehavior::NewSimulation && bFirstStep
1167             && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1168         {
1169             copy_mat(state->svir_prev, shake_vir);
1170             copy_mat(state->fvir_prev, force_vir);
1171         }
1172
1173         stopHandler->setSignal();
1174         resetHandler->setSignal(walltime_accounting);
1175
1176         if (bGStat || !PAR(cr))
1177         {
1178             /* In parallel we only have to check for checkpointing in steps
1179              * where we do global communication,
1180              *  otherwise the other nodes don't know.
1181              */
1182             checkpointHandler->setSignal(walltime_accounting);
1183         }
1184
1185         /* #########   START SECOND UPDATE STEP ################# */
1186
1187         /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen
1188            controlled in preprocessing */
1189
1190         if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1191         {
1192             gmx_bool bIfRandomize;
1193             bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
1194             /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1195             if (constr && bIfRandomize)
1196             {
1197                 constrain_velocities(constr, do_log, do_ene, step, state, nullptr, false, nullptr);
1198             }
1199         }
1200         /* Box is changed in update() when we do pressure coupling,
1201          * but we should still use the old box for energy corrections and when
1202          * writing it to the energy file, so it matches the trajectory files for
1203          * the same timestep above. Make a copy in a separate array.
1204          */
1205         copy_mat(state->box, lastbox);
1206
1207         dvdl_constr = 0;
1208
1209         if (!useGpuForUpdate)
1210         {
1211             wallcycle_start(wcycle, ewcUPDATE);
1212         }
1213         /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1214         if (bTrotter)
1215         {
1216             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1217             /* We can only do Berendsen coupling after we have summed
1218              * the kinetic energy or virial. Since the happens
1219              * in global_state after update, we should only do it at
1220              * step % nstlist = 1 with bGStatEveryStep=FALSE.
1221              */
1222         }
1223         else
1224         {
1225             update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1226             update_pcouple_before_coordinates(fplog, step, ir, state, pressureCouplingMu, M, bInitStep);
1227         }
1228
1229         if (EI_VV(ir->eI))
1230         {
1231             /* velocity half-step update */
1232             upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
1233                               M, etrtVELOCITY2, cr, constr != nullptr);
1234         }
1235
1236         /* Above, initialize just copies ekinh into ekin,
1237          * it doesn't copy position (for VV),
1238          * and entire integrator for MD.
1239          */
1240
1241         if (ir->eI == eiVVAK)
1242         {
1243             cbuf.resize(state->x.size());
1244             std::copy(state->x.begin(), state->x.end(), cbuf.begin());
1245         }
1246
1247         /* With leap-frog type integrators we compute the kinetic energy
1248          * at a whole time step as the average of the half-time step kinetic
1249          * energies of two subsequent steps. Therefore we need to compute the
1250          * half step kinetic energy also if we need energies at the next step.
1251          */
1252         const bool needHalfStepKineticEnergy =
1253                 (!EI_VV(ir->eI) && (do_per_step(step + 1, nstglobalcomm) || step_rel + 1 == ir->nsteps));
1254
1255         // Parrinello-Rahman requires the pressure to be availible before the update to compute
1256         // the velocity scaling matrix. Hence, it runs one step after the nstpcouple step.
1257         const bool doParrinelloRahman = (ir->epc == epcPARRINELLORAHMAN
1258                                          && do_per_step(step + ir->nstpcouple - 1, ir->nstpcouple));
1259
1260         if (useGpuForUpdate)
1261         {
1262             if (bNS && (bFirstStep || DOMAINDECOMP(cr)))
1263             {
1264                 integrator->set(stateGpu->getCoordinates(), stateGpu->getVelocities(),
1265                                 stateGpu->getForces(), top.idef, *mdatoms, ekind->ngtc);
1266
1267                 // Copy data to the GPU after buffers might have being reinitialized
1268                 stateGpu->copyVelocitiesToGpu(state->v, AtomLocality::Local);
1269                 stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
1270             }
1271
1272             if (simulationWork.useGpuPme && !runScheduleWork->simulationWork.useGpuPmePpCommunication
1273                 && !thisRankHasDuty(cr, DUTY_PME))
1274             {
1275                 // The PME forces were recieved to the host, so have to be copied
1276                 stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::All);
1277             }
1278             else if (!runScheduleWork->stepWork.useGpuFBufferOps)
1279             {
1280                 // The buffer ops were not offloaded this step, so the forces are on the
1281                 // host and have to be copied
1282                 stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::Local);
1283             }
1284
1285             const bool doTemperatureScaling =
1286                     (ir->etc != etcNO && do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
1287
1288             // This applies Leap-Frog, LINCS and SETTLE in succession
1289             integrator->integrate(stateGpu->getForcesReadyOnDeviceEvent(
1290                                           AtomLocality::Local, runScheduleWork->stepWork.useGpuFBufferOps),
1291                                   ir->delta_t, true, bCalcVir, shake_vir, doTemperatureScaling,
1292                                   ekind->tcstat, doParrinelloRahman, ir->nstpcouple * ir->delta_t, M);
1293
1294             // Copy velocities D2H after update if:
1295             // - Globals are computed this step (includes the energy output steps).
1296             // - Temperature is needed for the next step.
1297             if (bGStat || needHalfStepKineticEnergy)
1298             {
1299                 stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
1300                 stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
1301             }
1302         }
1303         else
1304         {
1305             /* With multiple time stepping we need to do an additional normal
1306              * update step to obtain the virial, as the actual MTS integration
1307              * using an acceleration where the slow forces are multiplied by mtsFactor.
1308              * Using that acceleration would result in a virial with the slow
1309              * force contribution would be a factor mtsFactor too large.
1310              */
1311             if (fr->useMts && bCalcVir && constr != nullptr)
1312             {
1313                 upd.update_for_constraint_virial(*ir, *mdatoms, *state, f.view().forceWithPadding(), *ekind);
1314
1315                 constrain_coordinates(constr, do_log, do_ene, step, state,
1316                                       upd.xp()->arrayRefWithPadding(), &dvdl_constr, bCalcVir, shake_vir);
1317             }
1318
1319             ArrayRefWithPadding<const RVec> forceCombined =
1320                     (fr->useMts && step % ir->mtsLevels[1].stepFactor == 0)
1321                             ? f.view().forceMtsCombinedWithPadding()
1322                             : f.view().forceWithPadding();
1323             upd.update_coords(*ir, step, mdatoms, state, forceCombined, fcdata, ekind, M,
1324                               etrtPOSITION, cr, constr != nullptr);
1325
1326             wallcycle_stop(wcycle, ewcUPDATE);
1327
1328             constrain_coordinates(constr, do_log, do_ene, step, state, upd.xp()->arrayRefWithPadding(),
1329                                   &dvdl_constr, bCalcVir && !fr->useMts, shake_vir);
1330
1331             upd.update_sd_second_half(*ir, step, &dvdl_constr, mdatoms, state, cr, nrnb, wcycle,
1332                                       constr, do_log, do_ene);
1333             upd.finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
1334         }
1335
1336         if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1337         {
1338             updatePrevStepPullCom(pull_work, state);
1339         }
1340
1341         if (ir->eI == eiVVAK)
1342         {
1343             /* erase F_EKIN and F_TEMP here? */
1344             /* just compute the kinetic energy at the half step to perform a trotter step */
1345             compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
1346                             makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle, enerd,
1347                             force_vir, shake_vir, total_vir, pres, constr, &nullSignaller, lastbox,
1348                             nullptr, &bSumEkinhOld, (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE);
1349             wallcycle_start(wcycle, ewcUPDATE);
1350             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1351             /* now we know the scaling, we can compute the positions again */
1352             std::copy(cbuf.begin(), cbuf.end(), state->x.begin());
1353
1354             upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
1355                               M, etrtPOSITION, cr, constr != nullptr);
1356             wallcycle_stop(wcycle, ewcUPDATE);
1357
1358             /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1359             /* are the small terms in the shake_vir here due
1360              * to numerical errors, or are they important
1361              * physically? I'm thinking they are just errors, but not completely sure.
1362              * For now, will call without actually constraining, constr=NULL*/
1363             upd.finish_update(*ir, mdatoms, state, wcycle, false);
1364         }
1365         if (EI_VV(ir->eI))
1366         {
1367             /* this factor or 2 correction is necessary
1368                because half of the constraint force is removed
1369                in the vv step, so we have to double it.  See
1370                the Issue #1255.  It is not yet clear
1371                if the factor of 2 is exact, or just a very
1372                good approximation, and this will be
1373                investigated.  The next step is to see if this
1374                can be done adding a dhdl contribution from the
1375                rattle step, but this is somewhat more
1376                complicated with the current code. Will be
1377                investigated, hopefully for 4.6.3. However,
1378                this current solution is much better than
1379                having it completely wrong.
1380              */
1381             enerd->term[F_DVDL_CONSTR] += 2 * dvdl_constr;
1382         }
1383         else
1384         {
1385             enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1386         }
1387
1388         if (vsite != nullptr)
1389         {
1390             wallcycle_start(wcycle, ewcVSITECONSTR);
1391             vsite->construct(state->x, ir->delta_t, state->v, state->box);
1392             wallcycle_stop(wcycle, ewcVSITECONSTR);
1393         }
1394
1395         /* ############## IF NOT VV, Calculate globals HERE  ############ */
1396         /* With Leap-Frog we can skip compute_globals at
1397          * non-communication steps, but we need to calculate
1398          * the kinetic energy one step before communication.
1399          */
1400         {
1401             // Organize to do inter-simulation signalling on steps if
1402             // and when algorithms require it.
1403             const bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1404
1405             if (bGStat || needHalfStepKineticEnergy || doInterSimSignal)
1406             {
1407                 // Copy coordinates when needed to stop the CM motion.
1408                 if (useGpuForUpdate && !EI_VV(ir->eI) && bStopCM)
1409                 {
1410                     stateGpu->copyCoordinatesFromGpu(state->x, AtomLocality::Local);
1411                     stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
1412                 }
1413                 // Since we're already communicating at this step, we
1414                 // can propagate intra-simulation signals. Note that
1415                 // check_nstglobalcomm has the responsibility for
1416                 // choosing the value of nstglobalcomm that is one way
1417                 // bGStat becomes true, so we can't get into a
1418                 // situation where e.g. checkpointing can't be
1419                 // signalled.
1420                 bool                doIntraSimSignal = true;
1421                 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1422
1423                 compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
1424                                 makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm,
1425                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, constr,
1426                                 &signaller, lastbox, &totalNumberOfBondedInteractions, &bSumEkinhOld,
1427                                 (bGStat ? CGLO_GSTAT : 0) | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1428                                         | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1429                                         | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1430                                         | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0) | CGLO_CONSTRAINT
1431                                         | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
1432                                                                                  : 0));
1433                 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1434                                                 top_global, &top, makeConstArrayRef(state->x),
1435                                                 state->box, &shouldCheckNumberOfBondedInteractions);
1436                 if (!EI_VV(ir->eI) && bStopCM)
1437                 {
1438                     process_and_stopcm_grp(fplog, &vcm, *mdatoms, makeArrayRef(state->x),
1439                                            makeArrayRef(state->v));
1440                     inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
1441
1442                     // TODO: The special case of removing CM motion should be dealt more gracefully
1443                     if (useGpuForUpdate)
1444                     {
1445                         stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
1446                         // Here we block until the H2D copy completes because event sync with the
1447                         // force kernels that use the coordinates on the next steps is not implemented
1448                         // (not because of a race on state->x being modified on the CPU while H2D is in progress).
1449                         stateGpu->waitCoordinatesCopiedToDevice(AtomLocality::Local);
1450                         // If the COM removal changed the velocities on the CPU, this has to be accounted for.
1451                         if (vcm.mode != ecmNO)
1452                         {
1453                             stateGpu->copyVelocitiesToGpu(state->v, AtomLocality::Local);
1454                         }
1455                     }
1456                 }
1457             }
1458         }
1459
1460         /* #############  END CALC EKIN AND PRESSURE ################# */
1461
1462         /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1463            the virial that should probably be addressed eventually. state->veta has better properies,
1464            but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1465            generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
1466
1467         if (ir->efep != efepNO && !EI_VV(ir->eI))
1468         {
1469             /* Sum up the foreign energy and dK/dl terms for md and sd.
1470                Currently done every step so that dH/dl is correct in the .edr */
1471             accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
1472         }
1473
1474         update_pcouple_after_coordinates(fplog, step, ir, mdatoms, pres, force_vir, shake_vir,
1475                                          pressureCouplingMu, state, nrnb, upd.deform(), !useGpuForUpdate);
1476
1477         const bool doBerendsenPressureCoupling =
1478                 (inputrec->epc == epcBERENDSEN && do_per_step(step, inputrec->nstpcouple));
1479         const bool doCRescalePressureCoupling =
1480                 (inputrec->epc == epcCRESCALE && do_per_step(step, inputrec->nstpcouple));
1481         if (useGpuForUpdate
1482             && (doBerendsenPressureCoupling || doCRescalePressureCoupling || doParrinelloRahman))
1483         {
1484             integrator->scaleCoordinates(pressureCouplingMu);
1485             if (doCRescalePressureCoupling)
1486             {
1487                 matrix pressureCouplingInvMu;
1488                 gmx::invertBoxMatrix(pressureCouplingMu, pressureCouplingInvMu);
1489                 integrator->scaleVelocities(pressureCouplingInvMu);
1490             }
1491             integrator->setPbc(PbcType::Xyz, state->box);
1492         }
1493
1494         /* ################# END UPDATE STEP 2 ################# */
1495         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
1496
1497         /* The coordinates (x) were unshifted in update */
1498         if (!bGStat)
1499         {
1500             /* We will not sum ekinh_old,
1501              * so signal that we still have to do it.
1502              */
1503             bSumEkinhOld = TRUE;
1504         }
1505
1506         if (bCalcEner)
1507         {
1508             /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
1509
1510             /* use the directly determined last velocity, not actually the averaged half steps */
1511             if (bTrotter && ir->eI == eiVV)
1512             {
1513                 enerd->term[F_EKIN] = last_ekin;
1514             }
1515             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1516
1517             if (integratorHasConservedEnergyQuantity(ir))
1518             {
1519                 if (EI_VV(ir->eI))
1520                 {
1521                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1522                 }
1523                 else
1524                 {
1525                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1526                 }
1527             }
1528             /* #########  END PREPARING EDR OUTPUT  ###########  */
1529         }
1530
1531         /* Output stuff */
1532         if (MASTER(cr))
1533         {
1534             if (fplog && do_log && bDoExpanded)
1535             {
1536                 /* only needed if doing expanded ensemble */
1537                 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals,
1538                                           ir->bSimTemp ? ir->simtempvals : nullptr,
1539                                           state_global->dfhist, state->fep_state, ir->nstlog, step);
1540             }
1541             if (bCalcEner)
1542             {
1543                 energyOutput.addDataAtEnergyStep(
1544                         bDoDHDL, bCalcEnerStep, t, mdatoms->tmass, enerd, ir->fepvals,
1545                         ir->expandedvals, lastbox,
1546                         PTCouplingArrays{ state->boxv, state->nosehoover_xi, state->nosehoover_vxi,
1547                                           state->nhpres_xi, state->nhpres_vxi },
1548                         state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr);
1549             }
1550             else
1551             {
1552                 energyOutput.recordNonEnergyStep();
1553             }
1554
1555             gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
1556             gmx_bool do_or = do_per_step(step, ir->nstorireout);
1557
1558             if (doSimulatedAnnealing)
1559             {
1560                 gmx::EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups,
1561                                                               &(ir->opts));
1562             }
1563             if (do_log || do_ene || do_dr || do_or)
1564             {
1565                 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
1566                                                    do_log ? fplog : nullptr, step, t,
1567                                                    fr->fcdata.get(), awh.get());
1568             }
1569             if (do_log && ir->bDoAwh && awh->hasFepLambdaDimension())
1570             {
1571                 const bool isInitialOutput = false;
1572                 printLambdaStateToLog(fplog, state->lambda, isInitialOutput);
1573             }
1574
1575             if (ir->bPull)
1576             {
1577                 pull_print_output(pull_work, step, t);
1578             }
1579
1580             if (do_per_step(step, ir->nstlog))
1581             {
1582                 if (fflush(fplog) != 0)
1583                 {
1584                     gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1585                 }
1586             }
1587         }
1588         if (bDoExpanded)
1589         {
1590             /* Have to do this part _after_ outputting the logfile and the edr file */
1591             /* Gets written into the state at the beginning of next loop*/
1592             state->fep_state = lamnew;
1593         }
1594         else if (ir->bDoAwh && awh->needForeignEnergyDifferences(step))
1595         {
1596             state->fep_state = awh->fepLambdaState();
1597         }
1598         /* Print the remaining wall clock time for the run */
1599         if (isMasterSimMasterRank(ms, MASTER(cr)) && (do_verbose || gmx_got_usr_signal()) && !bPMETunePrinting)
1600         {
1601             if (shellfc)
1602             {
1603                 fprintf(stderr, "\n");
1604             }
1605             print_time(stderr, walltime_accounting, step, ir, cr);
1606         }
1607
1608         /* Ion/water position swapping.
1609          * Not done in last step since trajectory writing happens before this call
1610          * in the MD loop and exchanges would be lost anyway. */
1611         bNeedRepartition = FALSE;
1612         if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep && do_per_step(step, ir->swap->nstswap))
1613         {
1614             bNeedRepartition =
1615                     do_swapcoords(cr, step, t, ir, swap, wcycle, as_rvec_array(state->x.data()),
1616                                   state->box, MASTER(cr) && mdrunOptions.verbose, bRerunMD);
1617
1618             if (bNeedRepartition && DOMAINDECOMP(cr))
1619             {
1620                 dd_collect_state(cr->dd, state, state_global);
1621             }
1622         }
1623
1624         /* Replica exchange */
1625         bExchanged = FALSE;
1626         if (bDoReplEx)
1627         {
1628             bExchanged = replica_exchange(fplog, cr, ms, repl_ex, state_global, enerd, state, step, t);
1629         }
1630
1631         if ((bExchanged || bNeedRepartition) && DOMAINDECOMP(cr))
1632         {
1633             dd_partition_system(fplog, mdlog, step, cr, TRUE, 1, state_global, *top_global, ir,
1634                                 imdSession, pull_work, state, &f, mdAtoms, &top, fr, vsite, constr,
1635                                 nrnb, wcycle, FALSE);
1636             shouldCheckNumberOfBondedInteractions = true;
1637             upd.setNumAtoms(state->natoms);
1638         }
1639
1640         bFirstStep = FALSE;
1641         bInitStep  = FALSE;
1642
1643         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1644         /* With all integrators, except VV, we need to retain the pressure
1645          * at the current step for coupling at the next step.
1646          */
1647         if ((state->flags & (1U << estPRES_PREV))
1648             && (bGStatEveryStep || (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1649         {
1650             /* Store the pressure in t_state for pressure coupling
1651              * at the next MD step.
1652              */
1653             copy_mat(pres, state->pres_prev);
1654         }
1655
1656         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
1657
1658         if ((membed != nullptr) && (!bLastStep))
1659         {
1660             rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1661         }
1662
1663         cycles = wallcycle_stop(wcycle, ewcSTEP);
1664         if (DOMAINDECOMP(cr) && wcycle)
1665         {
1666             dd_cycles_add(cr->dd, cycles, ddCyclStep);
1667         }
1668
1669         /* increase the MD step number */
1670         step++;
1671         step_rel++;
1672
1673 #if GMX_FAHCORE
1674         if (MASTER(cr))
1675         {
1676             fcReportProgress(ir->nsteps + ir->init_step, step);
1677         }
1678 #endif
1679
1680         resetHandler->resetCounters(step, step_rel, mdlog, fplog, cr, fr->nbv.get(), nrnb,
1681                                     fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
1682
1683         /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1684         imdSession->updateEnergyRecordAndSendPositionsAndEnergies(bInteractiveMDstep, step, bCalcEner);
1685     }
1686     /* End of main MD loop */
1687
1688     /* Closing TNG files can include compressing data. Therefore it is good to do that
1689      * before stopping the time measurements. */
1690     mdoutf_tng_close(outf);
1691
1692     /* Stop measuring walltime */
1693     walltime_accounting_end_time(walltime_accounting);
1694
1695     if (!thisRankHasDuty(cr, DUTY_PME))
1696     {
1697         /* Tell the PME only node to finish */
1698         gmx_pme_send_finish(cr);
1699     }
1700
1701     if (MASTER(cr))
1702     {
1703         if (ir->nstcalcenergy > 0)
1704         {
1705             energyOutput.printEnergyConservation(fplog, ir->simulation_part, EI_MD(ir->eI));
1706
1707             gmx::EnergyOutput::printAnnealingTemperatures(fplog, groups, &(ir->opts));
1708             energyOutput.printAverages(fplog, groups);
1709         }
1710     }
1711     done_mdoutf(outf);
1712
1713     if (bPMETune)
1714     {
1715         pme_loadbal_done(pme_loadbal, fplog, mdlog, fr->nbv->useGpu());
1716     }
1717
1718     done_shellfc(fplog, shellfc, step_rel);
1719
1720     if (useReplicaExchange && MASTER(cr))
1721     {
1722         print_replica_exchange_statistics(fplog, repl_ex);
1723     }
1724
1725     walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1726
1727     global_stat_destroy(gstat);
1728 }