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