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