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->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     step     = ir->init_step;
 795     step_rel = 0;
 796
 797     // TODO extract this to new multi-simulation module
 798     if (MASTER(cr) && MULTISIM(cr) && !useReplicaExchange)
 799     {
 800         if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
 801         {
 802             GMX_LOG(mdlog.warning).appendText(
 803                     "Note: The number of steps is not consistent across multi simulations,\n"
 804                     "but we are proceeding anyway!");
 805         }
 806         if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
 807         {
 808             GMX_LOG(mdlog.warning).appendText(
 809                     "Note: The initial step is not consistent across multi simulations,\n"
 810                     "but we are proceeding anyway!");
 811         }
 812     }
 813
 814     /* and stop now if we should */
 815     bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
 816     while (!bLastStep)
 817     {
 818
 819         /* Determine if this is a neighbor search step */
 820         bNStList = (ir->nstlist > 0  && step % ir->nstlist == 0);
 821
 822         if (bPMETune && bNStList)
 823         {
 824             /* PME grid + cut-off optimization with GPUs or PME nodes */
 825             pme_loadbal_do(pme_loadbal, cr,
 826                            (bVerbose && MASTER(cr)) ? stderr : nullptr,
 827                            fplog, mdlog,
 828                            ir, fr, state,
 829                            wcycle,
 830                            step, step_rel,
 831                            &bPMETunePrinting);
 832         }
 833
 834         wallcycle_start(wcycle, ewcSTEP);
 835
 836         if (bRerunMD)
 837         {
 838             if (rerun_fr.bStep)
 839             {
 840                 step     = rerun_fr.step;
 841                 step_rel = step - ir->init_step;
 842             }
 843             if (rerun_fr.bTime)
 844             {
 845                 t = rerun_fr.time;
 846             }
 847             else
 848             {
 849                 t = step;
 850             }
 851         }
 852         else
 853         {
 854             bLastStep = (step_rel == ir->nsteps);
 855             t         = t0 + step*ir->delta_t;
 856         }
 857
 858         // TODO Refactor this, so that nstfep does not need a default value of zero
 859         if (ir->efep != efepNO || ir->bSimTemp)
 860         {
 861             /* find and set the current lambdas.  If rerunning, we either read in a state, or a lambda value,
 862                requiring different logic. */
 863
 864             set_current_lambdas(step, ir->fepvals, bRerunMD, &rerun_fr, state_global, state, lam0);
 865             bDoDHDL      = do_per_step(step, ir->fepvals->nstdhdl);
 866             bDoFEP       = ((ir->efep != efepNO) && do_per_step(step, nstfep));
 867             bDoExpanded  = (do_per_step(step, ir->expandedvals->nstexpanded)
 868                             && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
 869         }
 870
 871         bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
 872                      do_per_step(step, replExParams.exchangeInterval));
 873
 874         if (bSimAnn)
 875         {
 876             update_annealing_target_temp(ir, t, upd);
 877         }
 878
 879         if (bRerunMD)
 880         {
 881             if (!DOMAINDECOMP(cr) || MASTER(cr))
 882             {
 883                 for (i = 0; i < state_global->natoms; i++)
 884                 {
 885                     copy_rvec(rerun_fr.x[i], state_global->x[i]);
 886                 }
 887                 if (rerun_fr.bV)
 888                 {
 889                     for (i = 0; i < state_global->natoms; i++)
 890                     {
 891                         copy_rvec(rerun_fr.v[i], state_global->v[i]);
 892                     }
 893                 }
 894                 else
 895                 {
 896                     for (i = 0; i < state_global->natoms; i++)
 897                     {
 898                         clear_rvec(state_global->v[i]);
 899                     }
 900                     if (bRerunWarnNoV)
 901                     {
 902                         fprintf(stderr, "\nWARNING: Some frames do not contain velocities.\n"
 903                                 "         Ekin, temperature and pressure are incorrect,\n"
 904                                 "         the virial will be incorrect when constraints are present.\n"
 905                                 "\n");
 906                         bRerunWarnNoV = FALSE;
 907                     }
 908                 }
 909             }
 910             copy_mat(rerun_fr.box, state_global->box);
 911             copy_mat(state_global->box, state->box);
 912
 913             if (vsite && (Flags & MD_RERUN_VSITE))
 914             {
 915                 if (DOMAINDECOMP(cr))
 916                 {
 917                     gmx_fatal(FARGS, "Vsite recalculation with -rerun is not implemented with domain decomposition, use a single rank");
 918                 }
 919                 if (graph)
 920                 {
 921                     /* Following is necessary because the graph may get out of sync
 922                      * with the coordinates if we only have every N'th coordinate set
 923                      */
 924                     mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
 925                     shift_self(graph, state->box, as_rvec_array(state->x.data()));
 926                 }
 927                 construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
 928                                  top->idef.iparams, top->idef.il,
 929                                  fr->ePBC, fr->bMolPBC, cr, state->box);
 930                 if (graph)
 931                 {
 932                     unshift_self(graph, state->box, as_rvec_array(state->x.data()));
 933                 }
 934             }
 935         }
 936
 937         /* Stop Center of Mass motion */
 938         bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
 939
 940         if (bRerunMD)
 941         {
 942             /* for rerun MD always do Neighbour Searching */
 943             bNS      = (bFirstStep || ir->nstlist != 0);
 944         }
 945         else
 946         {
 947             /* Determine whether or not to do Neighbour Searching */
 948             bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
 949         }
 950
 951         /* < 0 means stop at next step, > 0 means stop at next NS step */
 952         if ( (signals[eglsSTOPCOND].set < 0) ||
 953              ( (signals[eglsSTOPCOND].set > 0 ) && ( bNS || ir->nstlist == 0)))
 954         {
 955             bLastStep = TRUE;
 956         }
 957
 958         /* Determine whether or not to update the Born radii if doing GB */
 959         bBornRadii = bFirstStep;
 960         if (ir->implicit_solvent && (step % ir->nstgbradii == 0))
 961         {
 962             bBornRadii = TRUE;
 963         }
 964
 965         /* do_log triggers energy and virial calculation. Because this leads
 966          * to different code paths, forces can be different. Thus for exact
 967          * continuation we should avoid extra log output.
 968          * Note that the || bLastStep can result in non-exact continuation
 969          * beyond the last step. But we don't consider that to be an issue.
 970          */
 971         do_log     = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep || bRerunMD;
 972         do_verbose = bVerbose &&
 973             (step % stepout == 0 || bFirstStep || bLastStep || bRerunMD);
 974
 975         if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
 976         {
 977             if (bRerunMD)
 978             {
 979                 bMasterState = TRUE;
 980             }
 981             else
 982             {
 983                 bMasterState = FALSE;
 984                 /* Correct the new box if it is too skewed */
 985                 if (inputrecDynamicBox(ir))
 986                 {
 987                     if (correct_box(fplog, step, state->box, graph))
 988                     {
 989                         bMasterState = TRUE;
 990                     }
 991                 }
 992                 if (DOMAINDECOMP(cr) && bMasterState)
 993                 {
 994                     dd_collect_state(cr->dd, state, state_global);
 995                 }
 996             }
 997
 998             if (DOMAINDECOMP(cr))
 999             {
1000                 /* Repartition the domain decomposition */
1001                 dd_partition_system(fplog, step, cr,
1002                                     bMasterState, nstglobalcomm,
1003                                     state_global, top_global, ir,
1004                                     state, &f, mdatoms, top, fr,
1005                                     vsite, constr,
1006                                     nrnb, wcycle,
1007                                     do_verbose && !bPMETunePrinting);
1008                 shouldCheckNumberOfBondedInteractions = true;
1009                 update_realloc(upd, state->natoms);
1010             }
1011         }
1012
1013         if (MASTER(cr) && do_log)
1014         {
1015             print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
1016         }
1017
1018         if (ir->efep != efepNO)
1019         {
1020             update_mdatoms(mdatoms, state->lambda[efptMASS]);
1021         }
1022
1023         if ((bRerunMD && rerun_fr.bV) || bExchanged)
1024         {
1025
1026             /* We need the kinetic energy at minus the half step for determining
1027              * the full step kinetic energy and possibly for T-coupling.*/
1028             /* This may not be quite working correctly yet . . . . */
1029             compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1030                             wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1031                             constr, &nullSignaller, state->box,
1032                             &totalNumberOfBondedInteractions, &bSumEkinhOld,
1033                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
1034             checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
1035                                             top_global, top, state,
1036                                             &shouldCheckNumberOfBondedInteractions);
1037         }
1038         clear_mat(force_vir);
1039
1040         /* We write a checkpoint at this MD step when:
1041          * either at an NS step when we signalled through gs,
1042          * or at the last step (but not when we do not want confout),
1043          * but never at the first step or with rerun.
1044          */
1045         bCPT = (((signals[eglsCHKPT].set && (bNS || ir->nstlist == 0)) ||
1046                  (bLastStep && (Flags & MD_CONFOUT))) &&
1047                 step > ir->init_step && !bRerunMD);
1048         if (bCPT)
1049         {
1050             signals[eglsCHKPT].set = 0;
1051         }
1052
1053         /* Determine the energy and pressure:
1054          * at nstcalcenergy steps and at energy output steps (set below).
1055          */
1056         if (EI_VV(ir->eI) && (!bInitStep))
1057         {
1058             /* for vv, the first half of the integration actually corresponds
1059                to the previous step.  bCalcEner is only required to be evaluated on the 'next' step,
1060                but the virial needs to be calculated on both the current step and the 'next' step. Future
1061                reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
1062
1063             /* TODO: This is probably not what we want, we will write to energy file one step after nstcalcenergy steps. */
1064             bCalcEnerStep = do_per_step(step - 1, ir->nstcalcenergy);
1065             bCalcVir      = bCalcEnerStep ||
1066                 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
1067         }
1068         else
1069         {
1070             bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
1071             bCalcVir      = bCalcEnerStep ||
1072                 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
1073         }
1074         bCalcEner = bCalcEnerStep;
1075
1076         do_ene = (do_per_step(step, ir->nstenergy) || bLastStep || bRerunMD);
1077
1078         if (do_ene || do_log || bDoReplEx)
1079         {
1080             bCalcVir  = TRUE;
1081             bCalcEner = TRUE;
1082         }
1083
1084         /* Do we need global communication ? */
1085         bGStat = (bCalcVir || bCalcEner || bStopCM ||
1086                   do_per_step(step, nstglobalcomm) ||
1087                   (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
1088
1089         force_flags = (GMX_FORCE_STATECHANGED |
1090                        ((inputrecDynamicBox(ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1091                        GMX_FORCE_ALLFORCES |
1092                        (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
1093                        (bCalcEner ? GMX_FORCE_ENERGY : 0) |
1094                        (bDoFEP ? GMX_FORCE_DHDL : 0)
1095                        );
1096
1097         if (shellfc)
1098         {
1099             /* Now is the time to relax the shells */
1100             relax_shell_flexcon(fplog, cr, bVerbose, step,
1101                                 ir, bNS, force_flags, top,
1102                                 constr, enerd, fcd,
1103                                 state, &f, force_vir, mdatoms,
1104                                 nrnb, wcycle, graph, groups,
1105                                 shellfc, fr, bBornRadii, t, mu_tot,
1106                                 vsite);
1107         }
1108         else
1109         {
1110             /* The coordinates (x) are shifted (to get whole molecules)
1111              * in do_force.
1112              * This is parallellized as well, and does communication too.
1113              * Check comments in sim_util.c
1114              */
1115             do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
1116                      state->box, &state->x, &state->hist,
1117                      &f, force_vir, mdatoms, enerd, fcd,
1118                      state->lambda, graph,
1119                      fr, vsite, mu_tot, t, ed, bBornRadii,
1120                      (bNS ? GMX_FORCE_NS : 0) | force_flags);
1121         }
1122
1123         if (EI_VV(ir->eI) && !startingFromCheckpoint && !bRerunMD)
1124         /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1125         {
1126             rvec *vbuf = nullptr;
1127
1128             wallcycle_start(wcycle, ewcUPDATE);
1129             if (ir->eI == eiVV && bInitStep)
1130             {
1131                 /* if using velocity verlet with full time step Ekin,
1132                  * take the first half step only to compute the
1133                  * virial for the first step. From there,
1134                  * revert back to the initial coordinates
1135                  * so that the input is actually the initial step.
1136                  */
1137                 snew(vbuf, state->natoms);
1138                 copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
1139             }
1140             else
1141             {
1142                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1143                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
1144             }
1145
1146             update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
1147                           ekind, M, upd, etrtVELOCITY1,
1148                           cr, constr);
1149
1150             if (!bRerunMD || rerun_fr.bV || bForceUpdate)         /* Why is rerun_fr.bV here?  Unclear. */
1151             {
1152                 wallcycle_stop(wcycle, ewcUPDATE);
1153                 update_constraints(fplog, step, nullptr, ir, mdatoms,
1154                                    state, fr->bMolPBC, graph, &f,
1155                                    &top->idef, shake_vir,
1156                                    cr, nrnb, wcycle, upd, constr,
1157                                    TRUE, bCalcVir);
1158                 wallcycle_start(wcycle, ewcUPDATE);
1159             }
1160             else if (graph)
1161             {
1162                 /* Need to unshift here if a do_force has been
1163                    called in the previous step */
1164                 unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1165             }
1166             /* if VV, compute the pressure and constraints */
1167             /* For VV2, we strictly only need this if using pressure
1168              * control, but we really would like to have accurate pressures
1169              * printed out.
1170              * Think about ways around this in the future?
1171              * For now, keep this choice in comments.
1172              */
1173             /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
1174             /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
1175             bPres = TRUE;
1176             bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
1177             if (bCalcEner && ir->eI == eiVVAK)
1178             {
1179                 bSumEkinhOld = TRUE;
1180             }
1181             /* for vv, the first half of the integration actually corresponds to the previous step.
1182                So we need information from the last step in the first half of the integration */
1183             if (bGStat || do_per_step(step-1, nstglobalcomm))
1184             {
1185                 wallcycle_stop(wcycle, ewcUPDATE);
1186                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1187                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1188                                 constr, &nullSignaller, state->box,
1189                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1190                                 (bGStat ? CGLO_GSTAT : 0)
1191                                 | CGLO_ENERGY
1192                                 | (bTemp ? CGLO_TEMPERATURE : 0)
1193                                 | (bPres ? CGLO_PRESSURE : 0)
1194                                 | (bPres ? CGLO_CONSTRAINT : 0)
1195                                 | (bStopCM ? CGLO_STOPCM : 0)
1196                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1197                                 | CGLO_SCALEEKIN
1198                                 );
1199                 /* explanation of above:
1200                    a) We compute Ekin at the full time step
1201                    if 1) we are using the AveVel Ekin, and it's not the
1202                    initial step, or 2) if we are using AveEkin, but need the full
1203                    time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1204                    b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1205                    EkinAveVel because it's needed for the pressure */
1206                 checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
1207                                                 top_global, top, state,
1208                                                 &shouldCheckNumberOfBondedInteractions);
1209                 wallcycle_start(wcycle, ewcUPDATE);
1210             }
1211             /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1212             if (!bInitStep)
1213             {
1214                 if (bTrotter)
1215                 {
1216                     m_add(force_vir, shake_vir, total_vir);     /* we need the un-dispersion corrected total vir here */
1217                     trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
1218
1219                     /* TODO This is only needed when we're about to write
1220                      * a checkpoint, because we use it after the restart
1221                      * (in a kludge?). But what should we be doing if
1222                      * startingFromCheckpoint or bInitStep are true? */
1223                     if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1224                     {
1225                         copy_mat(shake_vir, state->svir_prev);
1226                         copy_mat(force_vir, state->fvir_prev);
1227                     }
1228                     if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
1229                     {
1230                         /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1231                         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
1232                         enerd->term[F_EKIN] = trace(ekind->ekin);
1233                     }
1234                 }
1235                 else if (bExchanged)
1236                 {
1237                     wallcycle_stop(wcycle, ewcUPDATE);
1238                     /* We need the kinetic energy at minus the half step for determining
1239                      * the full step kinetic energy and possibly for T-coupling.*/
1240                     /* This may not be quite working correctly yet . . . . */
1241                     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1242                                     wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1243                                     constr, &nullSignaller, state->box,
1244                                     nullptr, &bSumEkinhOld,
1245                                     CGLO_GSTAT | CGLO_TEMPERATURE);
1246                     wallcycle_start(wcycle, ewcUPDATE);
1247                 }
1248             }
1249             /* if it's the initial step, we performed this first step just to get the constraint virial */
1250             if (ir->eI == eiVV && bInitStep)
1251             {
1252                 copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
1253                 sfree(vbuf);
1254             }
1255             wallcycle_stop(wcycle, ewcUPDATE);
1256         }
1257
1258         /* compute the conserved quantity */
1259         if (EI_VV(ir->eI))
1260         {
1261             saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1262             if (ir->eI == eiVV)
1263             {
1264                 last_ekin = enerd->term[F_EKIN];
1265             }
1266             if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1267             {
1268                 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1269             }
1270             /* sum up the foreign energy and dhdl terms for vv.  currently done every step so that dhdl is correct in the .edr */
1271             if (ir->efep != efepNO && !bRerunMD)
1272             {
1273                 sum_dhdl(enerd, state->lambda, ir->fepvals);
1274             }
1275         }
1276
1277         /* ########  END FIRST UPDATE STEP  ############## */
1278         /* ########  If doing VV, we now have v(dt) ###### */
1279         if (bDoExpanded)
1280         {
1281             /* perform extended ensemble sampling in lambda - we don't
1282                actually move to the new state before outputting
1283                statistics, but if performing simulated tempering, we
1284                do update the velocities and the tau_t. */
1285
1286             lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
1287             /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1288             copy_df_history(state_global->dfhist, state->dfhist);
1289         }
1290
1291         /* Now we have the energies and forces corresponding to the
1292          * coordinates at time t. We must output all of this before
1293          * the update.
1294          */
1295         do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1296                                  ir, state, state_global, observablesHistory,
1297                                  top_global, fr,
1298                                  outf, mdebin, ekind, &f,
1299                                  &nchkpt,
1300                                  bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
1301                                  bSumEkinhOld);
1302         /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1303         bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
1304
1305         /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1306         if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1307         {
1308             copy_mat(state->svir_prev, shake_vir);
1309             copy_mat(state->fvir_prev, force_vir);
1310         }
1311
1312         elapsed_time = walltime_accounting_get_current_elapsed_time(walltime_accounting);
1313
1314         /* Check whether everything is still allright */
1315         if (((int)gmx_get_stop_condition() > handled_stop_condition)
1316 #if GMX_THREAD_MPI
1317             && MASTER(cr)
1318 #endif
1319             )
1320         {
1321             int nsteps_stop = -1;
1322
1323             /* this just makes signals[].sig compatible with the hack
1324                of sending signals around by MPI_Reduce together with
1325                other floats */
1326             if (gmx_get_stop_condition() == gmx_stop_cond_next_ns)
1327             {
1328                 signals[eglsSTOPCOND].sig = 1;
1329                 nsteps_stop               = std::max(ir->nstlist, 2*nstglobalcomm);
1330             }
1331             if (gmx_get_stop_condition() == gmx_stop_cond_next)
1332             {
1333                 signals[eglsSTOPCOND].sig = -1;
1334                 nsteps_stop               = nstglobalcomm + 1;
1335             }
1336             if (fplog)
1337             {
1338                 fprintf(fplog,
1339                         "\n\nReceived the %s signal, stopping within %d steps\n\n",
1340                         gmx_get_signal_name(), nsteps_stop);
1341                 fflush(fplog);
1342             }
1343             fprintf(stderr,
1344                     "\n\nReceived the %s signal, stopping within %d steps\n\n",
1345                     gmx_get_signal_name(), nsteps_stop);
1346             fflush(stderr);
1347             handled_stop_condition = (int)gmx_get_stop_condition();
1348         }
1349         else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1350                  (max_hours > 0 && elapsed_time > max_hours*60.0*60.0*0.99) &&
1351                  signals[eglsSTOPCOND].sig == 0 && signals[eglsSTOPCOND].set == 0)
1352         {
1353             /* Signal to terminate the run */
1354             signals[eglsSTOPCOND].sig = 1;
1355             if (fplog)
1356             {
1357                 fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
1358             }
1359             fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
1360         }
1361
1362         if (bResetCountersHalfMaxH && MASTER(cr) &&
1363             elapsed_time > max_hours*60.0*60.0*0.495)
1364         {
1365             /* Set flag that will communicate the signal to all ranks in the simulation */
1366             signals[eglsRESETCOUNTERS].sig = 1;
1367         }
1368
1369         /* In parallel we only have to check for checkpointing in steps
1370          * where we do global communication,
1371          *  otherwise the other nodes don't know.
1372          */
1373         if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1374                            cpt_period >= 0 &&
1375                            (cpt_period == 0 ||
1376                             elapsed_time >= nchkpt*cpt_period*60.0)) &&
1377             signals[eglsCHKPT].set == 0)
1378         {
1379             signals[eglsCHKPT].sig = 1;
1380         }
1381
1382         /* #########   START SECOND UPDATE STEP ################# */
1383
1384         /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1385            in preprocessing */
1386
1387         if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1388         {
1389             gmx_bool bIfRandomize;
1390             bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr);
1391             /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1392             if (constr && bIfRandomize)
1393             {
1394                 update_constraints(fplog, step, nullptr, ir, mdatoms,
1395                                    state, fr->bMolPBC, graph, &f,
1396                                    &top->idef, tmp_vir,
1397                                    cr, nrnb, wcycle, upd, constr,
1398                                    TRUE, bCalcVir);
1399             }
1400         }
1401         /* Box is changed in update() when we do pressure coupling,
1402          * but we should still use the old box for energy corrections and when
1403          * writing it to the energy file, so it matches the trajectory files for
1404          * the same timestep above. Make a copy in a separate array.
1405          */
1406         copy_mat(state->box, lastbox);
1407
1408         dvdl_constr = 0;
1409
1410         if (!bRerunMD || rerun_fr.bV || bForceUpdate)
1411         {
1412             wallcycle_start(wcycle, ewcUPDATE);
1413             /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1414             if (bTrotter)
1415             {
1416                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1417                 /* We can only do Berendsen coupling after we have summed
1418                  * the kinetic energy or virial. Since the happens
1419                  * in global_state after update, we should only do it at
1420                  * step % nstlist = 1 with bGStatEveryStep=FALSE.
1421                  */
1422             }
1423             else
1424             {
1425                 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1426                 update_pcouple_before_coordinates(fplog, step, ir, state,
1427                                                   parrinellorahmanMu, M,
1428                                                   bInitStep);
1429             }
1430
1431             if (EI_VV(ir->eI))
1432             {
1433                 /* velocity half-step update */
1434                 update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
1435                               ekind, M, upd, etrtVELOCITY2,
1436                               cr, constr);
1437             }
1438
1439             /* Above, initialize just copies ekinh into ekin,
1440              * it doesn't copy position (for VV),
1441              * and entire integrator for MD.
1442              */
1443
1444             if (ir->eI == eiVVAK)
1445             {
1446                 /* We probably only need md->homenr, not state->natoms */
1447                 if (state->natoms > cbuf_nalloc)
1448                 {
1449                     cbuf_nalloc = state->natoms;
1450                     srenew(cbuf, cbuf_nalloc);
1451                 }
1452                 copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1453             }
1454
1455             update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
1456                           ekind, M, upd, etrtPOSITION, cr, constr);
1457             wallcycle_stop(wcycle, ewcUPDATE);
1458
1459             update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
1460                                fr->bMolPBC, graph, &f,
1461                                &top->idef, shake_vir,
1462                                cr, nrnb, wcycle, upd, constr,
1463                                FALSE, bCalcVir);
1464
1465             if (ir->eI == eiVVAK)
1466             {
1467                 /* erase F_EKIN and F_TEMP here? */
1468                 /* just compute the kinetic energy at the half step to perform a trotter step */
1469                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1470                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1471                                 constr, &nullSignaller, lastbox,
1472                                 nullptr, &bSumEkinhOld,
1473                                 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1474                                 );
1475                 wallcycle_start(wcycle, ewcUPDATE);
1476                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1477                 /* now we know the scaling, we can compute the positions again again */
1478                 copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1479
1480                 update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
1481                               ekind, M, upd, etrtPOSITION, cr, constr);
1482                 wallcycle_stop(wcycle, ewcUPDATE);
1483
1484                 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1485                 /* are the small terms in the shake_vir here due
1486                  * to numerical errors, or are they important
1487                  * physically? I'm thinking they are just errors, but not completely sure.
1488                  * For now, will call without actually constraining, constr=NULL*/
1489                 update_constraints(fplog, step, nullptr, ir, mdatoms,
1490                                    state, fr->bMolPBC, graph, &f,
1491                                    &top->idef, tmp_vir,
1492                                    cr, nrnb, wcycle, upd, nullptr,
1493                                    FALSE, bCalcVir);
1494             }
1495             if (EI_VV(ir->eI))
1496             {
1497                 /* this factor or 2 correction is necessary
1498                    because half of the constraint force is removed
1499                    in the vv step, so we have to double it.  See
1500                    the Redmine issue #1255.  It is not yet clear
1501                    if the factor of 2 is exact, or just a very
1502                    good approximation, and this will be
1503                    investigated.  The next step is to see if this
1504                    can be done adding a dhdl contribution from the
1505                    rattle step, but this is somewhat more
1506                    complicated with the current code. Will be
1507                    investigated, hopefully for 4.6.3. However,
1508                    this current solution is much better than
1509                    having it completely wrong.
1510                  */
1511                 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1512             }
1513             else
1514             {
1515                 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1516             }
1517         }
1518         else if (graph)
1519         {
1520             /* Need to unshift here */
1521             unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1522         }
1523
1524         if (vsite != nullptr)
1525         {
1526             wallcycle_start(wcycle, ewcVSITECONSTR);
1527             if (graph != nullptr)
1528             {
1529                 shift_self(graph, state->box, as_rvec_array(state->x.data()));
1530             }
1531             construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
1532                              top->idef.iparams, top->idef.il,
1533                              fr->ePBC, fr->bMolPBC, cr, state->box);
1534
1535             if (graph != nullptr)
1536             {
1537                 unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1538             }
1539             wallcycle_stop(wcycle, ewcVSITECONSTR);
1540         }
1541
1542         /* ############## IF NOT VV, Calculate globals HERE  ############ */
1543         /* With Leap-Frog we can skip compute_globals at
1544          * non-communication steps, but we need to calculate
1545          * the kinetic energy one step before communication.
1546          */
1547         {
1548             // Organize to do inter-simulation signalling on steps if
1549             // and when algorithms require it.
1550             bool doInterSimSignal = (!bFirstStep && bDoReplEx) || bUsingEnsembleRestraints;
1551
1552             if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1553             {
1554                 // Since we're already communicating at this step, we
1555                 // can propagate intra-simulation signals. Note that
1556                 // check_nstglobalcomm has the responsibility for
1557                 // choosing the value of nstglobalcomm that is one way
1558                 // bGStat becomes true, so we can't get into a
1559                 // situation where e.g. checkpointing can't be
1560                 // signalled.
1561                 bool                doIntraSimSignal = true;
1562                 SimulationSignaller signaller(&signals, cr, doInterSimSignal, doIntraSimSignal);
1563
1564                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1565                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1566                                 constr, &signaller,
1567                                 lastbox,
1568                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1569                                 (bGStat ? CGLO_GSTAT : 0)
1570                                 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
1571                                 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1572                                 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1573                                 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1574                                 | CGLO_CONSTRAINT
1575                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1576                                 );
1577                 checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
1578                                                 top_global, top, state,
1579                                                 &shouldCheckNumberOfBondedInteractions);
1580             }
1581         }
1582
1583         /* #############  END CALC EKIN AND PRESSURE ################# */
1584
1585         /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1586            the virial that should probably be addressed eventually. state->veta has better properies,
1587            but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1588            generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
1589
1590         if (ir->efep != efepNO && (!EI_VV(ir->eI) || bRerunMD))
1591         {
1592             /* Sum up the foreign energy and dhdl terms for md and sd.
1593                Currently done every step so that dhdl is correct in the .edr */
1594             sum_dhdl(enerd, state->lambda, ir->fepvals);
1595         }
1596
1597         update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1598                                          pres, force_vir, shake_vir,
1599                                          parrinellorahmanMu,
1600                                          state, nrnb, upd);
1601
1602         /* ################# END UPDATE STEP 2 ################# */
1603         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
1604
1605         /* The coordinates (x) were unshifted in update */
1606         if (!bGStat)
1607         {
1608             /* We will not sum ekinh_old,
1609              * so signal that we still have to do it.
1610              */
1611             bSumEkinhOld = TRUE;
1612         }
1613
1614         if (bCalcEner)
1615         {
1616             /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
1617
1618             /* use the directly determined last velocity, not actually the averaged half steps */
1619             if (bTrotter && ir->eI == eiVV)
1620             {
1621                 enerd->term[F_EKIN] = last_ekin;
1622             }
1623             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1624
1625             if (integratorHasConservedEnergyQuantity(ir))
1626             {
1627                 if (EI_VV(ir->eI))
1628                 {
1629                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1630                 }
1631                 else
1632                 {
1633                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1634                 }
1635             }
1636             /* #########  END PREPARING EDR OUTPUT  ###########  */
1637         }
1638
1639         /* Output stuff */
1640         if (MASTER(cr))
1641         {
1642             if (fplog && do_log && bDoExpanded)
1643             {
1644                 /* only needed if doing expanded ensemble */
1645                 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1646                                           state_global->dfhist, state->fep_state, ir->nstlog, step);
1647             }
1648             if (bCalcEner)
1649             {
1650                 upd_mdebin(mdebin, bDoDHDL, bCalcEnerStep,
1651                            t, mdatoms->tmass, enerd, state,
1652                            ir->fepvals, ir->expandedvals, lastbox,
1653                            shake_vir, force_vir, total_vir, pres,
1654                            ekind, mu_tot, constr);
1655             }
1656             else
1657             {
1658                 upd_mdebin_step(mdebin);
1659             }
1660
1661             gmx_bool do_dr  = do_per_step(step, ir->nstdisreout);
1662             gmx_bool do_or  = do_per_step(step, ir->nstorireout);
1663
1664             print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
1665                        step, t,
1666                        eprNORMAL, mdebin, fcd, groups, &(ir->opts));
1667
1668             if (ir->bPull)
1669             {
1670                 pull_print_output(ir->pull_work, step, t);
1671             }
1672
1673             if (do_per_step(step, ir->nstlog))
1674             {
1675                 if (fflush(fplog) != 0)
1676                 {
1677                     gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1678                 }
1679             }
1680         }
1681         if (bDoExpanded)
1682         {
1683             /* Have to do this part _after_ outputting the logfile and the edr file */
1684             /* Gets written into the state at the beginning of next loop*/
1685             state->fep_state = lamnew;
1686         }
1687         /* Print the remaining wall clock time for the run */
1688         if (MULTIMASTER(cr) &&
1689             (do_verbose || gmx_got_usr_signal()) &&
1690             !bPMETunePrinting)
1691         {
1692             if (shellfc)
1693             {
1694                 fprintf(stderr, "\n");
1695             }
1696             print_time(stderr, walltime_accounting, step, ir, cr);
1697         }
1698
1699         /* Ion/water position swapping.
1700          * Not done in last step since trajectory writing happens before this call
1701          * in the MD loop and exchanges would be lost anyway. */
1702         bNeedRepartition = FALSE;
1703         if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1704             do_per_step(step, ir->swap->nstswap))
1705         {
1706             bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
1707                                              bRerunMD ? rerun_fr.x   : as_rvec_array(state->x.data()),
1708                                              bRerunMD ? rerun_fr.box : state->box,
1709                                              MASTER(cr) && bVerbose, bRerunMD);
1710
1711             if (bNeedRepartition && DOMAINDECOMP(cr))
1712             {
1713                 dd_collect_state(cr->dd, state, state_global);
1714             }
1715         }
1716
1717         /* Replica exchange */
1718         bExchanged = FALSE;
1719         if (bDoReplEx)
1720         {
1721             bExchanged = replica_exchange(fplog, cr, repl_ex,
1722                                           state_global, enerd,
1723                                           state, step, t);
1724         }
1725
1726         if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1727         {
1728             dd_partition_system(fplog, step, cr, TRUE, 1,
1729                                 state_global, top_global, ir,
1730                                 state, &f, mdatoms, top, fr,
1731                                 vsite, constr,
1732                                 nrnb, wcycle, FALSE);
1733             shouldCheckNumberOfBondedInteractions = true;
1734             update_realloc(upd, state->natoms);
1735         }
1736
1737         bFirstStep             = FALSE;
1738         bInitStep              = FALSE;
1739         startingFromCheckpoint = FALSE;
1740
1741         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1742         /* With all integrators, except VV, we need to retain the pressure
1743          * at the current step for coupling at the next step.
1744          */
1745         if ((state->flags & (1<<estPRES_PREV)) &&
1746             (bGStatEveryStep ||
1747              (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1748         {
1749             /* Store the pressure in t_state for pressure coupling
1750              * at the next MD step.
1751              */
1752             copy_mat(pres, state->pres_prev);
1753         }
1754
1755         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
1756
1757         if ( (membed != nullptr) && (!bLastStep) )
1758         {
1759             rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1760         }
1761
1762         if (bRerunMD)
1763         {
1764             if (MASTER(cr))
1765             {
1766                 /* read next frame from input trajectory */
1767                 bLastStep = !read_next_frame(oenv, status, &rerun_fr);
1768             }
1769
1770             if (PAR(cr))
1771             {
1772                 rerun_parallel_comm(cr, &rerun_fr, &bLastStep);
1773             }
1774         }
1775
1776         cycles = wallcycle_stop(wcycle, ewcSTEP);
1777         if (DOMAINDECOMP(cr) && wcycle)
1778         {
1779             dd_cycles_add(cr->dd, cycles, ddCyclStep);
1780         }
1781
1782         if (!bRerunMD || !rerun_fr.bStep)
1783         {
1784             /* increase the MD step number */
1785             step++;
1786             step_rel++;
1787         }
1788
1789         /* TODO make a counter-reset module */
1790         /* If it is time to reset counters, set a flag that remains
1791            true until counters actually get reset */
1792         if (step_rel == wcycle_get_reset_counters(wcycle) ||
1793             signals[eglsRESETCOUNTERS].set != 0)
1794         {
1795             if (pme_loadbal_is_active(pme_loadbal))
1796             {
1797                 /* Do not permit counter reset while PME load
1798                  * balancing is active. The only purpose for resetting
1799                  * counters is to measure reliable performance data,
1800                  * and that can't be done before balancing
1801                  * completes.
1802                  *
1803                  * TODO consider fixing this by delaying the reset
1804                  * until after load balancing completes,
1805                  * e.g. https://gerrit.gromacs.org/#/c/4964/2 */
1806                 gmx_fatal(FARGS, "PME tuning was still active when attempting to "
1807                           "reset mdrun counters at step %" GMX_PRId64 ". Try "
1808                           "resetting counters later in the run, e.g. with gmx "
1809                           "mdrun -resetstep.", step);
1810             }
1811             reset_all_counters(fplog, mdlog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
1812                                use_GPU(fr->nbv) ? fr->nbv : nullptr);
1813             wcycle_set_reset_counters(wcycle, -1);
1814             if (!(cr->duty & DUTY_PME))
1815             {
1816                 /* Tell our PME node to reset its counters */
1817                 gmx_pme_send_resetcounters(cr, step);
1818             }
1819             /* Correct max_hours for the elapsed time */
1820             max_hours                -= elapsed_time/(60.0*60.0);
1821             /* If mdrun -maxh -resethway was active, it can only trigger once */
1822             bResetCountersHalfMaxH    = FALSE; /* TODO move this to where signals[eglsRESETCOUNTERS].sig is set */
1823             /* Reset can only happen once, so clear the triggering flag. */
1824             signals[eglsRESETCOUNTERS].set = 0;
1825         }
1826
1827         /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1828         IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
1829
1830     }
1831     /* End of main MD loop */
1832
1833     /* Closing TNG files can include compressing data. Therefore it is good to do that
1834      * before stopping the time measurements. */
1835     mdoutf_tng_close(outf);
1836
1837     /* Stop measuring walltime */
1838     walltime_accounting_end(walltime_accounting);
1839
1840     if (bRerunMD && MASTER(cr))
1841     {
1842         close_trx(status);
1843     }
1844
1845     if (!(cr->duty & DUTY_PME))
1846     {
1847         /* Tell the PME only node to finish */
1848         gmx_pme_send_finish(cr);
1849     }
1850
1851     if (MASTER(cr))
1852     {
1853         if (ir->nstcalcenergy > 0 && !bRerunMD)
1854         {
1855             print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
1856                        eprAVER, mdebin, fcd, groups, &(ir->opts));
1857         }
1858     }
1859
1860     done_mdoutf(outf);
1861
1862     if (bPMETune)
1863     {
1864         pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
1865     }
1866
1867     done_shellfc(fplog, shellfc, step_rel);
1868
1869     if (useReplicaExchange && MASTER(cr))
1870     {
1871         print_replica_exchange_statistics(fplog, repl_ex);
1872     }
1873
1874     // Clean up swapcoords
1875     if (ir->eSwapCoords != eswapNO)
1876     {
1877         finish_swapcoords(ir->swap);
1878     }
1879
1880     /* Do essential dynamics cleanup if needed. Close .edo file */
1881     done_ed(&ed);
1882
1883     /* IMD cleanup, if bIMD is TRUE. */
1884     IMD_finalize(ir->bIMD, ir->imd);
1885
1886     walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1887     if (step_rel >= wcycle_get_reset_counters(wcycle) &&
1888         signals[eglsRESETCOUNTERS].set == 0 &&
1889         !bResetCountersHalfMaxH)
1890     {
1891         walltime_accounting_set_valid_finish(walltime_accounting);
1892     }
1893
1894     return 0;
1895 }