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
 334     StatePropagatorDataGpu *stateGpu = fr->stateGpu;
 335
 336     if (useGpuForUpdate)
 337     {
 338         GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr), "Domain decomposition is not supported with the GPU update.\n");
 339         GMX_RELEASE_ASSERT(useGpuForPme || (useGpuForNonbonded && simulationWork.useGpuBufferOps),
 340                            "Either PME or short-ranged non-bonded interaction tasks must run on the GPU to use GPU update.\n");
 341         GMX_RELEASE_ASSERT(ir->eI == eiMD, "Only the md integrator is supported with the GPU update.\n");
 342         GMX_RELEASE_ASSERT(ir->etc != etcNOSEHOOVER, "Nose-Hoover temperature coupling is not supported with the GPU update.\n");
 343         GMX_RELEASE_ASSERT(ir->epc == epcNO || ir->epc == epcPARRINELLORAHMAN, "Only Parrinello-Rahman pressure control is supported with the GPU update.\n");
 344         GMX_RELEASE_ASSERT(!mdatoms->haveVsites, "Virtual sites are not supported with the GPU update.\n");
 345         GMX_RELEASE_ASSERT(ed == nullptr, "Essential dynamics is not supported with the GPU update.\n");
 346         GMX_RELEASE_ASSERT(!ir->bPull && !ir->pull, "Pulling is not supported with the GPU update.\n");
 347         GMX_RELEASE_ASSERT(fcd->orires.nr == 0, "Orientation restraints are not supported with the GPU update.\n");
 348         GMX_RELEASE_ASSERT(fcd->disres.npair  == 0, "Distance restraints are not supported with the GPU update.\n");
 349         GMX_RELEASE_ASSERT(ir->efep == efepNO, "Free energy perturbations are not supported with the GPU update.");
 350
 351         if (constr != nullptr && constr->numConstraintsTotal() > 0)
 352         {
 353             GMX_LOG(mdlog.info).asParagraph().
 354                 appendText("Updating coordinates and applying constraints on the GPU.");
 355         }
 356         else
 357         {
 358             GMX_LOG(mdlog.info).asParagraph().
 359                 appendText("Updating coordinates on the GPU.");
 360         }
 361         integrator = std::make_unique<UpdateConstrainCuda>(*ir, *top_global, stateGpu->getUpdateStream(), stateGpu->xUpdatedOnDevice());
 362     }
 363
 364     if (useGpuForPme || (useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
 365     {
 366         changePinningPolicy(&state->x, PinningPolicy::PinnedIfSupported);
 367     }
 368     if ((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
 369     {
 370         changePinningPolicy(&f, PinningPolicy::PinnedIfSupported);
 371     }
 372     if (useGpuForUpdate)
 373     {
 374         changePinningPolicy(&state->v, PinningPolicy::PinnedIfSupported);
 375     }
 376 /*****************************************************************************************/
 377
 378     // NOTE: The global state is no longer used at this point.
 379     // But state_global is still used as temporary storage space for writing
 380     // the global state to file and potentially for replica exchange.
 381     // (Global topology should persist.)
 382
 383     update_mdatoms(mdatoms, state->lambda[efptMASS]);
 384
 385     if (ir->bExpanded)
 386     {
 387         /* Check nstexpanded here, because the grompp check was broken */
 388         if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
 389         {
 390             gmx_fatal(FARGS, "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
 391         }
 392         init_expanded_ensemble(startingBehavior != StartingBehavior::NewSimulation,
 393                                ir, state->dfhist);
 394     }
 395
 396     if (MASTER(cr))
 397     {
 398         EnergyElement::initializeEnergyHistory(
 399                 startingBehavior, observablesHistory, &energyOutput);
 400     }
 401
 402     preparePrevStepPullCom(ir, pull_work, mdatoms, state, state_global, cr,
 403                            startingBehavior != StartingBehavior::NewSimulation);
 404
 405     // TODO: Remove this by converting AWH into a ForceProvider
 406     auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingBehavior != StartingBehavior::NewSimulation,
 407                                 shellfc != nullptr,
 408                                 opt2fn("-awh", nfile, fnm), pull_work);
 409
 410     const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
 411     if (useReplicaExchange && MASTER(cr))
 412     {
 413         repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
 414                                         replExParams);
 415     }
 416     /* PME tuning is only supported in the Verlet scheme, with PME for
 417      * Coulomb. It is not supported with only LJ PME. */
 418     bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) &&
 419                 !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
 420
 421     pme_load_balancing_t *pme_loadbal      = nullptr;
 422     if (bPMETune)
 423     {
 424         pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box,
 425                          *fr->ic, *fr->nbv, fr->pmedata, fr->nbv->useGpu(),
 426                          &bPMETunePrinting);
 427     }
 428
 429     if (!ir->bContinuation)
 430     {
 431         if (state->flags & (1U << estV))
 432         {
 433             auto v = makeArrayRef(state->v);
 434             /* Set the velocities of vsites, shells and frozen atoms to zero */
 435             for (i = 0; i < mdatoms->homenr; i++)
 436             {
 437                 if (mdatoms->ptype[i] == eptVSite ||
 438                     mdatoms->ptype[i] == eptShell)
 439                 {
 440                     clear_rvec(v[i]);
 441                 }
 442                 else if (mdatoms->cFREEZE)
 443                 {
 444                     for (m = 0; m < DIM; m++)
 445                     {
 446                         if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 447                         {
 448                             v[i][m] = 0;
 449                         }
 450                     }
 451                 }
 452             }
 453         }
 454
 455         if (constr)
 456         {
 457             /* Constrain the initial coordinates and velocities */
 458             do_constrain_first(fplog, constr, ir, mdatoms,
 459                                state->natoms,
 460                                state->x.arrayRefWithPadding(),
 461                                state->v.arrayRefWithPadding(),
 462                                state->box, state->lambda[efptBONDED]);
 463         }
 464         if (vsite)
 465         {
 466             /* Construct the virtual sites for the initial configuration */
 467             construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
 468                              top.idef.iparams, top.idef.il,
 469                              fr->ePBC, fr->bMolPBC, cr, state->box);
 470         }
 471     }
 472
 473     if (ir->efep != efepNO)
 474     {
 475         /* Set free energy calculation frequency as the greatest common
 476          * denominator of nstdhdl and repl_ex_nst. */
 477         nstfep = ir->fepvals->nstdhdl;
 478         if (ir->bExpanded)
 479         {
 480             nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
 481         }
 482         if (useReplicaExchange)
 483         {
 484             nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
 485         }
 486     }
 487
 488     /* Be REALLY careful about what flags you set here. You CANNOT assume
 489      * this is the first step, since we might be restarting from a checkpoint,
 490      * and in that case we should not do any modifications to the state.
 491      */
 492     bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
 493
 494     // When restarting from a checkpoint, it can be appropriate to
 495     // initialize ekind from quantities in the checkpoint. Otherwise,
 496     // compute_globals must initialize ekind before the simulation
 497     // starts/restarts. However, only the master rank knows what was
 498     // found in the checkpoint file, so we have to communicate in
 499     // order to coordinate the restart.
 500     //
 501     // TODO Consider removing this communication if/when checkpoint
 502     // reading directly follows .tpr reading, because all ranks can
 503     // agree on hasReadEkinState at that time.
 504     bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
 505     if (PAR(cr))
 506     {
 507         gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr);
 508     }
 509     if (hasReadEkinState)
 510     {
 511         restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 512     }
 513
 514     unsigned int cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
 515                                | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
 516                                | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
 517                                | (hasReadEkinState ? CGLO_READEKIN : 0));
 518
 519     bSumEkinhOld = FALSE;
 520
 521     t_vcm vcm(top_global->groups, *ir);
 522     reportComRemovalInfo(fplog, vcm);
 523
 524     /* To minimize communication, compute_globals computes the COM velocity
 525      * and the kinetic energy for the velocities without COM motion removed.
 526      * Thus to get the kinetic energy without the COM contribution, we need
 527      * to call compute_globals twice.
 528      */
 529     for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
 530     {
 531         unsigned int cglo_flags_iteration = cglo_flags;
 532         if (bStopCM && cgloIteration == 0)
 533         {
 534             cglo_flags_iteration |= CGLO_STOPCM;
 535             cglo_flags_iteration &= ~CGLO_TEMPERATURE;
 536         }
 537         compute_globals(gstat, cr, ir, fr, ekind,
 538                         state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
 539                         mdatoms, nrnb, &vcm,
 540                         nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 541                         constr, &nullSignaller, state->box,
 542                         &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
 543                         | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
 544         if (cglo_flags_iteration & CGLO_STOPCM)
 545         {
 546             /* At initialization, do not pass x with acceleration-correction mode
 547              * to avoid (incorrect) correction of the initial coordinates.
 548              */
 549             rvec *xPtr = nullptr;
 550             if (vcm.mode != ecmLINEAR_ACCELERATION_CORRECTION)
 551             {
 552                 xPtr = state->x.rvec_array();
 553             }
 554             process_and_stopcm_grp(fplog, &vcm, *mdatoms, xPtr, state->v.rvec_array());
 555             inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
 556         }
 557     }
 558     checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
 559                                     top_global, &top, state->x.rvec_array(), state->box,
 560                                     &shouldCheckNumberOfBondedInteractions);
 561     if (ir->eI == eiVVAK)
 562     {
 563         /* a second call to get the half step temperature initialized as well */
 564         /* we do the same call as above, but turn the pressure off -- internally to
 565            compute_globals, this is recognized as a velocity verlet half-step
 566            kinetic energy calculation.  This minimized excess variables, but
 567            perhaps loses some logic?*/
 568
 569         compute_globals(gstat, cr, ir, fr, ekind,
 570                         state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
 571                         mdatoms, nrnb, &vcm,
 572                         nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 573                         constr, &nullSignaller, state->box,
 574                         nullptr, &bSumEkinhOld,
 575                         cglo_flags & ~CGLO_PRESSURE);
 576     }
 577
 578     /* Calculate the initial half step temperature, and save the ekinh_old */
 579     if (startingBehavior == StartingBehavior::NewSimulation)
 580     {
 581         for (i = 0; (i < ir->opts.ngtc); i++)
 582         {
 583             copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
 584         }
 585     }
 586
 587     /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
 588        temperature control */
 589     auto trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
 590
 591     if (MASTER(cr))
 592     {
 593         if (!ir->bContinuation)
 594         {
 595             if (constr && ir->eConstrAlg == econtLINCS)
 596             {
 597                 fprintf(fplog,
 598                         "RMS relative constraint deviation after constraining: %.2e\n",
 599                         constr->rmsd());
 600             }
 601             if (EI_STATE_VELOCITY(ir->eI))
 602             {
 603                 real temp = enerd->term[F_TEMP];
 604                 if (ir->eI != eiVV)
 605                 {
 606                     /* Result of Ekin averaged over velocities of -half
 607                      * and +half step, while we only have -half step here.
 608                      */
 609                     temp *= 2;
 610                 }
 611                 fprintf(fplog, "Initial temperature: %g K\n", temp);
 612             }
 613         }
 614
 615         char tbuf[20];
 616         fprintf(stderr, "starting mdrun '%s'\n",
 617                 *(top_global->name));
 618         if (ir->nsteps >= 0)
 619         {
 620             sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
 621         }
 622         else
 623         {
 624             sprintf(tbuf, "%s", "infinite");
 625         }
 626         if (ir->init_step > 0)
 627         {
 628             fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 629                     gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
 630                     gmx_step_str(ir->init_step, sbuf2),
 631                     ir->init_step*ir->delta_t);
 632         }
 633         else
 634         {
 635             fprintf(stderr, "%s steps, %s ps.\n",
 636                     gmx_step_str(ir->nsteps, sbuf), tbuf);
 637         }
 638         fprintf(fplog, "\n");
 639     }
 640
 641     walltime_accounting_start_time(walltime_accounting);
 642     wallcycle_start(wcycle, ewcRUN);
 643     print_start(fplog, cr, walltime_accounting, "mdrun");
 644
 645 #if GMX_FAHCORE
 646     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
 647     int chkpt_ret = fcCheckPointParallel( cr->nodeid,
 648                                           NULL, 0);
 649     if (chkpt_ret == 0)
 650     {
 651         gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
 652     }
 653 #endif
 654
 655     /***********************************************************
 656      *
 657      *             Loop over MD steps
 658      *
 659      ************************************************************/
 660
 661     bFirstStep       = TRUE;
 662     /* Skip the first Nose-Hoover integration when we get the state from tpx */
 663     bInitStep        = startingBehavior == StartingBehavior::NewSimulation || EI_VV(ir->eI);
 664     bSumEkinhOld     = FALSE;
 665     bExchanged       = FALSE;
 666     bNeedRepartition = FALSE;
 667
 668     bool simulationsShareState = false;
 669     int  nstSignalComm         = nstglobalcomm;
 670     {
 671         // TODO This implementation of ensemble orientation restraints is nasty because
 672         // a user can't just do multi-sim with single-sim orientation restraints.
 673         bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
 674         bool awhUsesMultiSim         = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
 675
 676         // Replica exchange, ensemble restraints and AWH need all
 677         // simulations to remain synchronized, so they need
 678         // checkpoints and stop conditions to act on the same step, so
 679         // the propagation of such signals must take place between
 680         // simulations, not just within simulations.
 681         // TODO: Make algorithm initializers set these flags.
 682         simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
 683
 684         if (simulationsShareState)
 685         {
 686             // Inter-simulation signal communication does not need to happen
 687             // often, so we use a minimum of 200 steps to reduce overhead.
 688             const int c_minimumInterSimulationSignallingInterval = 200;
 689             nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
 690         }
 691     }
 692
 693     auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
 694                 compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
 695                 MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
 696                 mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
 697
 698     auto checkpointHandler = std::make_unique<CheckpointHandler>(
 699                 compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]),
 700                 simulationsShareState, ir->nstlist == 0, MASTER(cr),
 701                 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
 702
 703     const bool resetCountersIsLocal = true;
 704     auto       resetHandler         = std::make_unique<ResetHandler>(
 705                 compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]), !resetCountersIsLocal,
 706                 ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
 707                 mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
 708
 709     const DDBalanceRegionHandler ddBalanceRegionHandler(cr);
 710
 711     step     = ir->init_step;
 712     step_rel = 0;
 713
 714     // TODO extract this to new multi-simulation module
 715     if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
 716     {
 717         if (!multisim_int_all_are_equal(ms, ir->nsteps))
 718         {
 719             GMX_LOG(mdlog.warning).appendText(
 720                     "Note: The number of steps is not consistent across multi simulations,\n"
 721                     "but we are proceeding anyway!");
 722         }
 723         if (!multisim_int_all_are_equal(ms, ir->init_step))
 724         {
 725             GMX_LOG(mdlog.warning).appendText(
 726                     "Note: The initial step is not consistent across multi simulations,\n"
 727                     "but we are proceeding anyway!");
 728         }
 729     }
 730
 731     /* and stop now if we should */
 732     bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
 733     while (!bLastStep)
 734     {
 735
 736         /* Determine if this is a neighbor search step */
 737         bNStList = (ir->nstlist > 0  && step % ir->nstlist == 0);
 738
 739         if (bPMETune && bNStList)
 740         {
 741             /* PME grid + cut-off optimization with GPUs or PME nodes */
 742             pme_loadbal_do(pme_loadbal, cr,
 743                            (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
 744                            fplog, mdlog,
 745                            *ir, fr, state->box, state->x,
 746                            wcycle,
 747                            step, step_rel,
 748                            &bPMETunePrinting);
 749         }
 750
 751         wallcycle_start(wcycle, ewcSTEP);
 752
 753         bLastStep = (step_rel == ir->nsteps);
 754         t         = t0 + step*ir->delta_t;
 755
 756         // TODO Refactor this, so that nstfep does not need a default value of zero
 757         if (ir->efep != efepNO || ir->bSimTemp)
 758         {
 759             /* find and set the current lambdas */
 760             setCurrentLambdasLocal(step, ir->fepvals, lam0, state->lambda, state->fep_state);
 761
 762             bDoDHDL      = do_per_step(step, ir->fepvals->nstdhdl);
 763             bDoFEP       = ((ir->efep != efepNO) && do_per_step(step, nstfep));
 764             bDoExpanded  = (do_per_step(step, ir->expandedvals->nstexpanded)
 765                             && (ir->bExpanded) && (step > 0) &&
 766                             (startingBehavior == StartingBehavior::NewSimulation));
 767         }
 768
 769         bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
 770                      do_per_step(step, replExParams.exchangeInterval));
 771
 772         if (doSimulatedAnnealing)
 773         {
 774             update_annealing_target_temp(ir, t, &upd);
 775         }
 776
 777         /* Stop Center of Mass motion */
 778         bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
 779
 780         /* Determine whether or not to do Neighbour Searching */
 781         bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
 782
 783         bLastStep = bLastStep || stopHandler->stoppingAfterCurrentStep(bNS);
 784
 785         /* do_log triggers energy and virial calculation. Because this leads
 786          * to different code paths, forces can be different. Thus for exact
 787          * continuation we should avoid extra log output.
 788          * Note that the || bLastStep can result in non-exact continuation
 789          * beyond the last step. But we don't consider that to be an issue.
 790          */
 791         do_log     =
 792             (do_per_step(step, ir->nstlog) ||
 793              (bFirstStep && startingBehavior == StartingBehavior::NewSimulation) ||
 794              bLastStep);
 795         do_verbose = mdrunOptions.verbose &&
 796             (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
 797
 798         // Copy velocities from the GPU when needed:
 799         // - On search steps to keep copy on host (device buffers are reinitialized).
 800         // - When needed for the output.
 801         if (useGpuForUpdate && !bFirstStep)
 802         {
 803             if (bNS || do_per_step(step, ir->nstvout))
 804             {
 805                 stateGpu->copyVelocitiesFromGpu(state->v, StatePropagatorDataGpu::AtomLocality::Local);
 806                 stateGpu->waitVelocitiesReadyOnHost(StatePropagatorDataGpu::AtomLocality::Local);
 807             }
 808         }
 809
 810
 811         if (bNS && !(bFirstStep && ir->bContinuation))
 812         {
 813             bMasterState = FALSE;
 814             /* Correct the new box if it is too skewed */
 815             if (inputrecDynamicBox(ir))
 816             {
 817                 if (correct_box(fplog, step, state->box, graph))
 818                 {
 819                     bMasterState = TRUE;
 820                 }
 821             }
 822             if (DOMAINDECOMP(cr) && bMasterState)
 823             {
 824                 dd_collect_state(cr->dd, state, state_global);
 825             }
 826
 827             if (DOMAINDECOMP(cr))
 828             {
 829                 /* Repartition the domain decomposition */
 830                 dd_partition_system(fplog, mdlog, step, cr,
 831                                     bMasterState, nstglobalcomm,
 832                                     state_global, *top_global, ir, imdSession,
 833                                     pull_work,
 834                                     state, &f, mdAtoms, &top, fr,
 835                                     vsite, constr,
 836                                     nrnb, wcycle,
 837                                     do_verbose && !bPMETunePrinting);
 838                 shouldCheckNumberOfBondedInteractions = true;
 839                 upd.setNumAtoms(state->natoms);
 840             }
 841         }
 842
 843         if (MASTER(cr) && do_log)
 844         {
 845             energyOutput.printHeader(fplog, step, t); /* can we improve the information printed here? */
 846         }
 847
 848         if (ir->efep != efepNO)
 849         {
 850             update_mdatoms(mdatoms, state->lambda[efptMASS]);
 851         }
 852
 853         if (bExchanged)
 854         {
 855
 856             /* We need the kinetic energy at minus the half step for determining
 857              * the full step kinetic energy and possibly for T-coupling.*/
 858             /* This may not be quite working correctly yet . . . . */
 859             compute_globals(gstat, cr, ir, fr, ekind,
 860                             state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
 861                             mdatoms, nrnb, &vcm,
 862                             wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
 863                             constr, &nullSignaller, state->box,
 864                             &totalNumberOfBondedInteractions, &bSumEkinhOld,
 865                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
 866             checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
 867                                             top_global, &top, state->x.rvec_array(), state->box,
 868                                             &shouldCheckNumberOfBondedInteractions);
 869         }
 870         clear_mat(force_vir);
 871
 872         checkpointHandler->decideIfCheckpointingThisStep(bNS, bFirstStep, bLastStep);
 873
 874         /* Determine the energy and pressure:
 875          * at nstcalcenergy steps and at energy output steps (set below).
 876          */
 877         if (EI_VV(ir->eI) && (!bInitStep))
 878         {
 879             bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
 880             bCalcVir      = bCalcEnerStep ||
 881                 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
 882         }
 883         else
 884         {
 885             bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
 886             bCalcVir      = bCalcEnerStep ||
 887                 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
 888         }
 889         bCalcEner = bCalcEnerStep;
 890
 891         do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
 892
 893         if (do_ene || do_log || bDoReplEx)
 894         {
 895             bCalcVir  = TRUE;
 896             bCalcEner = TRUE;
 897         }
 898
 899         /* Do we need global communication ? */
 900         bGStat = (bCalcVir || bCalcEner || bStopCM ||
 901                   do_per_step(step, nstglobalcomm) ||
 902                   (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
 903
 904         force_flags = (GMX_FORCE_STATECHANGED |
 905                        ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0) |
 906                        GMX_FORCE_ALLFORCES |
 907                        (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
 908                        (bCalcEner ? GMX_FORCE_ENERGY : 0) |
 909                        (bDoFEP ? GMX_FORCE_DHDL : 0)
 910                        );
 911
 912         if (shellfc)
 913         {
 914             /* Now is the time to relax the shells */
 915             relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
 916                                 enforcedRotation, step,
 917                                 ir, imdSession, pull_work, bNS, force_flags, &top,
 918                                 constr, enerd, fcd,
 919                                 state->natoms,
 920                                 state->x.arrayRefWithPadding(),
 921                                 state->v.arrayRefWithPadding(),
 922                                 state->box,
 923                                 state->lambda,
 924                                 &state->hist,
 925                                 f.arrayRefWithPadding(), force_vir, mdatoms,
 926                                 nrnb, wcycle, graph,
 927                                 shellfc, fr, runScheduleWork, t, mu_tot,
 928                                 vsite,
 929                                 ddBalanceRegionHandler);
 930         }
 931         else
 932         {
 933             /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
 934                (or the AWH update will be performed twice for one step when continuing). It would be best to
 935                call this update function from do_md_trajectory_writing but that would occur after do_force.
 936                One would have to divide the update_awh function into one function applying the AWH force
 937                and one doing the AWH bias update. The update AWH bias function could then be called after
 938                do_md_trajectory_writing (then containing update_awh_history).
 939                The checkpointing will in the future probably moved to the start of the md loop which will
 940                rid of this issue. */
 941             if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
 942             {
 943                 awh->updateHistory(state_global->awhHistory.get());
 944             }
 945
 946             /* The coordinates (x) are shifted (to get whole molecules)
 947              * in do_force.
 948              * This is parallellized as well, and does communication too.
 949              * Check comments in sim_util.c
 950              */
 951             do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession,
 952                      pull_work,
 953                      step, nrnb, wcycle, &top,
 954                      state->box, state->x.arrayRefWithPadding(), &state->hist,
 955                      f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
 956                      state->lambda, graph,
 957                      fr, runScheduleWork, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
 958                      (bNS ? GMX_FORCE_NS : 0) | force_flags,
 959                      ddBalanceRegionHandler);
 960         }
 961
 962         // VV integrators do not need the following velocity half step
 963         // if it is the first step after starting from a checkpoint.
 964         // That is, the half step is needed on all other steps, and
 965         // also the first step when starting from a .tpr file.
 966         if (EI_VV(ir->eI) && (!bFirstStep || startingBehavior == StartingBehavior::NewSimulation))
 967         /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
 968         {
 969             rvec *vbuf = nullptr;
 970
 971             wallcycle_start(wcycle, ewcUPDATE);
 972             if (ir->eI == eiVV && bInitStep)
 973             {
 974                 /* if using velocity verlet with full time step Ekin,
 975                  * take the first half step only to compute the
 976                  * virial for the first step. From there,
 977                  * revert back to the initial coordinates
 978                  * so that the input is actually the initial step.
 979                  */
 980                 snew(vbuf, state->natoms);
 981                 copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
 982             }
 983             else
 984             {
 985                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
 986                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
 987             }
 988
 989             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
 990                           ekind, M, &upd, etrtVELOCITY1,
 991                           cr, constr);
 992
 993             wallcycle_stop(wcycle, ewcUPDATE);
 994             constrain_velocities(step, nullptr,
 995                                  state,
 996                                  shake_vir,
 997                                  constr,
 998                                  bCalcVir, do_log, do_ene);
 999             wallcycle_start(wcycle, ewcUPDATE);
1000             /* if VV, compute the pressure and constraints */
1001             /* For VV2, we strictly only need this if using pressure
1002              * control, but we really would like to have accurate pressures
1003              * printed out.
1004              * Think about ways around this in the future?
1005              * For now, keep this choice in comments.
1006              */
1007             /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
1008             /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
1009             bPres = TRUE;
1010             bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
1011             if (bCalcEner && ir->eI == eiVVAK)
1012             {
1013                 bSumEkinhOld = TRUE;
1014             }
1015             /* for vv, the first half of the integration actually corresponds to the previous step.
1016                So we need information from the last step in the first half of the integration */
1017             if (bGStat || do_per_step(step-1, nstglobalcomm))
1018             {
1019                 wallcycle_stop(wcycle, ewcUPDATE);
1020                 compute_globals(gstat, cr, ir, fr, ekind,
1021                                 state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
1022                                 mdatoms, nrnb, &vcm,
1023                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1024                                 constr, &nullSignaller, state->box,
1025                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1026                                 (bGStat ? CGLO_GSTAT : 0)
1027                                 | (bCalcEner ? CGLO_ENERGY : 0)
1028                                 | (bTemp ? CGLO_TEMPERATURE : 0)
1029                                 | (bPres ? CGLO_PRESSURE : 0)
1030                                 | (bPres ? CGLO_CONSTRAINT : 0)
1031                                 | (bStopCM ? CGLO_STOPCM : 0)
1032                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1033                                 | CGLO_SCALEEKIN
1034                                 );
1035                 /* explanation of above:
1036                    a) We compute Ekin at the full time step
1037                    if 1) we are using the AveVel Ekin, and it's not the
1038                    initial step, or 2) if we are using AveEkin, but need the full
1039                    time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1040                    b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1041                    EkinAveVel because it's needed for the pressure */
1042                 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1043                                                 top_global, &top, state->x.rvec_array(), state->box,
1044                                                 &shouldCheckNumberOfBondedInteractions);
1045                 if (bStopCM)
1046                 {
1047                     process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(), state->v.rvec_array());
1048                     inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
1049                 }
1050                 wallcycle_start(wcycle, ewcUPDATE);
1051             }
1052             /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1053             if (!bInitStep)
1054             {
1055                 if (bTrotter)
1056                 {
1057                     m_add(force_vir, shake_vir, total_vir);     /* we need the un-dispersion corrected total vir here */
1058                     trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
1059
1060                     /* TODO This is only needed when we're about to write
1061                      * a checkpoint, because we use it after the restart
1062                      * (in a kludge?). But what should we be doing if
1063                      * the startingBehavior is NewSimulation or bInitStep are true? */
1064                     if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1065                     {
1066                         copy_mat(shake_vir, state->svir_prev);
1067                         copy_mat(force_vir, state->fvir_prev);
1068                     }
1069                     if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
1070                     {
1071                         /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1072                         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
1073                         enerd->term[F_EKIN] = trace(ekind->ekin);
1074                     }
1075                 }
1076                 else if (bExchanged)
1077                 {
1078                     wallcycle_stop(wcycle, ewcUPDATE);
1079                     /* We need the kinetic energy at minus the half step for determining
1080                      * the full step kinetic energy and possibly for T-coupling.*/
1081                     /* This may not be quite working correctly yet . . . . */
1082                     compute_globals(gstat, cr, ir, fr, ekind,
1083                                     state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
1084                                     mdatoms, nrnb, &vcm,
1085                                     wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1086                                     constr, &nullSignaller, state->box,
1087                                     nullptr, &bSumEkinhOld,
1088                                     CGLO_GSTAT | CGLO_TEMPERATURE);
1089                     wallcycle_start(wcycle, ewcUPDATE);
1090                 }
1091             }
1092             /* if it's the initial step, we performed this first step just to get the constraint virial */
1093             if (ir->eI == eiVV && bInitStep)
1094             {
1095                 copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
1096                 sfree(vbuf);
1097             }
1098             wallcycle_stop(wcycle, ewcUPDATE);
1099         }
1100
1101         /* compute the conserved quantity */
1102         if (EI_VV(ir->eI))
1103         {
1104             saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1105             if (ir->eI == eiVV)
1106             {
1107                 last_ekin = enerd->term[F_EKIN];
1108             }
1109             if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1110             {
1111                 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1112             }
1113             /* sum up the foreign energy and dhdl terms for vv.  currently done every step so that dhdl is correct in the .edr */
1114             if (ir->efep != efepNO)
1115             {
1116                 sum_dhdl(enerd, state->lambda, *ir->fepvals);
1117             }
1118         }
1119
1120         /* ########  END FIRST UPDATE STEP  ############## */
1121         /* ########  If doing VV, we now have v(dt) ###### */
1122         if (bDoExpanded)
1123         {
1124             /* perform extended ensemble sampling in lambda - we don't
1125                actually move to the new state before outputting
1126                statistics, but if performing simulated tempering, we
1127                do update the velocities and the tau_t. */
1128
1129             lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v.rvec_array(), mdatoms);
1130             /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1131             if (MASTER(cr))
1132             {
1133                 copy_df_history(state_global->dfhist, state->dfhist);
1134             }
1135         }
1136
1137         /* Now we have the energies and forces corresponding to the
1138          * coordinates at time t. We must output all of this before
1139          * the update.
1140          */
1141         do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1142                                  ir, state, state_global, observablesHistory,
1143                                  top_global, fr,
1144                                  outf, energyOutput, ekind, f,
1145                                  checkpointHandler->isCheckpointingStep(),
1146                                  bRerunMD, bLastStep,
1147                                  mdrunOptions.writeConfout,
1148                                  bSumEkinhOld);
1149         /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1150         bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x.rvec_array(), t);
1151
1152         /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1153         if (startingBehavior != StartingBehavior::NewSimulation &&
1154             (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1155         {
1156             copy_mat(state->svir_prev, shake_vir);
1157             copy_mat(state->fvir_prev, force_vir);
1158         }
1159
1160         stopHandler->setSignal();
1161         resetHandler->setSignal(walltime_accounting);
1162
1163         if (bGStat || !PAR(cr))
1164         {
1165             /* In parallel we only have to check for checkpointing in steps
1166              * where we do global communication,
1167              *  otherwise the other nodes don't know.
1168              */
1169             checkpointHandler->setSignal(walltime_accounting);
1170         }
1171
1172         /* #########   START SECOND UPDATE STEP ################# */
1173
1174         /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1175            in preprocessing */
1176
1177         if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1178         {
1179             gmx_bool bIfRandomize;
1180             bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
1181             /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1182             if (constr && bIfRandomize)
1183             {
1184                 constrain_velocities(step, nullptr,
1185                                      state,
1186                                      tmp_vir,
1187                                      constr,
1188                                      bCalcVir, do_log, do_ene);
1189             }
1190         }
1191         /* Box is changed in update() when we do pressure coupling,
1192          * but we should still use the old box for energy corrections and when
1193          * writing it to the energy file, so it matches the trajectory files for
1194          * the same timestep above. Make a copy in a separate array.
1195          */
1196         copy_mat(state->box, lastbox);
1197
1198         dvdl_constr = 0;
1199
1200         wallcycle_start(wcycle, ewcUPDATE);
1201         /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1202         if (bTrotter)
1203         {
1204             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1205             /* We can only do Berendsen coupling after we have summed
1206              * the kinetic energy or virial. Since the happens
1207              * in global_state after update, we should only do it at
1208              * step % nstlist = 1 with bGStatEveryStep=FALSE.
1209              */
1210         }
1211         else
1212         {
1213             update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1214             update_pcouple_before_coordinates(fplog, step, ir, state,
1215                                               parrinellorahmanMu, M,
1216                                               bInitStep);
1217         }
1218
1219         if (EI_VV(ir->eI))
1220         {
1221             /* velocity half-step update */
1222             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1223                           ekind, M, &upd, etrtVELOCITY2,
1224                           cr, constr);
1225         }
1226
1227         /* Above, initialize just copies ekinh into ekin,
1228          * it doesn't copy position (for VV),
1229          * and entire integrator for MD.
1230          */
1231
1232         if (ir->eI == eiVVAK)
1233         {
1234             cbuf.resize(state->x.size());
1235             std::copy(state->x.begin(), state->x.end(), cbuf.begin());
1236         }
1237
1238         /* With leap-frog type integrators we compute the kinetic energy
1239          * at a whole time step as the average of the half-time step kinetic
1240          * energies of two subsequent steps. Therefore we need to compute the
1241          * half step kinetic energy also if we need energies at the next step.
1242          */
1243         const bool needHalfStepKineticEnergy = (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm));
1244
1245         if (useGpuForUpdate)
1246         {
1247             if (bNS)
1248             {
1249                 integrator->set(stateGpu->getCoordinates(), stateGpu->getVelocities(), stateGpu->getForces(),
1250                                 top.idef, *mdatoms, ekind->ngtc);
1251                 t_pbc pbc;
1252                 set_pbc(&pbc, epbcXYZ, state->box);
1253                 integrator->setPbc(&pbc);
1254
1255                 // Copy data to the GPU after buffers might have being reinitialized
1256                 stateGpu->copyVelocitiesToGpu(state->v, StatePropagatorDataGpu::AtomLocality::Local);
1257             }
1258
1259             stateGpu->copyCoordinatesToGpu(ArrayRef<RVec>(state->x), StatePropagatorDataGpu::AtomLocality::All);
1260             stateGpu->copyForcesToGpu(ArrayRef<RVec>(f), StatePropagatorDataGpu::AtomLocality::All);
1261
1262             // TODO: Use StepWorkload fields.
1263             bool useGpuFBufferOps = simulationWork.useGpuBufferOps && !(bCalcVir || bCalcEner);
1264
1265             bool doTempCouple       = (ir->etc != etcNO && do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
1266             bool doParrinelloRahman = (ir->epc == epcPARRINELLORAHMAN && do_per_step(step + ir->nstpcouple - 1, ir->nstpcouple));
1267
1268             // This applies Leap-Frog, LINCS and SETTLE in succession
1269             integrator->integrate(stateGpu->getForcesReadyOnDeviceEvent(StatePropagatorDataGpu::AtomLocality::Local, useGpuFBufferOps),
1270                                   ir->delta_t, true, bCalcVir, shake_vir,
1271                                   doTempCouple, ekind->tcstat,
1272                                   doParrinelloRahman, ir->nstpcouple*ir->delta_t, M);
1273             stateGpu->copyCoordinatesFromGpu(ArrayRef<RVec>(state->x), StatePropagatorDataGpu::AtomLocality::All);
1274
1275             // Copy velocities D2H after update if:
1276             // - Globals are computed this step (includes the energy output steps).
1277             // - Temperature is needed for the next step.
1278             if (bGStat || needHalfStepKineticEnergy)
1279             {
1280                 stateGpu->copyVelocitiesFromGpu(state->v, StatePropagatorDataGpu::AtomLocality::Local);
1281                 stateGpu->waitVelocitiesReadyOnHost(StatePropagatorDataGpu::AtomLocality::Local);
1282             }
1283
1284             // TODO: replace with stateGpu->waitForCopyCoordinatesFromGpu(...)
1285             integrator->waitCoordinatesReadyOnDevice();
1286         }
1287         else
1288         {
1289             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1290                           ekind, M, &upd, etrtPOSITION, cr, constr);
1291
1292             wallcycle_stop(wcycle, ewcUPDATE);
1293
1294             constrain_coordinates(step, &dvdl_constr, state,
1295                                   shake_vir,
1296                                   &upd, constr,
1297                                   bCalcVir, do_log, do_ene);
1298
1299             update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
1300                                   cr, nrnb, wcycle, &upd, constr, do_log, do_ene);
1301             finish_update(ir, mdatoms,
1302                           state, graph,
1303                           nrnb, wcycle, &upd, constr);
1304         }
1305
1306         if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1307         {
1308             updatePrevStepPullCom(pull_work, state);
1309         }
1310
1311         if (ir->eI == eiVVAK)
1312         {
1313             /* erase F_EKIN and F_TEMP here? */
1314             /* just compute the kinetic energy at the half step to perform a trotter step */
1315             compute_globals(gstat, cr, ir, fr, ekind,
1316                             state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
1317                             mdatoms, nrnb, &vcm,
1318                             wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1319                             constr, &nullSignaller, lastbox,
1320                             nullptr, &bSumEkinhOld,
1321                             (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1322                             );
1323             wallcycle_start(wcycle, ewcUPDATE);
1324             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1325             /* now we know the scaling, we can compute the positions again */
1326             std::copy(cbuf.begin(), cbuf.end(), state->x.begin());
1327
1328             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1329                           ekind, M, &upd, etrtPOSITION, cr, constr);
1330             wallcycle_stop(wcycle, ewcUPDATE);
1331
1332             /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1333             /* are the small terms in the shake_vir here due
1334              * to numerical errors, or are they important
1335              * physically? I'm thinking they are just errors, but not completely sure.
1336              * For now, will call without actually constraining, constr=NULL*/
1337             finish_update(ir, mdatoms,
1338                           state, graph,
1339                           nrnb, wcycle, &upd, nullptr);
1340         }
1341         if (EI_VV(ir->eI))
1342         {
1343             /* this factor or 2 correction is necessary
1344                because half of the constraint force is removed
1345                in the vv step, so we have to double it.  See
1346                the Redmine issue #1255.  It is not yet clear
1347                if the factor of 2 is exact, or just a very
1348                good approximation, and this will be
1349                investigated.  The next step is to see if this
1350                can be done adding a dhdl contribution from the
1351                rattle step, but this is somewhat more
1352                complicated with the current code. Will be
1353                investigated, hopefully for 4.6.3. However,
1354                this current solution is much better than
1355                having it completely wrong.
1356              */
1357             enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1358         }
1359         else
1360         {
1361             enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1362         }
1363
1364         if (vsite != nullptr)
1365         {
1366             wallcycle_start(wcycle, ewcVSITECONSTR);
1367             if (graph != nullptr)
1368             {
1369                 shift_self(graph, state->box, state->x.rvec_array());
1370             }
1371             construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
1372                              top.idef.iparams, top.idef.il,
1373                              fr->ePBC, fr->bMolPBC, cr, state->box);
1374
1375             if (graph != nullptr)
1376             {
1377                 unshift_self(graph, state->box, state->x.rvec_array());
1378             }
1379             wallcycle_stop(wcycle, ewcVSITECONSTR);
1380         }
1381
1382         /* ############## IF NOT VV, Calculate globals HERE  ############ */
1383         /* With Leap-Frog we can skip compute_globals at
1384          * non-communication steps, but we need to calculate
1385          * the kinetic energy one step before communication.
1386          */
1387         {
1388             // Organize to do inter-simulation signalling on steps if
1389             // and when algorithms require it.
1390             bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1391
1392             if (bGStat || needHalfStepKineticEnergy || doInterSimSignal)
1393             {
1394                 // Since we're already communicating at this step, we
1395                 // can propagate intra-simulation signals. Note that
1396                 // check_nstglobalcomm has the responsibility for
1397                 // choosing the value of nstglobalcomm that is one way
1398                 // bGStat becomes true, so we can't get into a
1399                 // situation where e.g. checkpointing can't be
1400                 // signalled.
1401                 bool                doIntraSimSignal = true;
1402                 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1403
1404                 compute_globals(gstat, cr, ir, fr, ekind,
1405                                 state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
1406                                 mdatoms, nrnb, &vcm,
1407                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1408                                 constr, &signaller,
1409                                 lastbox,
1410                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1411                                 (bGStat ? CGLO_GSTAT : 0)
1412                                 | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1413                                 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1414                                 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1415                                 | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
1416                                 | CGLO_CONSTRAINT
1417                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1418                                 );
1419                 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1420                                                 top_global, &top, state->x.rvec_array(), state->box,
1421                                                 &shouldCheckNumberOfBondedInteractions);
1422                 if (!EI_VV(ir->eI) && bStopCM)
1423                 {
1424                     process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(), state->v.rvec_array());
1425                     inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
1426                 }
1427             }
1428         }
1429
1430         /* #############  END CALC EKIN AND PRESSURE ################# */
1431
1432         /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1433            the virial that should probably be addressed eventually. state->veta has better properies,
1434            but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1435            generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
1436
1437         if (ir->efep != efepNO && !EI_VV(ir->eI))
1438         {
1439             /* Sum up the foreign energy and dhdl terms for md and sd.
1440                Currently done every step so that dhdl is correct in the .edr */
1441             sum_dhdl(enerd, state->lambda, *ir->fepvals);
1442         }
1443
1444         update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1445                                          pres, force_vir, shake_vir,
1446                                          parrinellorahmanMu,
1447                                          state, nrnb, &upd);
1448
1449         /* ################# END UPDATE STEP 2 ################# */
1450         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
1451
1452         /* The coordinates (x) were unshifted in update */
1453         if (!bGStat)
1454         {
1455             /* We will not sum ekinh_old,
1456              * so signal that we still have to do it.
1457              */
1458             bSumEkinhOld = TRUE;
1459         }
1460
1461         if (bCalcEner)
1462         {
1463             /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
1464
1465             /* use the directly determined last velocity, not actually the averaged half steps */
1466             if (bTrotter && ir->eI == eiVV)
1467             {
1468                 enerd->term[F_EKIN] = last_ekin;
1469             }
1470             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1471
1472             if (integratorHasConservedEnergyQuantity(ir))
1473             {
1474                 if (EI_VV(ir->eI))
1475                 {
1476                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1477                 }
1478                 else
1479                 {
1480                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1481                 }
1482             }
1483             /* #########  END PREPARING EDR OUTPUT  ###########  */
1484         }
1485
1486         /* Output stuff */
1487         if (MASTER(cr))
1488         {
1489             if (fplog && do_log && bDoExpanded)
1490             {
1491                 /* only needed if doing expanded ensemble */
1492                 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1493                                           state_global->dfhist, state->fep_state, ir->nstlog, step);
1494             }
1495             if (bCalcEner)
1496             {
1497                 energyOutput.addDataAtEnergyStep(bDoDHDL, bCalcEnerStep,
1498                                                  t, mdatoms->tmass, enerd, state,
1499                                                  ir->fepvals, ir->expandedvals, lastbox,
1500                                                  shake_vir, force_vir, total_vir, pres,
1501                                                  ekind, mu_tot, constr);
1502             }
1503             else
1504             {
1505                 energyOutput.recordNonEnergyStep();
1506             }
1507
1508             gmx_bool do_dr  = do_per_step(step, ir->nstdisreout);
1509             gmx_bool do_or  = do_per_step(step, ir->nstorireout);
1510
1511             if (doSimulatedAnnealing)
1512             {
1513                 energyOutput.printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
1514             }
1515             if (do_log || do_ene || do_dr || do_or)
1516             {
1517                 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
1518                                                    do_log ? fplog : nullptr,
1519                                                    step, t,
1520                                                    fcd, awh.get());
1521             }
1522
1523             if (ir->bPull)
1524             {
1525                 pull_print_output(pull_work, step, t);
1526             }
1527
1528             if (do_per_step(step, ir->nstlog))
1529             {
1530                 if (fflush(fplog) != 0)
1531                 {
1532                     gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1533                 }
1534             }
1535         }
1536         if (bDoExpanded)
1537         {
1538             /* Have to do this part _after_ outputting the logfile and the edr file */
1539             /* Gets written into the state at the beginning of next loop*/
1540             state->fep_state = lamnew;
1541         }
1542         /* Print the remaining wall clock time for the run */
1543         if (isMasterSimMasterRank(ms, MASTER(cr)) &&
1544             (do_verbose || gmx_got_usr_signal()) &&
1545             !bPMETunePrinting)
1546         {
1547             if (shellfc)
1548             {
1549                 fprintf(stderr, "\n");
1550             }
1551             print_time(stderr, walltime_accounting, step, ir, cr);
1552         }
1553
1554         /* Ion/water position swapping.
1555          * Not done in last step since trajectory writing happens before this call
1556          * in the MD loop and exchanges would be lost anyway. */
1557         bNeedRepartition = FALSE;
1558         if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1559             do_per_step(step, ir->swap->nstswap))
1560         {
1561             bNeedRepartition = do_swapcoords(cr, step, t, ir, swap, wcycle,
1562                                              as_rvec_array(state->x.data()),
1563                                              state->box,
1564                                              MASTER(cr) && mdrunOptions.verbose,
1565                                              bRerunMD);
1566
1567             if (bNeedRepartition && DOMAINDECOMP(cr))
1568             {
1569                 dd_collect_state(cr->dd, state, state_global);
1570             }
1571         }
1572
1573         /* Replica exchange */
1574         bExchanged = FALSE;
1575         if (bDoReplEx)
1576         {
1577             bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
1578                                           state_global, enerd,
1579                                           state, step, t);
1580         }
1581
1582         if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1583         {
1584             dd_partition_system(fplog, mdlog, step, cr, TRUE, 1,
1585                                 state_global, *top_global, ir, imdSession,
1586                                 pull_work,
1587                                 state, &f, mdAtoms, &top, fr,
1588                                 vsite, constr,
1589                                 nrnb, wcycle, FALSE);
1590             shouldCheckNumberOfBondedInteractions = true;
1591             upd.setNumAtoms(state->natoms);
1592         }
1593
1594         bFirstStep             = FALSE;
1595         bInitStep              = FALSE;
1596
1597         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1598         /* With all integrators, except VV, we need to retain the pressure
1599          * at the current step for coupling at the next step.
1600          */
1601         if ((state->flags & (1U<<estPRES_PREV)) &&
1602             (bGStatEveryStep ||
1603              (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1604         {
1605             /* Store the pressure in t_state for pressure coupling
1606              * at the next MD step.
1607              */
1608             copy_mat(pres, state->pres_prev);
1609         }
1610
1611         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
1612
1613         if ( (membed != nullptr) && (!bLastStep) )
1614         {
1615             rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1616         }
1617
1618         cycles = wallcycle_stop(wcycle, ewcSTEP);
1619         if (DOMAINDECOMP(cr) && wcycle)
1620         {
1621             dd_cycles_add(cr->dd, cycles, ddCyclStep);
1622         }
1623
1624         /* increase the MD step number */
1625         step++;
1626         step_rel++;
1627
1628         resetHandler->resetCounters(
1629                 step, step_rel, mdlog, fplog, cr, fr->nbv.get(),
1630                 nrnb, fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
1631
1632         /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1633         imdSession->updateEnergyRecordAndSendPositionsAndEnergies(bInteractiveMDstep, step, bCalcEner);
1634
1635     }
1636     /* End of main MD loop */
1637
1638     /* Closing TNG files can include compressing data. Therefore it is good to do that
1639      * before stopping the time measurements. */
1640     mdoutf_tng_close(outf);
1641
1642     /* Stop measuring walltime */
1643     walltime_accounting_end_time(walltime_accounting);
1644
1645     if (!thisRankHasDuty(cr, DUTY_PME))
1646     {
1647         /* Tell the PME only node to finish */
1648         gmx_pme_send_finish(cr);
1649     }
1650
1651     if (MASTER(cr))
1652     {
1653         if (ir->nstcalcenergy > 0)
1654         {
1655             energyOutput.printAnnealingTemperatures(fplog, groups, &(ir->opts));
1656             energyOutput.printAverages(fplog, groups);
1657         }
1658     }
1659     done_mdoutf(outf);
1660
1661     if (bPMETune)
1662     {
1663         pme_loadbal_done(pme_loadbal, fplog, mdlog, fr->nbv->useGpu());
1664     }
1665
1666     done_shellfc(fplog, shellfc, step_rel);
1667
1668     if (useReplicaExchange && MASTER(cr))
1669     {
1670         print_replica_exchange_statistics(fplog, repl_ex);
1671     }
1672
1673     walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1674
1675     global_stat_destroy(gstat);
1676
1677 }