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