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