src/programs/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, 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 #include "gmxpre.h"
  38
  39 #include "md.h"
  40
  41 #include "config.h"
  42
  43 #include <math.h>
  44 #include <stdio.h>
  45 #include <stdlib.h>
  46
  47 #include <algorithm>
  48 #include <memory>
  49
  50 #include "thread_mpi/threads.h"
  51
  52 #include "gromacs/commandline/filenm.h"
  53 #include "gromacs/domdec/domdec.h"
  54 #include "gromacs/domdec/domdec_network.h"
  55 #include "gromacs/domdec/domdec_struct.h"
  56 #include "gromacs/essentialdynamics/edsam.h"
  57 #include "gromacs/ewald/pme.h"
  58 #include "gromacs/ewald/pme-load-balancing.h"
  59 #include "gromacs/fileio/trxio.h"
  60 #include "gromacs/gmxlib/network.h"
  61 #include "gromacs/gmxlib/nrnb.h"
  62 #include "gromacs/gpu_utils/gpu_utils.h"
  63 #include "gromacs/imd/imd.h"
  64 #include "gromacs/listed-forces/manage-threading.h"
  65 #include "gromacs/math/functions.h"
  66 #include "gromacs/math/utilities.h"
  67 #include "gromacs/math/vec.h"
  68 #include "gromacs/math/vectypes.h"
  69 #include "gromacs/mdlib/compute_io.h"
  70 #include "gromacs/mdlib/constr.h"
  71 #include "gromacs/mdlib/ebin.h"
  72 #include "gromacs/mdlib/force.h"
  73 #include "gromacs/mdlib/force_flags.h"
  74 #include "gromacs/mdlib/forcerec.h"
  75 #include "gromacs/mdlib/md_support.h"
  76 #include "gromacs/mdlib/mdatoms.h"
  77 #include "gromacs/mdlib/mdebin.h"
  78 #include "gromacs/mdlib/mdoutf.h"
  79 #include "gromacs/mdlib/mdrun.h"
  80 #include "gromacs/mdlib/mdsetup.h"
  81 #include "gromacs/mdlib/nb_verlet.h"
  82 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
  83 #include "gromacs/mdlib/ns.h"
  84 #include "gromacs/mdlib/shellfc.h"
  85 #include "gromacs/mdlib/sighandler.h"
  86 #include "gromacs/mdlib/sim_util.h"
  87 #include "gromacs/mdlib/simulationsignal.h"
  88 #include "gromacs/mdlib/tgroup.h"
  89 #include "gromacs/mdlib/trajectory_writing.h"
  90 #include "gromacs/mdlib/update.h"
  91 #include "gromacs/mdlib/vcm.h"
  92 #include "gromacs/mdlib/vsite.h"
  93 #include "gromacs/mdtypes/commrec.h"
  94 #include "gromacs/mdtypes/df_history.h"
  95 #include "gromacs/mdtypes/energyhistory.h"
  96 #include "gromacs/mdtypes/fcdata.h"
  97 #include "gromacs/mdtypes/forcerec.h"
  98 #include "gromacs/mdtypes/group.h"
  99 #include "gromacs/mdtypes/inputrec.h"
 100 #include "gromacs/mdtypes/interaction_const.h"
 101 #include "gromacs/mdtypes/md_enums.h"
 102 #include "gromacs/mdtypes/mdatom.h"
 103 #include "gromacs/mdtypes/observableshistory.h"
 104 #include "gromacs/mdtypes/state.h"
 105 #include "gromacs/pbcutil/mshift.h"
 106 #include "gromacs/pbcutil/pbc.h"
 107 #include "gromacs/pulling/pull.h"
 108 #include "gromacs/swap/swapcoords.h"
 109 #include "gromacs/timing/wallcycle.h"
 110 #include "gromacs/timing/walltime_accounting.h"
 111 #include "gromacs/topology/atoms.h"
 112 #include "gromacs/topology/idef.h"
 113 #include "gromacs/topology/mtop_util.h"
 114 #include "gromacs/topology/topology.h"
 115 #include "gromacs/trajectory/trajectoryframe.h"
 116 #include "gromacs/utility/basedefinitions.h"
 117 #include "gromacs/utility/cstringutil.h"
 118 #include "gromacs/utility/fatalerror.h"
 119 #include "gromacs/utility/logger.h"
 120 #include "gromacs/utility/real.h"
 121 #include "gromacs/utility/smalloc.h"
 122
 123 #include "deform.h"
 124 #include "membed.h"
 125 #include "repl_ex.h"
 126
 127 #ifdef GMX_FAHCORE
 128 #include "corewrap.h"
 129 #endif
 130
 131 using gmx::SimulationSignaller;
 132
 133 /*! \brief Check whether bonded interactions are missing, if appropriate
 134  *
 135  * \param[in]    fplog                                  Log file pointer
 136  * \param[in]    cr                                     Communication object
 137  * \param[in]    totalNumberOfBondedInteractions        Result of the global reduction over the number of bonds treated in each domain
 138  * \param[in]    top_global                             Global topology for the error message
 139  * \param[in]    top_local                              Local topology for the error message
 140  * \param[in]    state                                  Global state for the error message
 141  * \param[inout] shouldCheckNumberOfBondedInteractions  Whether we should do the check.
 142  *
 143  * \return Nothing, except that shouldCheckNumberOfBondedInteractions
 144  * is always set to false after exit.
 145  */
 146 static void checkNumberOfBondedInteractions(FILE *fplog, t_commrec *cr, int totalNumberOfBondedInteractions,
 147                                             gmx_mtop_t *top_global, gmx_localtop_t *top_local, t_state *state,
 148                                             bool *shouldCheckNumberOfBondedInteractions)
 149 {
 150     if (*shouldCheckNumberOfBondedInteractions)
 151     {
 152         if (totalNumberOfBondedInteractions != cr->dd->nbonded_global)
 153         {
 154             dd_print_missing_interactions(fplog, cr, totalNumberOfBondedInteractions, top_global, top_local, state); // Does not return
 155         }
 156         *shouldCheckNumberOfBondedInteractions = false;
 157     }
 158 }
 159
 160 static void reset_all_counters(FILE *fplog, const gmx::MDLogger &mdlog, t_commrec *cr,
 161                                gmx_int64_t step,
 162                                gmx_int64_t *step_rel, t_inputrec *ir,
 163                                gmx_wallcycle_t wcycle, t_nrnb *nrnb,
 164                                gmx_walltime_accounting_t walltime_accounting,
 165                                struct nonbonded_verlet_t *nbv)
 166 {
 167     char sbuf[STEPSTRSIZE];
 168
 169     /* Reset all the counters related to performance over the run */
 170     GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
 171             "step %s: resetting all time and cycle counters",
 172             gmx_step_str(step, sbuf));
 173
 174     if (use_GPU(nbv))
 175     {
 176         nbnxn_gpu_reset_timings(nbv);
 177         resetGpuProfiler();
 178     }
 179
 180     wallcycle_stop(wcycle, ewcRUN);
 181     wallcycle_reset_all(wcycle);
 182     if (DOMAINDECOMP(cr))
 183     {
 184         reset_dd_statistics_counters(cr->dd);
 185     }
 186     init_nrnb(nrnb);
 187     ir->init_step += *step_rel;
 188     ir->nsteps    -= *step_rel;
 189     *step_rel      = 0;
 190     wallcycle_start(wcycle, ewcRUN);
 191     walltime_accounting_start(walltime_accounting);
 192     print_date_and_time(fplog, cr->nodeid, "Restarted time", gmx_gettime());
 193 }
 194
 195 /*! \libinternal
 196     \copydoc integrator_t (FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
 197                            int nfile, const t_filenm fnm[],
 198                            const gmx_output_env_t *oenv, gmx_bool bVerbose,
 199                            int nstglobalcomm,
 200                            gmx_vsite_t *vsite, gmx_constr_t constr,
 201                            int stepout,
 202                            gmx::IMDOutputProvider *outputProvider,
 203                            t_inputrec *inputrec,
 204                            gmx_mtop_t *top_global, t_fcdata *fcd,
 205                            t_state *state_global,
 206                            t_mdatoms *mdatoms,
 207                            t_nrnb *nrnb, gmx_wallcycle_t wcycle,
 208                            t_forcerec *fr,
 209                            int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
 210                            real cpt_period, real max_hours,
 211                            int imdport,
 212                            unsigned long Flags,
 213                            gmx_walltime_accounting_t walltime_accounting)
 214  */
 215 double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
 216                   int nfile, const t_filenm fnm[],
 217                   const gmx_output_env_t *oenv, gmx_bool bVerbose,
 218                   int nstglobalcomm,
 219                   gmx_vsite_t *vsite, gmx_constr_t constr,
 220                   int stepout, gmx::IMDOutputProvider *outputProvider,
 221                   t_inputrec *ir,
 222                   gmx_mtop_t *top_global,
 223                   t_fcdata *fcd,
 224                   t_state *state_global,
 225                   ObservablesHistory *observablesHistory,
 226                   t_mdatoms *mdatoms,
 227                   t_nrnb *nrnb, gmx_wallcycle_t wcycle,
 228                   t_forcerec *fr,
 229                   const ReplicaExchangeParameters &replExParams,
 230                   gmx_membed_t *membed,
 231                   real cpt_period, real max_hours,
 232                   int imdport,
 233                   unsigned long Flags,
 234                   gmx_walltime_accounting_t walltime_accounting)
 235 {
 236     gmx_mdoutf_t    outf = nullptr;
 237     gmx_int64_t     step, step_rel;
 238     double          elapsed_time;
 239     double          t, t0, lam0[efptNR];
 240     gmx_bool        bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
 241     gmx_bool        bNS, bNStList, bSimAnn, bStopCM, bRerunMD,
 242                     bFirstStep, startingFromCheckpoint, bInitStep, bLastStep = FALSE,
 243                     bBornRadii, bUsingEnsembleRestraints;
 244     gmx_bool          bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
 245     gmx_bool          do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE,
 246                       bForceUpdate = FALSE, bCPT;
 247     gmx_bool          bMasterState;
 248     int               force_flags, cglo_flags;
 249     tensor            force_vir, shake_vir, total_vir, tmp_vir, pres;
 250     int               i, m;
 251     t_trxstatus      *status;
 252     rvec              mu_tot;
 253     t_vcm            *vcm;
 254     matrix            parrinellorahmanMu, M;
 255     t_trxframe        rerun_fr;
 256     gmx_repl_ex_t     repl_ex = nullptr;
 257     int               nchkpt  = 1;
 258     gmx_localtop_t   *top;
 259     t_mdebin         *mdebin   = nullptr;
 260     gmx_enerdata_t   *enerd;
 261     PaddedRVecVector  f {};
 262     gmx_global_stat_t gstat;
 263     gmx_update_t     *upd   = nullptr;
 264     t_graph          *graph = nullptr;
 265     gmx_groups_t     *groups;
 266     gmx_ekindata_t   *ekind;
 267     gmx_shellfc_t    *shellfc;
 268     gmx_bool          bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
 269     gmx_bool          bResetCountersHalfMaxH = FALSE;
 270     gmx_bool          bTemp, bPres, bTrotter;
 271     real              dvdl_constr;
 272     rvec             *cbuf        = nullptr;
 273     int               cbuf_nalloc = 0;
 274     matrix            lastbox;
 275     int               lamnew  = 0;
 276     /* for FEP */
 277     int               nstfep = 0;
 278     double            cycles;
 279     real              saved_conserved_quantity = 0;
 280     real              last_ekin                = 0;
 281     t_extmass         MassQ;
 282     int             **trotter_seq;
 283     char              sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
 284     int               handled_stop_condition = gmx_stop_cond_none; /* compare to get_stop_condition*/
 285
 286
 287     /* PME load balancing data for GPU kernels */
 288     pme_load_balancing_t *pme_loadbal      = nullptr;
 289     gmx_bool              bPMETune         = FALSE;
 290     gmx_bool              bPMETunePrinting = FALSE;
 291
 292     /* Interactive MD */
 293     gmx_bool          bIMDstep = FALSE;
 294
 295 #ifdef GMX_FAHCORE
 296     /* Temporary addition for FAHCORE checkpointing */
 297     int chkpt_ret;
 298 #endif
 299     /* Domain decomposition could incorrectly miss a bonded
 300        interaction, but checking for that requires a global
 301        communication stage, which does not otherwise happen in DD
 302        code. So we do that alongside the first global energy reduction
 303        after a new DD is made. These variables handle whether the
 304        check happens, and the result it returns. */
 305     bool              shouldCheckNumberOfBondedInteractions = false;
 306     int               totalNumberOfBondedInteractions       = -1;
 307
 308     SimulationSignals signals;
 309     // Most global communnication stages don't propagate mdrun
 310     // signals, and will use this object to achieve that.
 311     SimulationSignaller nullSignaller(nullptr, nullptr, false, false);
 312
 313     /* Check for special mdrun options */
 314     bRerunMD = (Flags & MD_RERUN);
 315     if (Flags & MD_RESETCOUNTERSHALFWAY)
 316     {
 317         if (ir->nsteps > 0)
 318         {
 319             /* Signal to reset the counters half the simulation steps. */
 320             wcycle_set_reset_counters(wcycle, ir->nsteps/2);
 321         }
 322         /* Signal to reset the counters halfway the simulation time. */
 323         bResetCountersHalfMaxH = (max_hours > 0);
 324     }
 325
 326     /* md-vv uses averaged full step velocities for T-control
 327        md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
 328        md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
 329     bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
 330
 331     if (bRerunMD)
 332     {
 333         /* Since we don't know if the frames read are related in any way,
 334          * rebuild the neighborlist at every step.
 335          */
 336         ir->nstlist       = 1;
 337         ir->nstcalcenergy = 1;
 338         nstglobalcomm     = 1;
 339     }
 340
 341     nstglobalcomm   = check_nstglobalcomm(mdlog, nstglobalcomm, ir);
 342     bGStatEveryStep = (nstglobalcomm == 1);
 343
 344     if (bRerunMD)
 345     {
 346         ir->nstxout_compressed = 0;
 347     }
 348     groups = &top_global->groups;
 349
 350     gmx_edsam *ed = nullptr;
 351     if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
 352     {
 353         /* Initialize essential dynamics sampling */
 354         ed = init_edsam(opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
 355                         top_global, ir, cr, constr,
 356                         as_rvec_array(state_global->x.data()),
 357                         state_global->box, observablesHistory,
 358                         oenv, Flags & MD_APPENDFILES);
 359     }
 360
 361     if (ir->eSwapCoords != eswapNO)
 362     {
 363         /* Initialize ion swapping code */
 364         init_swapcoords(fplog, bVerbose, ir, opt2fn_master("-swap", nfile, fnm, cr),
 365                         top_global, as_rvec_array(state_global->x.data()), state_global->box, observablesHistory, cr, oenv, Flags);
 366     }
 367
 368     /* Initial values */
 369     init_md(fplog, cr, outputProvider, ir, oenv, &t, &t0, state_global->lambda,
 370             &(state_global->fep_state), lam0,
 371             nrnb, top_global, &upd,
 372             nfile, fnm, &outf, &mdebin,
 373             force_vir, shake_vir, mu_tot, &bSimAnn, &vcm, Flags, wcycle);
 374
 375     clear_mat(total_vir);
 376     clear_mat(pres);
 377     /* Energy terms and groups */
 378     snew(enerd, 1);
 379     init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
 380                   enerd);
 381
 382     /* Kinetic energy data */
 383     snew(ekind, 1);
 384     init_ekindata(fplog, top_global, &(ir->opts), ekind);
 385     /* Copy the cos acceleration to the groups struct */
 386     ekind->cosacc.cos_accel = ir->cos_accel;
 387
 388     gstat = global_stat_init(ir);
 389
 390     /* Check for polarizable models and flexible constraints */
 391     shellfc = init_shell_flexcon(fplog,
 392                                  top_global, n_flexible_constraints(constr),
 393                                  ir->nstcalcenergy, DOMAINDECOMP(cr));
 394
 395     if (shellfc && ir->nstcalcenergy != 1)
 396     {
 397         gmx_fatal(FARGS, "You have nstcalcenergy set to a value (%d) that is different from 1.\nThis is not supported in combinations with shell particles.\nPlease make a new tpr file.", ir->nstcalcenergy);
 398     }
 399     if (shellfc && DOMAINDECOMP(cr))
 400     {
 401         gmx_fatal(FARGS, "Shell particles are not implemented with domain decomposition, use a single rank");
 402     }
 403
 404     if (inputrecDeform(ir))
 405     {
 406         tMPI_Thread_mutex_lock(&deform_init_box_mutex);
 407         set_deform_reference_box(upd,
 408                                  deform_init_init_step_tpx,
 409                                  deform_init_box_tpx);
 410         tMPI_Thread_mutex_unlock(&deform_init_box_mutex);
 411     }
 412
 413     {
 414         double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1);
 415         if ((io > 2000) && MASTER(cr))
 416         {
 417             fprintf(stderr,
 418                     "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
 419                     io);
 420         }
 421     }
 422
 423     std::unique_ptr<t_state> stateInstance;
 424     t_state *                state;
 425
 426     if (DOMAINDECOMP(cr))
 427     {
 428         top = dd_init_local_top(top_global);
 429
 430         stateInstance = std::unique_ptr<t_state>(new t_state);
 431         state         = stateInstance.get();
 432         dd_init_local_state(cr->dd, state_global, state);
 433     }
 434     else
 435     {
 436         state_change_natoms(state_global, state_global->natoms);
 437         /* We need to allocate one element extra, since we might use
 438          * (unaligned) 4-wide SIMD loads to access rvec entries.
 439          */
 440         f.resize(state_global->natoms + 1);
 441         /* Copy the pointer to the global state */
 442         state = state_global;
 443
 444         snew(top, 1);
 445         mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
 446                                   &graph, mdatoms, vsite, shellfc);
 447
 448         update_realloc(upd, state->natoms);
 449     }
 450
 451     /* Set up interactive MD (IMD) */
 452     init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, as_rvec_array(state_global->x.data()),
 453              nfile, fnm, oenv, imdport, Flags);
 454
 455     if (DOMAINDECOMP(cr))
 456     {
 457         /* Distribute the charge groups over the nodes from the master node */
 458         dd_partition_system(fplog, ir->init_step, cr, TRUE, 1,
 459                             state_global, top_global, ir,
 460                             state, &f, mdatoms, top, fr,
 461                             vsite, constr,
 462                             nrnb, nullptr, FALSE);
 463         shouldCheckNumberOfBondedInteractions = true;
 464         update_realloc(upd, state->natoms);
 465     }
 466
 467     update_mdatoms(mdatoms, state->lambda[efptMASS]);
 468
 469     startingFromCheckpoint = Flags & MD_STARTFROMCPT;
 470
 471     if (ir->bExpanded)
 472     {
 473         init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
 474     }
 475
 476     if (MASTER(cr))
 477     {
 478         if (startingFromCheckpoint)
 479         {
 480             /* Update mdebin with energy history if appending to output files */
 481             if (Flags & MD_APPENDFILES)
 482             {
 483                 restore_energyhistory_from_state(mdebin, observablesHistory->energyHistory.get());
 484             }
 485             else if (observablesHistory->energyHistory.get() != nullptr)
 486             {
 487                 /* We might have read an energy history from checkpoint.
 488                  * As we are not appending, we want to restart the statistics.
 489                  * Free the allocated memory and reset the counts.
 490                  */
 491                 observablesHistory->energyHistory = {};
 492             }
 493         }
 494         if (observablesHistory->energyHistory.get() == nullptr)
 495         {
 496             observablesHistory->energyHistory = std::unique_ptr<energyhistory_t>(new energyhistory_t {});
 497         }
 498         /* Set the initial energy history in state by updating once */
 499         update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
 500     }
 501
 502     /* Initialize constraints */
 503     if (constr && !DOMAINDECOMP(cr))
 504     {
 505         set_constraints(constr, top, ir, mdatoms, cr);
 506     }
 507
 508     const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
 509     if (useReplicaExchange && MASTER(cr))
 510     {
 511         repl_ex = init_replica_exchange(fplog, cr->ms, state_global, ir,
 512                                         replExParams);
 513     }
 514
 515     /* PME tuning is only supported with PME for Coulomb. Is is not supported
 516      * with only LJ PME, or for reruns.
 517      */
 518     bPMETune = ((Flags & MD_TUNEPME) && EEL_PME(fr->eeltype) && !bRerunMD &&
 519                 !(Flags & MD_REPRODUCIBLE));
 520     if (bPMETune)
 521     {
 522         pme_loadbal_init(&pme_loadbal, cr, mdlog, ir, state->box,
 523                          fr->ic, fr->nbv->listParams.get(), fr->pmedata, use_GPU(fr->nbv),
 524                          &bPMETunePrinting);
 525     }
 526
 527     if (!ir->bContinuation && !bRerunMD)
 528     {
 529         if (state->flags & (1 << estV))
 530         {
 531             /* Set the velocities of vsites, shells and frozen atoms to zero */
 532             for (i = 0; i < mdatoms->homenr; i++)
 533             {
 534                 if (mdatoms->ptype[i] == eptVSite ||
 535                     mdatoms->ptype[i] == eptShell)
 536                 {
 537                     clear_rvec(state->v[i]);
 538                 }
 539                 else if (mdatoms->cFREEZE)
 540                 {
 541                     for (m = 0; m < DIM; m++)
 542                     {
 543                         if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 544                         {
 545                             state->v[i][m] = 0;
 546                         }
 547                     }
 548                 }
 549             }
 550         }
 551
 552         if (constr)
 553         {
 554             /* Constrain the initial coordinates and velocities */
 555             do_constrain_first(fplog, constr, ir, mdatoms, state,
 556                                cr, nrnb, fr, top);
 557         }
 558         if (vsite)
 559         {
 560             /* Construct the virtual sites for the initial configuration */
 561             construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, nullptr,
 562                              top->idef.iparams, top->idef.il,
 563                              fr->ePBC, fr->bMolPBC, cr, state->box);
 564         }
 565     }
 566
 567     if (ir->efep != efepNO)
 568     {
 569         /* Set free energy calculation frequency as the greatest common
 570          * denominator of nstdhdl and repl_ex_nst. */
 571         nstfep = ir->fepvals->nstdhdl;
 572         if (ir->bExpanded)
 573         {
 574             nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
 575         }
 576         if (useReplicaExchange)
 577         {
 578             nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
 579         }
 580     }
 581
 582     /* Be REALLY careful about what flags you set here. You CANNOT assume
 583      * this is the first step, since we might be restarting from a checkpoint,
 584      * and in that case we should not do any modifications to the state.
 585      */
 586     bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
 587
 588     if (Flags & MD_READ_EKIN)
 589     {
 590         restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 591     }
 592
 593     cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
 594                   | (bStopCM ? CGLO_STOPCM : 0)
 595                   | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
 596                   | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
 597                   | ((Flags & MD_READ_EKIN) ? CGLO_READEKIN : 0));
 598
 599     bSumEkinhOld = FALSE;
 600     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 601                     nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 602                     constr, &nullSignaller, state->box,
 603                     &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags
 604                     | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
 605     checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
 606                                     top_global, top, state,
 607                                     &shouldCheckNumberOfBondedInteractions);
 608     if (ir->eI == eiVVAK)
 609     {
 610         /* a second call to get the half step temperature initialized as well */
 611         /* we do the same call as above, but turn the pressure off -- internally to
 612            compute_globals, this is recognized as a velocity verlet half-step
 613            kinetic energy calculation.  This minimized excess variables, but
 614            perhaps loses some logic?*/
 615
 616         compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 617                         nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 618                         constr, &nullSignaller, state->box,
 619                         nullptr, &bSumEkinhOld,
 620                         cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
 621     }
 622
 623     /* Calculate the initial half step temperature, and save the ekinh_old */
 624     if (!(Flags & MD_STARTFROMCPT))
 625     {
 626         for (i = 0; (i < ir->opts.ngtc); i++)
 627         {
 628             copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
 629         }
 630     }
 631
 632     /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
 633        temperature control */
 634     trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
 635
 636     if (MASTER(cr))
 637     {
 638         if (!ir->bContinuation)
 639         {
 640             if (constr && ir->eConstrAlg == econtLINCS)
 641             {
 642                 fprintf(fplog,
 643                         "RMS relative constraint deviation after constraining: %.2e\n",
 644                         constr_rmsd(constr));
 645             }
 646             if (EI_STATE_VELOCITY(ir->eI))
 647             {
 648                 real temp = enerd->term[F_TEMP];
 649                 if (ir->eI != eiVV)
 650                 {
 651                     /* Result of Ekin averaged over velocities of -half
 652                      * and +half step, while we only have -half step here.
 653                      */
 654                     temp *= 2;
 655                 }
 656                 fprintf(fplog, "Initial temperature: %g K\n", temp);
 657             }
 658         }
 659
 660         if (bRerunMD)
 661         {
 662             fprintf(stderr, "starting md rerun '%s', reading coordinates from"
 663                     " input trajectory '%s'\n\n",
 664                     *(top_global->name), opt2fn("-rerun", nfile, fnm));
 665             if (bVerbose)
 666             {
 667                 fprintf(stderr, "Calculated time to finish depends on nsteps from "
 668                         "run input file,\nwhich may not correspond to the time "
 669                         "needed to process input trajectory.\n\n");
 670             }
 671         }
 672         else
 673         {
 674             char tbuf[20];
 675             fprintf(stderr, "starting mdrun '%s'\n",
 676                     *(top_global->name));
 677             if (ir->nsteps >= 0)
 678             {
 679                 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
 680             }
 681             else
 682             {
 683                 sprintf(tbuf, "%s", "infinite");
 684             }
 685             if (ir->init_step > 0)
 686             {
 687                 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 688                         gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
 689                         gmx_step_str(ir->init_step, sbuf2),
 690                         ir->init_step*ir->delta_t);
 691             }
 692             else
 693             {
 694                 fprintf(stderr, "%s steps, %s ps.\n",
 695                         gmx_step_str(ir->nsteps, sbuf), tbuf);
 696             }
 697         }
 698         fprintf(fplog, "\n");
 699     }
 700
 701     walltime_accounting_start(walltime_accounting);
 702     wallcycle_start(wcycle, ewcRUN);
 703     print_start(fplog, cr, walltime_accounting, "mdrun");
 704
 705     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
 706 #ifdef GMX_FAHCORE
 707     chkpt_ret = fcCheckPointParallel( cr->nodeid,
 708                                       NULL, 0);
 709     if (chkpt_ret == 0)
 710     {
 711         gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
 712     }
 713 #endif
 714
 715     /***********************************************************
 716      *
 717      *             Loop over MD steps
 718      *
 719      ************************************************************/
 720
 721     /* if rerunMD then read coordinates and velocities from input trajectory */
 722     if (bRerunMD)
 723     {
 724         if (getenv("GMX_FORCE_UPDATE"))
 725         {
 726             bForceUpdate = TRUE;
 727         }
 728
 729         rerun_fr.natoms = 0;
 730         if (MASTER(cr))
 731         {
 732             bLastStep = !read_first_frame(oenv, &status,
 733                                           opt2fn("-rerun", nfile, fnm),
 734                                           &rerun_fr, TRX_NEED_X | TRX_READ_V);
 735             if (rerun_fr.natoms != top_global->natoms)
 736             {
 737                 gmx_fatal(FARGS,
 738                           "Number of atoms in trajectory (%d) does not match the "
 739                           "run input file (%d)\n",
 740                           rerun_fr.natoms, top_global->natoms);
 741             }
 742             if (ir->ePBC != epbcNONE)
 743             {
 744                 if (!rerun_fr.bBox)
 745                 {
 746                     gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f does not contain a box, while pbc is used", rerun_fr.step, rerun_fr.time);
 747                 }
 748                 if (max_cutoff2(ir->ePBC, rerun_fr.box) < gmx::square(fr->rlist))
 749                 {
 750                     gmx_fatal(FARGS, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr.step, rerun_fr.time);
 751                 }
 752             }
 753         }
 754
 755         if (PAR(cr))
 756         {
 757             rerun_parallel_comm(cr, &rerun_fr, &bLastStep);
 758         }
 759
 760         if (ir->ePBC != epbcNONE)
 761         {
 762             /* Set the shift vectors.
 763              * Necessary here when have a static box different from the tpr box.
 764              */
 765             calc_shifts(rerun_fr.box, fr->shift_vec);
 766         }
 767     }
 768
 769     /* loop over MD steps or if rerunMD to end of input trajectory */
 770     bFirstStep = TRUE;
 771     /* Skip the first Nose-Hoover integration when we get the state from tpx */
 772     bInitStep        = !startingFromCheckpoint || EI_VV(ir->eI);
 773     bSumEkinhOld     = FALSE;
 774     bExchanged       = FALSE;
 775     bNeedRepartition = FALSE;
 776     // TODO This implementation of ensemble orientation restraints is nasty because
 777     // a user can't just do multi-sim with single-sim orientation restraints.
 778     bUsingEnsembleRestraints = (fcd->disres.nsystems > 1) || (cr->ms && fcd->orires.nr);
 779
 780     {
 781         // Replica exchange and ensemble restraints need all
 782         // simulations to remain synchronized, so they need
 783         // checkpoints and stop conditions to act on the same step, so
 784         // the propagation of such signals must take place between
 785         // simulations, not just within simulations.
 786         bool checkpointIsLocal    = !useReplicaExchange && !bUsingEnsembleRestraints;
 787         bool stopConditionIsLocal = !useReplicaExchange && !bUsingEnsembleRestraints;
 788         bool resetCountersIsLocal = true;
 789         signals[eglsCHKPT]         = SimulationSignal(checkpointIsLocal);
 790         signals[eglsSTOPCOND]      = SimulationSignal(stopConditionIsLocal);
 791         signals[eglsRESETCOUNTERS] = SimulationSignal(resetCountersIsLocal);
 792     }
 793
 794     DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion   = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
 795     DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion  = (DOMAINDECOMP(cr) ? DdCloseBalanceRegionAfterForceComputation::yes : DdCloseBalanceRegionAfterForceComputation::no);
 796
 797     step     = ir->init_step;
 798     step_rel = 0;
 799
 800     // TODO extract this to new multi-simulation module
 801     if (MASTER(cr) && MULTISIM(cr) && !useReplicaExchange)
 802     {
 803         if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
 804         {
 805             GMX_LOG(mdlog.warning).appendText(
 806                     "Note: The number of steps is not consistent across multi simulations,\n"
 807                     "but we are proceeding anyway!");
 808         }
 809         if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
 810         {
 811             GMX_LOG(mdlog.warning).appendText(
 812                     "Note: The initial step is not consistent across multi simulations,\n"
 813                     "but we are proceeding anyway!");
 814         }
 815     }
 816
 817     /* and stop now if we should */
 818     bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
 819     while (!bLastStep)
 820     {
 821
 822         /* Determine if this is a neighbor search step */
 823         bNStList = (ir->nstlist > 0  && step % ir->nstlist == 0);
 824
 825         if (bPMETune && bNStList)
 826         {
 827             /* PME grid + cut-off optimization with GPUs or PME nodes */
 828             pme_loadbal_do(pme_loadbal, cr,
 829                            (bVerbose && MASTER(cr)) ? stderr : nullptr,
 830                            fplog, mdlog,
 831                            ir, fr, state,
 832                            wcycle,
 833                            step, step_rel,
 834                            &bPMETunePrinting);
 835         }
 836
 837         wallcycle_start(wcycle, ewcSTEP);
 838
 839         if (bRerunMD)
 840         {
 841             if (rerun_fr.bStep)
 842             {
 843                 step     = rerun_fr.step;
 844                 step_rel = step - ir->init_step;
 845             }
 846             if (rerun_fr.bTime)
 847             {
 848                 t = rerun_fr.time;
 849             }
 850             else
 851             {
 852                 t = step;
 853             }
 854         }
 855         else
 856         {
 857             bLastStep = (step_rel == ir->nsteps);
 858             t         = t0 + step*ir->delta_t;
 859         }
 860
 861         // TODO Refactor this, so that nstfep does not need a default value of zero
 862         if (ir->efep != efepNO || ir->bSimTemp)
 863         {
 864             /* find and set the current lambdas.  If rerunning, we either read in a state, or a lambda value,
 865                requiring different logic. */
 866
 867             set_current_lambdas(step, ir->fepvals, bRerunMD, &rerun_fr, state_global, state, lam0);
 868             bDoDHDL      = do_per_step(step, ir->fepvals->nstdhdl);
 869             bDoFEP       = ((ir->efep != efepNO) && do_per_step(step, nstfep));
 870             bDoExpanded  = (do_per_step(step, ir->expandedvals->nstexpanded)
 871                             && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
 872         }
 873
 874         bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
 875                      do_per_step(step, replExParams.exchangeInterval));
 876
 877         if (bSimAnn)
 878         {
 879             update_annealing_target_temp(ir, t, upd);
 880         }
 881
 882         if (bRerunMD)
 883         {
 884             if (!DOMAINDECOMP(cr) || MASTER(cr))
 885             {
 886                 for (i = 0; i < state_global->natoms; i++)
 887                 {
 888                     copy_rvec(rerun_fr.x[i], state_global->x[i]);
 889                 }
 890                 if (rerun_fr.bV)
 891                 {
 892                     for (i = 0; i < state_global->natoms; i++)
 893                     {
 894                         copy_rvec(rerun_fr.v[i], state_global->v[i]);
 895                     }
 896                 }
 897                 else
 898                 {
 899                     for (i = 0; i < state_global->natoms; i++)
 900                     {
 901                         clear_rvec(state_global->v[i]);
 902                     }
 903                     if (bRerunWarnNoV)
 904                     {
 905                         fprintf(stderr, "\nWARNING: Some frames do not contain velocities.\n"
 906                                 "         Ekin, temperature and pressure are incorrect,\n"
 907                                 "         the virial will be incorrect when constraints are present.\n"
 908                                 "\n");
 909                         bRerunWarnNoV = FALSE;
 910                     }
 911                 }
 912             }
 913             copy_mat(rerun_fr.box, state_global->box);
 914             copy_mat(state_global->box, state->box);
 915
 916             if (vsite && (Flags & MD_RERUN_VSITE))
 917             {
 918                 if (DOMAINDECOMP(cr))
 919                 {
 920                     gmx_fatal(FARGS, "Vsite recalculation with -rerun is not implemented with domain decomposition, use a single rank");
 921                 }
 922                 if (graph)
 923                 {
 924                     /* Following is necessary because the graph may get out of sync
 925                      * with the coordinates if we only have every N'th coordinate set
 926                      */
 927                     mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
 928                     shift_self(graph, state->box, as_rvec_array(state->x.data()));
 929                 }
 930                 construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
 931                                  top->idef.iparams, top->idef.il,
 932                                  fr->ePBC, fr->bMolPBC, cr, state->box);
 933                 if (graph)
 934                 {
 935                     unshift_self(graph, state->box, as_rvec_array(state->x.data()));
 936                 }
 937             }
 938         }
 939
 940         /* Stop Center of Mass motion */
 941         bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
 942
 943         if (bRerunMD)
 944         {
 945             /* for rerun MD always do Neighbour Searching */
 946             bNS      = (bFirstStep || ir->nstlist != 0);
 947         }
 948         else
 949         {
 950             /* Determine whether or not to do Neighbour Searching */
 951             bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
 952         }
 953
 954         /* < 0 means stop at next step, > 0 means stop at next NS step */
 955         if ( (signals[eglsSTOPCOND].set < 0) ||
 956              ( (signals[eglsSTOPCOND].set > 0 ) && ( bNS || ir->nstlist == 0)))
 957         {
 958             bLastStep = TRUE;
 959         }
 960
 961         /* Determine whether or not to update the Born radii if doing GB */
 962         bBornRadii = bFirstStep;
 963         if (ir->implicit_solvent && (step % ir->nstgbradii == 0))
 964         {
 965             bBornRadii = TRUE;
 966         }
 967
 968         /* do_log triggers energy and virial calculation. Because this leads
 969          * to different code paths, forces can be different. Thus for exact
 970          * continuation we should avoid extra log output.
 971          * Note that the || bLastStep can result in non-exact continuation
 972          * beyond the last step. But we don't consider that to be an issue.
 973          */
 974         do_log     = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep || bRerunMD;
 975         do_verbose = bVerbose &&
 976             (step % stepout == 0 || bFirstStep || bLastStep || bRerunMD);
 977
 978         if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
 979         {
 980             if (bRerunMD)
 981             {
 982                 bMasterState = TRUE;
 983             }
 984             else
 985             {
 986                 bMasterState = FALSE;
 987                 /* Correct the new box if it is too skewed */
 988                 if (inputrecDynamicBox(ir))
 989                 {
 990                     if (correct_box(fplog, step, state->box, graph))
 991                     {
 992                         bMasterState = TRUE;
 993                     }
 994                 }
 995                 if (DOMAINDECOMP(cr) && bMasterState)
 996                 {
 997                     dd_collect_state(cr->dd, state, state_global);
 998                 }
 999             }
1000
1001             if (DOMAINDECOMP(cr))
1002             {
1003                 /* Repartition the domain decomposition */
1004                 dd_partition_system(fplog, step, cr,
1005                                     bMasterState, nstglobalcomm,
1006                                     state_global, top_global, ir,
1007                                     state, &f, mdatoms, top, fr,
1008                                     vsite, constr,
1009                                     nrnb, wcycle,
1010                                     do_verbose && !bPMETunePrinting);
1011                 shouldCheckNumberOfBondedInteractions = true;
1012                 update_realloc(upd, state->natoms);
1013             }
1014         }
1015
1016         if (MASTER(cr) && do_log)
1017         {
1018             print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
1019         }
1020
1021         if (ir->efep != efepNO)
1022         {
1023             update_mdatoms(mdatoms, state->lambda[efptMASS]);
1024         }
1025
1026         if ((bRerunMD && rerun_fr.bV) || bExchanged)
1027         {
1028
1029             /* We need the kinetic energy at minus the half step for determining
1030              * the full step kinetic energy and possibly for T-coupling.*/
1031             /* This may not be quite working correctly yet . . . . */
1032             compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1033                             wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1034                             constr, &nullSignaller, state->box,
1035                             &totalNumberOfBondedInteractions, &bSumEkinhOld,
1036                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
1037             checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
1038                                             top_global, top, state,
1039                                             &shouldCheckNumberOfBondedInteractions);
1040         }
1041         clear_mat(force_vir);
1042
1043         /* We write a checkpoint at this MD step when:
1044          * either at an NS step when we signalled through gs,
1045          * or at the last step (but not when we do not want confout),
1046          * but never at the first step or with rerun.
1047          */
1048         bCPT = (((signals[eglsCHKPT].set && (bNS || ir->nstlist == 0)) ||
1049                  (bLastStep && (Flags & MD_CONFOUT))) &&
1050                 step > ir->init_step && !bRerunMD);
1051         if (bCPT)
1052         {
1053             signals[eglsCHKPT].set = 0;
1054         }
1055
1056         /* Determine the energy and pressure:
1057          * at nstcalcenergy steps and at energy output steps (set below).
1058          */
1059         if (EI_VV(ir->eI) && (!bInitStep))
1060         {
1061             /* for vv, the first half of the integration actually corresponds
1062                to the previous step.  bCalcEner is only required to be evaluated on the 'next' step,
1063                but the virial needs to be calculated on both the current step and the 'next' step. Future
1064                reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
1065
1066             /* TODO: This is probably not what we want, we will write to energy file one step after nstcalcenergy steps. */
1067             bCalcEnerStep = do_per_step(step - 1, ir->nstcalcenergy);
1068             bCalcVir      = bCalcEnerStep ||
1069                 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
1070         }
1071         else
1072         {
1073             bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
1074             bCalcVir      = bCalcEnerStep ||
1075                 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
1076         }
1077         bCalcEner = bCalcEnerStep;
1078
1079         do_ene = (do_per_step(step, ir->nstenergy) || bLastStep || bRerunMD);
1080
1081         if (do_ene || do_log || bDoReplEx)
1082         {
1083             bCalcVir  = TRUE;
1084             bCalcEner = TRUE;
1085         }
1086
1087         /* Do we need global communication ? */
1088         bGStat = (bCalcVir || bCalcEner || bStopCM ||
1089                   do_per_step(step, nstglobalcomm) ||
1090                   (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
1091
1092         force_flags = (GMX_FORCE_STATECHANGED |
1093                        ((inputrecDynamicBox(ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1094                        GMX_FORCE_ALLFORCES |
1095                        (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
1096                        (bCalcEner ? GMX_FORCE_ENERGY : 0) |
1097                        (bDoFEP ? GMX_FORCE_DHDL : 0)
1098                        );
1099
1100         if (shellfc)
1101         {
1102             /* Now is the time to relax the shells */
1103             relax_shell_flexcon(fplog, cr, bVerbose, step,
1104                                 ir, bNS, force_flags, top,
1105                                 constr, enerd, fcd,
1106                                 state, &f, force_vir, mdatoms,
1107                                 nrnb, wcycle, graph, groups,
1108                                 shellfc, fr, bBornRadii, t, mu_tot,
1109                                 vsite,
1110                                 ddOpenBalanceRegion, ddCloseBalanceRegion);
1111         }
1112         else
1113         {
1114             /* The coordinates (x) are shifted (to get whole molecules)
1115              * in do_force.
1116              * This is parallellized as well, and does communication too.
1117              * Check comments in sim_util.c
1118              */
1119             do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
1120                      state->box, &state->x, &state->hist,
1121                      &f, force_vir, mdatoms, enerd, fcd,
1122                      state->lambda, graph,
1123                      fr, vsite, mu_tot, t, ed, bBornRadii,
1124                      (bNS ? GMX_FORCE_NS : 0) | force_flags,
1125                      ddOpenBalanceRegion, ddCloseBalanceRegion);
1126         }
1127
1128         if (EI_VV(ir->eI) && !startingFromCheckpoint && !bRerunMD)
1129         /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1130         {
1131             rvec *vbuf = nullptr;
1132
1133             wallcycle_start(wcycle, ewcUPDATE);
1134             if (ir->eI == eiVV && bInitStep)
1135             {
1136                 /* if using velocity verlet with full time step Ekin,
1137                  * take the first half step only to compute the
1138                  * virial for the first step. From there,
1139                  * revert back to the initial coordinates
1140                  * so that the input is actually the initial step.
1141                  */
1142                 snew(vbuf, state->natoms);
1143                 copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
1144             }
1145             else
1146             {
1147                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1148                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
1149             }
1150
1151             update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
1152                           ekind, M, upd, etrtVELOCITY1,
1153                           cr, constr);
1154
1155             if (!bRerunMD || rerun_fr.bV || bForceUpdate)         /* Why is rerun_fr.bV here?  Unclear. */
1156             {
1157                 wallcycle_stop(wcycle, ewcUPDATE);
1158                 update_constraints(fplog, step, nullptr, ir, mdatoms,
1159                                    state, fr->bMolPBC, graph, &f,
1160                                    &top->idef, shake_vir,
1161                                    cr, nrnb, wcycle, upd, constr,
1162                                    TRUE, bCalcVir);
1163                 wallcycle_start(wcycle, ewcUPDATE);
1164             }
1165             else if (graph)
1166             {
1167                 /* Need to unshift here if a do_force has been
1168                    called in the previous step */
1169                 unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1170             }
1171             /* if VV, compute the pressure and constraints */
1172             /* For VV2, we strictly only need this if using pressure
1173              * control, but we really would like to have accurate pressures
1174              * printed out.
1175              * Think about ways around this in the future?
1176              * For now, keep this choice in comments.
1177              */
1178             /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
1179             /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
1180             bPres = TRUE;
1181             bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
1182             if (bCalcEner && ir->eI == eiVVAK)
1183             {
1184                 bSumEkinhOld = TRUE;
1185             }
1186             /* for vv, the first half of the integration actually corresponds to the previous step.
1187                So we need information from the last step in the first half of the integration */
1188             if (bGStat || do_per_step(step-1, nstglobalcomm))
1189             {
1190                 wallcycle_stop(wcycle, ewcUPDATE);
1191                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1192                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1193                                 constr, &nullSignaller, state->box,
1194                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1195                                 (bGStat ? CGLO_GSTAT : 0)
1196                                 | CGLO_ENERGY
1197                                 | (bTemp ? CGLO_TEMPERATURE : 0)
1198                                 | (bPres ? CGLO_PRESSURE : 0)
1199                                 | (bPres ? CGLO_CONSTRAINT : 0)
1200                                 | (bStopCM ? CGLO_STOPCM : 0)
1201                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1202                                 | CGLO_SCALEEKIN
1203                                 );
1204                 /* explanation of above:
1205                    a) We compute Ekin at the full time step
1206                    if 1) we are using the AveVel Ekin, and it's not the
1207                    initial step, or 2) if we are using AveEkin, but need the full
1208                    time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1209                    b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1210                    EkinAveVel because it's needed for the pressure */
1211                 checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
1212                                                 top_global, top, state,
1213                                                 &shouldCheckNumberOfBondedInteractions);
1214                 wallcycle_start(wcycle, ewcUPDATE);
1215             }
1216             /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1217             if (!bInitStep)
1218             {
1219                 if (bTrotter)
1220                 {
1221                     m_add(force_vir, shake_vir, total_vir);     /* we need the un-dispersion corrected total vir here */
1222                     trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
1223
1224                     /* TODO This is only needed when we're about to write
1225                      * a checkpoint, because we use it after the restart
1226                      * (in a kludge?). But what should we be doing if
1227                      * startingFromCheckpoint or bInitStep are true? */
1228                     if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1229                     {
1230                         copy_mat(shake_vir, state->svir_prev);
1231                         copy_mat(force_vir, state->fvir_prev);
1232                     }
1233                     if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
1234                     {
1235                         /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1236                         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
1237                         enerd->term[F_EKIN] = trace(ekind->ekin);
1238                     }
1239                 }
1240                 else if (bExchanged)
1241                 {
1242                     wallcycle_stop(wcycle, ewcUPDATE);
1243                     /* We need the kinetic energy at minus the half step for determining
1244                      * the full step kinetic energy and possibly for T-coupling.*/
1245                     /* This may not be quite working correctly yet . . . . */
1246                     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1247                                     wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1248                                     constr, &nullSignaller, state->box,
1249                                     nullptr, &bSumEkinhOld,
1250                                     CGLO_GSTAT | CGLO_TEMPERATURE);
1251                     wallcycle_start(wcycle, ewcUPDATE);
1252                 }
1253             }
1254             /* if it's the initial step, we performed this first step just to get the constraint virial */
1255             if (ir->eI == eiVV && bInitStep)
1256             {
1257                 copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
1258                 sfree(vbuf);
1259             }
1260             wallcycle_stop(wcycle, ewcUPDATE);
1261         }
1262
1263         /* compute the conserved quantity */
1264         if (EI_VV(ir->eI))
1265         {
1266             saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1267             if (ir->eI == eiVV)
1268             {
1269                 last_ekin = enerd->term[F_EKIN];
1270             }
1271             if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1272             {
1273                 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1274             }
1275             /* sum up the foreign energy and dhdl terms for vv.  currently done every step so that dhdl is correct in the .edr */
1276             if (ir->efep != efepNO && !bRerunMD)
1277             {
1278                 sum_dhdl(enerd, state->lambda, ir->fepvals);
1279             }
1280         }
1281
1282         /* ########  END FIRST UPDATE STEP  ############## */
1283         /* ########  If doing VV, we now have v(dt) ###### */
1284         if (bDoExpanded)
1285         {
1286             /* perform extended ensemble sampling in lambda - we don't
1287                actually move to the new state before outputting
1288                statistics, but if performing simulated tempering, we
1289                do update the velocities and the tau_t. */
1290
1291             lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
1292             /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1293             copy_df_history(state_global->dfhist, state->dfhist);
1294         }
1295
1296         /* Now we have the energies and forces corresponding to the
1297          * coordinates at time t. We must output all of this before
1298          * the update.
1299          */
1300         do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1301                                  ir, state, state_global, observablesHistory,
1302                                  top_global, fr,
1303                                  outf, mdebin, ekind, &f,
1304                                  &nchkpt,
1305                                  bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
1306                                  bSumEkinhOld);
1307         /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1308         bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
1309
1310         /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1311         if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1312         {
1313             copy_mat(state->svir_prev, shake_vir);
1314             copy_mat(state->fvir_prev, force_vir);
1315         }
1316
1317         elapsed_time = walltime_accounting_get_current_elapsed_time(walltime_accounting);
1318
1319         /* Check whether everything is still allright */
1320         if (((int)gmx_get_stop_condition() > handled_stop_condition)
1321 #if GMX_THREAD_MPI
1322             && MASTER(cr)
1323 #endif
1324             )
1325         {
1326             int nsteps_stop = -1;
1327
1328             /* this just makes signals[].sig compatible with the hack
1329                of sending signals around by MPI_Reduce together with
1330                other floats */
1331             if (gmx_get_stop_condition() == gmx_stop_cond_next_ns)
1332             {
1333                 signals[eglsSTOPCOND].sig = 1;
1334                 nsteps_stop               = std::max(ir->nstlist, 2*nstglobalcomm);
1335             }
1336             if (gmx_get_stop_condition() == gmx_stop_cond_next)
1337             {
1338                 signals[eglsSTOPCOND].sig = -1;
1339                 nsteps_stop               = nstglobalcomm + 1;
1340             }
1341             if (fplog)
1342             {
1343                 fprintf(fplog,
1344                         "\n\nReceived the %s signal, stopping within %d steps\n\n",
1345                         gmx_get_signal_name(), nsteps_stop);
1346                 fflush(fplog);
1347             }
1348             fprintf(stderr,
1349                     "\n\nReceived the %s signal, stopping within %d steps\n\n",
1350                     gmx_get_signal_name(), nsteps_stop);
1351             fflush(stderr);
1352             handled_stop_condition = (int)gmx_get_stop_condition();
1353         }
1354         else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1355                  (max_hours > 0 && elapsed_time > max_hours*60.0*60.0*0.99) &&
1356                  signals[eglsSTOPCOND].sig == 0 && signals[eglsSTOPCOND].set == 0)
1357         {
1358             /* Signal to terminate the run */
1359             signals[eglsSTOPCOND].sig = 1;
1360             if (fplog)
1361             {
1362                 fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
1363             }
1364             fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
1365         }
1366
1367         if (bResetCountersHalfMaxH && MASTER(cr) &&
1368             elapsed_time > max_hours*60.0*60.0*0.495)
1369         {
1370             /* Set flag that will communicate the signal to all ranks in the simulation */
1371             signals[eglsRESETCOUNTERS].sig = 1;
1372         }
1373
1374         /* In parallel we only have to check for checkpointing in steps
1375          * where we do global communication,
1376          *  otherwise the other nodes don't know.
1377          */
1378         if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1379                            cpt_period >= 0 &&
1380                            (cpt_period == 0 ||
1381                             elapsed_time >= nchkpt*cpt_period*60.0)) &&
1382             signals[eglsCHKPT].set == 0)
1383         {
1384             signals[eglsCHKPT].sig = 1;
1385         }
1386
1387         /* #########   START SECOND UPDATE STEP ################# */
1388
1389         /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1390            in preprocessing */
1391
1392         if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1393         {
1394             gmx_bool bIfRandomize;
1395             bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr);
1396             /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1397             if (constr && bIfRandomize)
1398             {
1399                 update_constraints(fplog, step, nullptr, ir, mdatoms,
1400                                    state, fr->bMolPBC, graph, &f,
1401                                    &top->idef, tmp_vir,
1402                                    cr, nrnb, wcycle, upd, constr,
1403                                    TRUE, bCalcVir);
1404             }
1405         }
1406         /* Box is changed in update() when we do pressure coupling,
1407          * but we should still use the old box for energy corrections and when
1408          * writing it to the energy file, so it matches the trajectory files for
1409          * the same timestep above. Make a copy in a separate array.
1410          */
1411         copy_mat(state->box, lastbox);
1412
1413         dvdl_constr = 0;
1414
1415         if (!bRerunMD || rerun_fr.bV || bForceUpdate)
1416         {
1417             wallcycle_start(wcycle, ewcUPDATE);
1418             /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1419             if (bTrotter)
1420             {
1421                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1422                 /* We can only do Berendsen coupling after we have summed
1423                  * the kinetic energy or virial. Since the happens
1424                  * in global_state after update, we should only do it at
1425                  * step % nstlist = 1 with bGStatEveryStep=FALSE.
1426                  */
1427             }
1428             else
1429             {
1430                 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1431                 update_pcouple_before_coordinates(fplog, step, ir, state,
1432                                                   parrinellorahmanMu, M,
1433                                                   bInitStep);
1434             }
1435
1436             if (EI_VV(ir->eI))
1437             {
1438                 /* velocity half-step update */
1439                 update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
1440                               ekind, M, upd, etrtVELOCITY2,
1441                               cr, constr);
1442             }
1443
1444             /* Above, initialize just copies ekinh into ekin,
1445              * it doesn't copy position (for VV),
1446              * and entire integrator for MD.
1447              */
1448
1449             if (ir->eI == eiVVAK)
1450             {
1451                 /* We probably only need md->homenr, not state->natoms */
1452                 if (state->natoms > cbuf_nalloc)
1453                 {
1454                     cbuf_nalloc = state->natoms;
1455                     srenew(cbuf, cbuf_nalloc);
1456                 }
1457                 copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1458             }
1459
1460             update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
1461                           ekind, M, upd, etrtPOSITION, cr, constr);
1462             wallcycle_stop(wcycle, ewcUPDATE);
1463
1464             update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
1465                                fr->bMolPBC, graph, &f,
1466                                &top->idef, shake_vir,
1467                                cr, nrnb, wcycle, upd, constr,
1468                                FALSE, bCalcVir);
1469
1470             if (ir->eI == eiVVAK)
1471             {
1472                 /* erase F_EKIN and F_TEMP here? */
1473                 /* just compute the kinetic energy at the half step to perform a trotter step */
1474                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1475                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1476                                 constr, &nullSignaller, lastbox,
1477                                 nullptr, &bSumEkinhOld,
1478                                 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1479                                 );
1480                 wallcycle_start(wcycle, ewcUPDATE);
1481                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1482                 /* now we know the scaling, we can compute the positions again again */
1483                 copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1484
1485                 update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
1486                               ekind, M, upd, etrtPOSITION, cr, constr);
1487                 wallcycle_stop(wcycle, ewcUPDATE);
1488
1489                 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1490                 /* are the small terms in the shake_vir here due
1491                  * to numerical errors, or are they important
1492                  * physically? I'm thinking they are just errors, but not completely sure.
1493                  * For now, will call without actually constraining, constr=NULL*/
1494                 update_constraints(fplog, step, nullptr, ir, mdatoms,
1495                                    state, fr->bMolPBC, graph, &f,
1496                                    &top->idef, tmp_vir,
1497                                    cr, nrnb, wcycle, upd, nullptr,
1498                                    FALSE, bCalcVir);
1499             }
1500             if (EI_VV(ir->eI))
1501             {
1502                 /* this factor or 2 correction is necessary
1503                    because half of the constraint force is removed
1504                    in the vv step, so we have to double it.  See
1505                    the Redmine issue #1255.  It is not yet clear
1506                    if the factor of 2 is exact, or just a very
1507                    good approximation, and this will be
1508                    investigated.  The next step is to see if this
1509                    can be done adding a dhdl contribution from the
1510                    rattle step, but this is somewhat more
1511                    complicated with the current code. Will be
1512                    investigated, hopefully for 4.6.3. However,
1513                    this current solution is much better than
1514                    having it completely wrong.
1515                  */
1516                 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1517             }
1518             else
1519             {
1520                 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1521             }
1522         }
1523         else if (graph)
1524         {
1525             /* Need to unshift here */
1526             unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1527         }
1528
1529         if (vsite != nullptr)
1530         {
1531             wallcycle_start(wcycle, ewcVSITECONSTR);
1532             if (graph != nullptr)
1533             {
1534                 shift_self(graph, state->box, as_rvec_array(state->x.data()));
1535             }
1536             construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
1537                              top->idef.iparams, top->idef.il,
1538                              fr->ePBC, fr->bMolPBC, cr, state->box);
1539
1540             if (graph != nullptr)
1541             {
1542                 unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1543             }
1544             wallcycle_stop(wcycle, ewcVSITECONSTR);
1545         }
1546
1547         /* ############## IF NOT VV, Calculate globals HERE  ############ */
1548         /* With Leap-Frog we can skip compute_globals at
1549          * non-communication steps, but we need to calculate
1550          * the kinetic energy one step before communication.
1551          */
1552         {
1553             // Organize to do inter-simulation signalling on steps if
1554             // and when algorithms require it.
1555             bool doInterSimSignal = (!bFirstStep && bDoReplEx) || bUsingEnsembleRestraints;
1556
1557             if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1558             {
1559                 // Since we're already communicating at this step, we
1560                 // can propagate intra-simulation signals. Note that
1561                 // check_nstglobalcomm has the responsibility for
1562                 // choosing the value of nstglobalcomm that is one way
1563                 // bGStat becomes true, so we can't get into a
1564                 // situation where e.g. checkpointing can't be
1565                 // signalled.
1566                 bool                doIntraSimSignal = true;
1567                 SimulationSignaller signaller(&signals, cr, doInterSimSignal, doIntraSimSignal);
1568
1569                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1570                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1571                                 constr, &signaller,
1572                                 lastbox,
1573                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1574                                 (bGStat ? CGLO_GSTAT : 0)
1575                                 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
1576                                 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1577                                 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1578                                 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1579                                 | CGLO_CONSTRAINT
1580                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1581                                 );
1582                 checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
1583                                                 top_global, top, state,
1584                                                 &shouldCheckNumberOfBondedInteractions);
1585             }
1586         }
1587
1588         /* #############  END CALC EKIN AND PRESSURE ################# */
1589
1590         /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1591            the virial that should probably be addressed eventually. state->veta has better properies,
1592            but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1593            generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
1594
1595         if (ir->efep != efepNO && (!EI_VV(ir->eI) || bRerunMD))
1596         {
1597             /* Sum up the foreign energy and dhdl terms for md and sd.
1598                Currently done every step so that dhdl is correct in the .edr */
1599             sum_dhdl(enerd, state->lambda, ir->fepvals);
1600         }
1601
1602         update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1603                                          pres, force_vir, shake_vir,
1604                                          parrinellorahmanMu,
1605                                          state, nrnb, upd);
1606
1607         /* ################# END UPDATE STEP 2 ################# */
1608         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
1609
1610         /* The coordinates (x) were unshifted in update */
1611         if (!bGStat)
1612         {
1613             /* We will not sum ekinh_old,
1614              * so signal that we still have to do it.
1615              */
1616             bSumEkinhOld = TRUE;
1617         }
1618
1619         if (bCalcEner)
1620         {
1621             /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
1622
1623             /* use the directly determined last velocity, not actually the averaged half steps */
1624             if (bTrotter && ir->eI == eiVV)
1625             {
1626                 enerd->term[F_EKIN] = last_ekin;
1627             }
1628             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1629
1630             if (integratorHasConservedEnergyQuantity(ir))
1631             {
1632                 if (EI_VV(ir->eI))
1633                 {
1634                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1635                 }
1636                 else
1637                 {
1638                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1639                 }
1640             }
1641             /* #########  END PREPARING EDR OUTPUT  ###########  */
1642         }
1643
1644         /* Output stuff */
1645         if (MASTER(cr))
1646         {
1647             if (fplog && do_log && bDoExpanded)
1648             {
1649                 /* only needed if doing expanded ensemble */
1650                 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1651                                           state_global->dfhist, state->fep_state, ir->nstlog, step);
1652             }
1653             if (bCalcEner)
1654             {
1655                 upd_mdebin(mdebin, bDoDHDL, bCalcEnerStep,
1656                            t, mdatoms->tmass, enerd, state,
1657                            ir->fepvals, ir->expandedvals, lastbox,
1658                            shake_vir, force_vir, total_vir, pres,
1659                            ekind, mu_tot, constr);
1660             }
1661             else
1662             {
1663                 upd_mdebin_step(mdebin);
1664             }
1665
1666             gmx_bool do_dr  = do_per_step(step, ir->nstdisreout);
1667             gmx_bool do_or  = do_per_step(step, ir->nstorireout);
1668
1669             print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
1670                        step, t,
1671                        eprNORMAL, mdebin, fcd, groups, &(ir->opts));
1672
1673             if (ir->bPull)
1674             {
1675                 pull_print_output(ir->pull_work, step, t);
1676             }
1677
1678             if (do_per_step(step, ir->nstlog))
1679             {
1680                 if (fflush(fplog) != 0)
1681                 {
1682                     gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1683                 }
1684             }
1685         }
1686         if (bDoExpanded)
1687         {
1688             /* Have to do this part _after_ outputting the logfile and the edr file */
1689             /* Gets written into the state at the beginning of next loop*/
1690             state->fep_state = lamnew;
1691         }
1692         /* Print the remaining wall clock time for the run */
1693         if (MULTIMASTER(cr) &&
1694             (do_verbose || gmx_got_usr_signal()) &&
1695             !bPMETunePrinting)
1696         {
1697             if (shellfc)
1698             {
1699                 fprintf(stderr, "\n");
1700             }
1701             print_time(stderr, walltime_accounting, step, ir, cr);
1702         }
1703
1704         /* Ion/water position swapping.
1705          * Not done in last step since trajectory writing happens before this call
1706          * in the MD loop and exchanges would be lost anyway. */
1707         bNeedRepartition = FALSE;
1708         if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1709             do_per_step(step, ir->swap->nstswap))
1710         {
1711             bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
1712                                              bRerunMD ? rerun_fr.x   : as_rvec_array(state->x.data()),
1713                                              bRerunMD ? rerun_fr.box : state->box,
1714                                              MASTER(cr) && bVerbose, bRerunMD);
1715
1716             if (bNeedRepartition && DOMAINDECOMP(cr))
1717             {
1718                 dd_collect_state(cr->dd, state, state_global);
1719             }
1720         }
1721
1722         /* Replica exchange */
1723         bExchanged = FALSE;
1724         if (bDoReplEx)
1725         {
1726             bExchanged = replica_exchange(fplog, cr, repl_ex,
1727                                           state_global, enerd,
1728                                           state, step, t);
1729         }
1730
1731         if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1732         {
1733             dd_partition_system(fplog, step, cr, TRUE, 1,
1734                                 state_global, top_global, ir,
1735                                 state, &f, mdatoms, top, fr,
1736                                 vsite, constr,
1737                                 nrnb, wcycle, FALSE);
1738             shouldCheckNumberOfBondedInteractions = true;
1739             update_realloc(upd, state->natoms);
1740         }
1741
1742         bFirstStep             = FALSE;
1743         bInitStep              = FALSE;
1744         startingFromCheckpoint = FALSE;
1745
1746         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1747         /* With all integrators, except VV, we need to retain the pressure
1748          * at the current step for coupling at the next step.
1749          */
1750         if ((state->flags & (1<<estPRES_PREV)) &&
1751             (bGStatEveryStep ||
1752              (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1753         {
1754             /* Store the pressure in t_state for pressure coupling
1755              * at the next MD step.
1756              */
1757             copy_mat(pres, state->pres_prev);
1758         }
1759
1760         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
1761
1762         if ( (membed != nullptr) && (!bLastStep) )
1763         {
1764             rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1765         }
1766
1767         if (bRerunMD)
1768         {
1769             if (MASTER(cr))
1770             {
1771                 /* read next frame from input trajectory */
1772                 bLastStep = !read_next_frame(oenv, status, &rerun_fr);
1773             }
1774
1775             if (PAR(cr))
1776             {
1777                 rerun_parallel_comm(cr, &rerun_fr, &bLastStep);
1778             }
1779         }
1780
1781         cycles = wallcycle_stop(wcycle, ewcSTEP);
1782         if (DOMAINDECOMP(cr) && wcycle)
1783         {
1784             dd_cycles_add(cr->dd, cycles, ddCyclStep);
1785         }
1786
1787         if (!bRerunMD || !rerun_fr.bStep)
1788         {
1789             /* increase the MD step number */
1790             step++;
1791             step_rel++;
1792         }
1793
1794         /* TODO make a counter-reset module */
1795         /* If it is time to reset counters, set a flag that remains
1796            true until counters actually get reset */
1797         if (step_rel == wcycle_get_reset_counters(wcycle) ||
1798             signals[eglsRESETCOUNTERS].set != 0)
1799         {
1800             if (pme_loadbal_is_active(pme_loadbal))
1801             {
1802                 /* Do not permit counter reset while PME load
1803                  * balancing is active. The only purpose for resetting
1804                  * counters is to measure reliable performance data,
1805                  * and that can't be done before balancing
1806                  * completes.
1807                  *
1808                  * TODO consider fixing this by delaying the reset
1809                  * until after load balancing completes,
1810                  * e.g. https://gerrit.gromacs.org/#/c/4964/2 */
1811                 gmx_fatal(FARGS, "PME tuning was still active when attempting to "
1812                           "reset mdrun counters at step %" GMX_PRId64 ". Try "
1813                           "resetting counters later in the run, e.g. with gmx "
1814                           "mdrun -resetstep.", step);
1815             }
1816             reset_all_counters(fplog, mdlog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
1817                                use_GPU(fr->nbv) ? fr->nbv : nullptr);
1818             wcycle_set_reset_counters(wcycle, -1);
1819             if (!(cr->duty & DUTY_PME))
1820             {
1821                 /* Tell our PME node to reset its counters */
1822                 gmx_pme_send_resetcounters(cr, step);
1823             }
1824             /* Correct max_hours for the elapsed time */
1825             max_hours                -= elapsed_time/(60.0*60.0);
1826             /* If mdrun -maxh -resethway was active, it can only trigger once */
1827             bResetCountersHalfMaxH    = FALSE; /* TODO move this to where signals[eglsRESETCOUNTERS].sig is set */
1828             /* Reset can only happen once, so clear the triggering flag. */
1829             signals[eglsRESETCOUNTERS].set = 0;
1830         }
1831
1832         /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1833         IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
1834
1835     }
1836     /* End of main MD loop */
1837
1838     /* Closing TNG files can include compressing data. Therefore it is good to do that
1839      * before stopping the time measurements. */
1840     mdoutf_tng_close(outf);
1841
1842     /* Stop measuring walltime */
1843     walltime_accounting_end(walltime_accounting);
1844
1845     if (bRerunMD && MASTER(cr))
1846     {
1847         close_trx(status);
1848     }
1849
1850     if (!(cr->duty & DUTY_PME))
1851     {
1852         /* Tell the PME only node to finish */
1853         gmx_pme_send_finish(cr);
1854     }
1855
1856     if (MASTER(cr))
1857     {
1858         if (ir->nstcalcenergy > 0 && !bRerunMD)
1859         {
1860             print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
1861                        eprAVER, mdebin, fcd, groups, &(ir->opts));
1862         }
1863     }
1864
1865     done_mdoutf(outf);
1866
1867     if (bPMETune)
1868     {
1869         pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
1870     }
1871
1872     done_shellfc(fplog, shellfc, step_rel);
1873
1874     if (useReplicaExchange && MASTER(cr))
1875     {
1876         print_replica_exchange_statistics(fplog, repl_ex);
1877     }
1878
1879     // Clean up swapcoords
1880     if (ir->eSwapCoords != eswapNO)
1881     {
1882         finish_swapcoords(ir->swap);
1883     }
1884
1885     /* Do essential dynamics cleanup if needed. Close .edo file */
1886     done_ed(&ed);
1887
1888     /* IMD cleanup, if bIMD is TRUE. */
1889     IMD_finalize(ir->bIMD, ir->imd);
1890
1891     walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1892     if (step_rel >= wcycle_get_reset_counters(wcycle) &&
1893         signals[eglsRESETCOUNTERS].set == 0 &&
1894         !bResetCountersHalfMaxH)
1895     {
1896         walltime_accounting_set_valid_finish(walltime_accounting);
1897     }
1898
1899     return 0;
1900 }