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