src/gromacs/mdrun/md.cpp

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