src/gromacs/mdrun/md.cpp

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