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                         || (!haveFreeEnergyType(*ir, efptBONDED) && !haveFreeEnergyType(*ir, efptMASS)),
 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         wallcycle_start(wcycle, ewcUPDATE);
1210         /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1211         if (bTrotter)
1212         {
1213             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1214             /* We can only do Berendsen coupling after we have summed
1215              * the kinetic energy or virial. Since the happens
1216              * in global_state after update, we should only do it at
1217              * step % nstlist = 1 with bGStatEveryStep=FALSE.
1218              */
1219         }
1220         else
1221         {
1222             update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1223             update_pcouple_before_coordinates(fplog, step, ir, state, pressureCouplingMu, M, bInitStep);
1224         }
1225
1226         if (EI_VV(ir->eI))
1227         {
1228             /* velocity half-step update */
1229             upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
1230                               M, etrtVELOCITY2, cr, constr != nullptr);
1231         }
1232
1233         /* Above, initialize just copies ekinh into ekin,
1234          * it doesn't copy position (for VV),
1235          * and entire integrator for MD.
1236          */
1237
1238         if (ir->eI == eiVVAK)
1239         {
1240             cbuf.resize(state->x.size());
1241             std::copy(state->x.begin(), state->x.end(), cbuf.begin());
1242         }
1243
1244         /* With leap-frog type integrators we compute the kinetic energy
1245          * at a whole time step as the average of the half-time step kinetic
1246          * energies of two subsequent steps. Therefore we need to compute the
1247          * half step kinetic energy also if we need energies at the next step.
1248          */
1249         const bool needHalfStepKineticEnergy =
1250                 (!EI_VV(ir->eI) && (do_per_step(step + 1, nstglobalcomm) || step_rel + 1 == ir->nsteps));
1251
1252         // Parrinello-Rahman requires the pressure to be availible before the update to compute
1253         // the velocity scaling matrix. Hence, it runs one step after the nstpcouple step.
1254         const bool doParrinelloRahman = (ir->epc == epcPARRINELLORAHMAN
1255                                          && do_per_step(step + ir->nstpcouple - 1, ir->nstpcouple));
1256
1257         if (useGpuForUpdate)
1258         {
1259             wallcycle_stop(wcycle, ewcUPDATE);
1260
1261             if (bNS && (bFirstStep || DOMAINDECOMP(cr)))
1262             {
1263                 integrator->set(stateGpu->getCoordinates(), stateGpu->getVelocities(),
1264                                 stateGpu->getForces(), top.idef, *mdatoms, ekind->ngtc);
1265
1266                 // Copy data to the GPU after buffers might have being reinitialized
1267                 stateGpu->copyVelocitiesToGpu(state->v, AtomLocality::Local);
1268                 stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
1269             }
1270
1271             if (simulationWork.useGpuPme && !runScheduleWork->simulationWork.useGpuPmePpCommunication
1272                 && !thisRankHasDuty(cr, DUTY_PME))
1273             {
1274                 // The PME forces were recieved to the host, so have to be copied
1275                 stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::All);
1276             }
1277             else if (!runScheduleWork->stepWork.useGpuFBufferOps)
1278             {
1279                 // The buffer ops were not offloaded this step, so the forces are on the
1280                 // host and have to be copied
1281                 stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::Local);
1282             }
1283
1284             const bool doTemperatureScaling =
1285                     (ir->etc != etcNO && do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
1286
1287             // This applies Leap-Frog, LINCS and SETTLE in succession
1288             integrator->integrate(stateGpu->getForcesReadyOnDeviceEvent(
1289                                           AtomLocality::Local, runScheduleWork->stepWork.useGpuFBufferOps),
1290                                   ir->delta_t, true, bCalcVir, shake_vir, doTemperatureScaling,
1291                                   ekind->tcstat, doParrinelloRahman, ir->nstpcouple * ir->delta_t, M);
1292
1293             // Copy velocities D2H after update if:
1294             // - Globals are computed this step (includes the energy output steps).
1295             // - Temperature is needed for the next step.
1296             if (bGStat || needHalfStepKineticEnergy)
1297             {
1298                 stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
1299                 stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
1300             }
1301         }
1302         else
1303         {
1304             /* With multiple time stepping we need to do an additional normal
1305              * update step to obtain the virial, as the actual MTS integration
1306              * using an acceleration where the slow forces are multiplied by mtsFactor.
1307              * Using that acceleration would result in a virial with the slow
1308              * force contribution would be a factor mtsFactor too large.
1309              */
1310             if (fr->useMts && bCalcVir && constr != nullptr)
1311             {
1312                 upd.update_for_constraint_virial(*ir, *mdatoms, *state, f.view().forceWithPadding(), *ekind);
1313
1314                 constrain_coordinates(constr, do_log, do_ene, step, state,
1315                                       upd.xp()->arrayRefWithPadding(), &dvdl_constr, bCalcVir, shake_vir);
1316             }
1317
1318             ArrayRefWithPadding<const RVec> forceCombined =
1319                     (fr->useMts && step % ir->mtsLevels[1].stepFactor == 0)
1320                             ? f.view().forceMtsCombinedWithPadding()
1321                             : f.view().forceWithPadding();
1322             upd.update_coords(*ir, step, mdatoms, state, forceCombined, fcdata, ekind, M,
1323                               etrtPOSITION, cr, constr != nullptr);
1324
1325             wallcycle_stop(wcycle, ewcUPDATE);
1326
1327             constrain_coordinates(constr, do_log, do_ene, step, state, upd.xp()->arrayRefWithPadding(),
1328                                   &dvdl_constr, bCalcVir && !fr->useMts, shake_vir);
1329
1330             upd.update_sd_second_half(*ir, step, &dvdl_constr, mdatoms, state, cr, nrnb, wcycle,
1331                                       constr, do_log, do_ene);
1332             upd.finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
1333         }
1334
1335         if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1336         {
1337             updatePrevStepPullCom(pull_work, state);
1338         }
1339
1340         if (ir->eI == eiVVAK)
1341         {
1342             /* erase F_EKIN and F_TEMP here? */
1343             /* just compute the kinetic energy at the half step to perform a trotter step */
1344             compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
1345                             makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle, enerd,
1346                             force_vir, shake_vir, total_vir, pres, constr, &nullSignaller, lastbox,
1347                             nullptr, &bSumEkinhOld, (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE);
1348             wallcycle_start(wcycle, ewcUPDATE);
1349             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1350             /* now we know the scaling, we can compute the positions again */
1351             std::copy(cbuf.begin(), cbuf.end(), state->x.begin());
1352
1353             upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
1354                               M, etrtPOSITION, cr, constr != nullptr);
1355             wallcycle_stop(wcycle, ewcUPDATE);
1356
1357             /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1358             /* are the small terms in the shake_vir here due
1359              * to numerical errors, or are they important
1360              * physically? I'm thinking they are just errors, but not completely sure.
1361              * For now, will call without actually constraining, constr=NULL*/
1362             upd.finish_update(*ir, mdatoms, state, wcycle, false);
1363         }
1364         if (EI_VV(ir->eI))
1365         {
1366             /* this factor or 2 correction is necessary
1367                because half of the constraint force is removed
1368                in the vv step, so we have to double it.  See
1369                the Issue #1255.  It is not yet clear
1370                if the factor of 2 is exact, or just a very
1371                good approximation, and this will be
1372                investigated.  The next step is to see if this
1373                can be done adding a dhdl contribution from the
1374                rattle step, but this is somewhat more
1375                complicated with the current code. Will be
1376                investigated, hopefully for 4.6.3. However,
1377                this current solution is much better than
1378                having it completely wrong.
1379              */
1380             enerd->term[F_DVDL_CONSTR] += 2 * dvdl_constr;
1381         }
1382         else
1383         {
1384             enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1385         }
1386
1387         if (vsite != nullptr)
1388         {
1389             wallcycle_start(wcycle, ewcVSITECONSTR);
1390             vsite->construct(state->x, ir->delta_t, state->v, state->box);
1391             wallcycle_stop(wcycle, ewcVSITECONSTR);
1392         }
1393
1394         /* ############## IF NOT VV, Calculate globals HERE  ############ */
1395         /* With Leap-Frog we can skip compute_globals at
1396          * non-communication steps, but we need to calculate
1397          * the kinetic energy one step before communication.
1398          */
1399         {
1400             // Organize to do inter-simulation signalling on steps if
1401             // and when algorithms require it.
1402             const bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1403
1404             if (bGStat || needHalfStepKineticEnergy || doInterSimSignal)
1405             {
1406                 // Copy coordinates when needed to stop the CM motion.
1407                 if (useGpuForUpdate && !EI_VV(ir->eI) && bStopCM)
1408                 {
1409                     stateGpu->copyCoordinatesFromGpu(state->x, AtomLocality::Local);
1410                     stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
1411                 }
1412                 // Since we're already communicating at this step, we
1413                 // can propagate intra-simulation signals. Note that
1414                 // check_nstglobalcomm has the responsibility for
1415                 // choosing the value of nstglobalcomm that is one way
1416                 // bGStat becomes true, so we can't get into a
1417                 // situation where e.g. checkpointing can't be
1418                 // signalled.
1419                 bool                doIntraSimSignal = true;
1420                 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1421
1422                 compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
1423                                 makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm,
1424                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, constr,
1425                                 &signaller, lastbox, &totalNumberOfBondedInteractions, &bSumEkinhOld,
1426                                 (bGStat ? CGLO_GSTAT : 0) | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1427                                         | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1428                                         | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1429                                         | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0) | CGLO_CONSTRAINT
1430                                         | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
1431                                                                                  : 0));
1432                 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1433                                                 top_global, &top, makeConstArrayRef(state->x),
1434                                                 state->box, &shouldCheckNumberOfBondedInteractions);
1435                 if (!EI_VV(ir->eI) && bStopCM)
1436                 {
1437                     process_and_stopcm_grp(fplog, &vcm, *mdatoms, makeArrayRef(state->x),
1438                                            makeArrayRef(state->v));
1439                     inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
1440
1441                     // TODO: The special case of removing CM motion should be dealt more gracefully
1442                     if (useGpuForUpdate)
1443                     {
1444                         stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
1445                         // Here we block until the H2D copy completes because event sync with the
1446                         // force kernels that use the coordinates on the next steps is not implemented
1447                         // (not because of a race on state->x being modified on the CPU while H2D is in progress).
1448                         stateGpu->waitCoordinatesCopiedToDevice(AtomLocality::Local);
1449                         // If the COM removal changed the velocities on the CPU, this has to be accounted for.
1450                         if (vcm.mode != ecmNO)
1451                         {
1452                             stateGpu->copyVelocitiesToGpu(state->v, AtomLocality::Local);
1453                         }
1454                     }
1455                 }
1456             }
1457         }
1458
1459         /* #############  END CALC EKIN AND PRESSURE ################# */
1460
1461         /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1462            the virial that should probably be addressed eventually. state->veta has better properies,
1463            but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1464            generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
1465
1466         if (ir->efep != efepNO && !EI_VV(ir->eI))
1467         {
1468             /* Sum up the foreign energy and dK/dl terms for md and sd.
1469                Currently done every step so that dH/dl is correct in the .edr */
1470             accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
1471         }
1472
1473         update_pcouple_after_coordinates(fplog, step, ir, mdatoms, pres, force_vir, shake_vir,
1474                                          pressureCouplingMu, state, nrnb, upd.deform(), !useGpuForUpdate);
1475
1476         const bool doBerendsenPressureCoupling =
1477                 (inputrec->epc == epcBERENDSEN && do_per_step(step, inputrec->nstpcouple));
1478         const bool doCRescalePressureCoupling =
1479                 (inputrec->epc == epcCRESCALE && do_per_step(step, inputrec->nstpcouple));
1480         if (useGpuForUpdate
1481             && (doBerendsenPressureCoupling || doCRescalePressureCoupling || doParrinelloRahman))
1482         {
1483             integrator->scaleCoordinates(pressureCouplingMu);
1484             if (doCRescalePressureCoupling)
1485             {
1486                 matrix pressureCouplingInvMu;
1487                 gmx::invertBoxMatrix(pressureCouplingMu, pressureCouplingInvMu);
1488                 integrator->scaleVelocities(pressureCouplingInvMu);
1489             }
1490             integrator->setPbc(PbcType::Xyz, state->box);
1491         }
1492
1493         /* ################# END UPDATE STEP 2 ################# */
1494         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
1495
1496         /* The coordinates (x) were unshifted in update */
1497         if (!bGStat)
1498         {
1499             /* We will not sum ekinh_old,
1500              * so signal that we still have to do it.
1501              */
1502             bSumEkinhOld = TRUE;
1503         }
1504
1505         if (bCalcEner)
1506         {
1507             /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
1508
1509             /* use the directly determined last velocity, not actually the averaged half steps */
1510             if (bTrotter && ir->eI == eiVV)
1511             {
1512                 enerd->term[F_EKIN] = last_ekin;
1513             }
1514             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1515
1516             if (integratorHasConservedEnergyQuantity(ir))
1517             {
1518                 if (EI_VV(ir->eI))
1519                 {
1520                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1521                 }
1522                 else
1523                 {
1524                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1525                 }
1526             }
1527             /* #########  END PREPARING EDR OUTPUT  ###########  */
1528         }
1529
1530         /* Output stuff */
1531         if (MASTER(cr))
1532         {
1533             if (fplog && do_log && bDoExpanded)
1534             {
1535                 /* only needed if doing expanded ensemble */
1536                 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals,
1537                                           ir->bSimTemp ? ir->simtempvals : nullptr,
1538                                           state_global->dfhist, state->fep_state, ir->nstlog, step);
1539             }
1540             if (bCalcEner)
1541             {
1542                 energyOutput.addDataAtEnergyStep(
1543                         bDoDHDL, bCalcEnerStep, t, mdatoms->tmass, enerd, ir->fepvals,
1544                         ir->expandedvals, lastbox,
1545                         PTCouplingArrays{ state->boxv, state->nosehoover_xi, state->nosehoover_vxi,
1546                                           state->nhpres_xi, state->nhpres_vxi },
1547                         state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr);
1548             }
1549             else
1550             {
1551                 energyOutput.recordNonEnergyStep();
1552             }
1553
1554             gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
1555             gmx_bool do_or = do_per_step(step, ir->nstorireout);
1556
1557             if (doSimulatedAnnealing)
1558             {
1559                 gmx::EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups,
1560                                                               &(ir->opts));
1561             }
1562             if (do_log || do_ene || do_dr || do_or)
1563             {
1564                 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
1565                                                    do_log ? fplog : nullptr, step, t,
1566                                                    fr->fcdata.get(), awh.get());
1567             }
1568             if (do_log && ir->bDoAwh && awh->hasFepLambdaDimension())
1569             {
1570                 const bool isInitialOutput = false;
1571                 printLambdaStateToLog(fplog, state->lambda, isInitialOutput);
1572             }
1573
1574             if (ir->bPull)
1575             {
1576                 pull_print_output(pull_work, step, t);
1577             }
1578
1579             if (do_per_step(step, ir->nstlog))
1580             {
1581                 if (fflush(fplog) != 0)
1582                 {
1583                     gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1584                 }
1585             }
1586         }
1587         if (bDoExpanded)
1588         {
1589             /* Have to do this part _after_ outputting the logfile and the edr file */
1590             /* Gets written into the state at the beginning of next loop*/
1591             state->fep_state = lamnew;
1592         }
1593         else if (ir->bDoAwh && awh->needForeignEnergyDifferences(step))
1594         {
1595             state->fep_state = awh->fepLambdaState();
1596         }
1597         /* Print the remaining wall clock time for the run */
1598         if (isMasterSimMasterRank(ms, MASTER(cr)) && (do_verbose || gmx_got_usr_signal()) && !bPMETunePrinting)
1599         {
1600             if (shellfc)
1601             {
1602                 fprintf(stderr, "\n");
1603             }
1604             print_time(stderr, walltime_accounting, step, ir, cr);
1605         }
1606
1607         /* Ion/water position swapping.
1608          * Not done in last step since trajectory writing happens before this call
1609          * in the MD loop and exchanges would be lost anyway. */
1610         bNeedRepartition = FALSE;
1611         if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep && do_per_step(step, ir->swap->nstswap))
1612         {
1613             bNeedRepartition =
1614                     do_swapcoords(cr, step, t, ir, swap, wcycle, as_rvec_array(state->x.data()),
1615                                   state->box, MASTER(cr) && mdrunOptions.verbose, bRerunMD);
1616
1617             if (bNeedRepartition && DOMAINDECOMP(cr))
1618             {
1619                 dd_collect_state(cr->dd, state, state_global);
1620             }
1621         }
1622
1623         /* Replica exchange */
1624         bExchanged = FALSE;
1625         if (bDoReplEx)
1626         {
1627             bExchanged = replica_exchange(fplog, cr, ms, repl_ex, state_global, enerd, state, step, t);
1628         }
1629
1630         if ((bExchanged || bNeedRepartition) && DOMAINDECOMP(cr))
1631         {
1632             dd_partition_system(fplog, mdlog, step, cr, TRUE, 1, state_global, *top_global, ir,
1633                                 imdSession, pull_work, state, &f, mdAtoms, &top, fr, vsite, constr,
1634                                 nrnb, wcycle, FALSE);
1635             shouldCheckNumberOfBondedInteractions = true;
1636             upd.setNumAtoms(state->natoms);
1637         }
1638
1639         bFirstStep = FALSE;
1640         bInitStep  = FALSE;
1641
1642         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1643         /* With all integrators, except VV, we need to retain the pressure
1644          * at the current step for coupling at the next step.
1645          */
1646         if ((state->flags & (1U << estPRES_PREV))
1647             && (bGStatEveryStep || (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1648         {
1649             /* Store the pressure in t_state for pressure coupling
1650              * at the next MD step.
1651              */
1652             copy_mat(pres, state->pres_prev);
1653         }
1654
1655         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
1656
1657         if ((membed != nullptr) && (!bLastStep))
1658         {
1659             rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1660         }
1661
1662         cycles = wallcycle_stop(wcycle, ewcSTEP);
1663         if (DOMAINDECOMP(cr) && wcycle)
1664         {
1665             dd_cycles_add(cr->dd, cycles, ddCyclStep);
1666         }
1667
1668         /* increase the MD step number */
1669         step++;
1670         step_rel++;
1671
1672 #if GMX_FAHCORE
1673         if (MASTER(cr))
1674         {
1675             fcReportProgress(ir->nsteps + ir->init_step, step);
1676         }
1677 #endif
1678
1679         resetHandler->resetCounters(step, step_rel, mdlog, fplog, cr, fr->nbv.get(), nrnb,
1680                                     fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
1681
1682         /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1683         imdSession->updateEnergyRecordAndSendPositionsAndEnergies(bInteractiveMDstep, step, bCalcEner);
1684     }
1685     /* End of main MD loop */
1686
1687     /* Closing TNG files can include compressing data. Therefore it is good to do that
1688      * before stopping the time measurements. */
1689     mdoutf_tng_close(outf);
1690
1691     /* Stop measuring walltime */
1692     walltime_accounting_end_time(walltime_accounting);
1693
1694     if (!thisRankHasDuty(cr, DUTY_PME))
1695     {
1696         /* Tell the PME only node to finish */
1697         gmx_pme_send_finish(cr);
1698     }
1699
1700     if (MASTER(cr))
1701     {
1702         if (ir->nstcalcenergy > 0)
1703         {
1704             energyOutput.printEnergyConservation(fplog, ir->simulation_part, EI_MD(ir->eI));
1705
1706             gmx::EnergyOutput::printAnnealingTemperatures(fplog, groups, &(ir->opts));
1707             energyOutput.printAverages(fplog, groups);
1708         }
1709     }
1710     done_mdoutf(outf);
1711
1712     if (bPMETune)
1713     {
1714         pme_loadbal_done(pme_loadbal, fplog, mdlog, fr->nbv->useGpu());
1715     }
1716
1717     done_shellfc(fplog, shellfc, step_rel);
1718
1719     if (useReplicaExchange && MASTER(cr))
1720     {
1721         print_replica_exchange_statistics(fplog, repl_ex);
1722     }
1723
1724     walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1725
1726     global_stat_destroy(gstat);
1727 }