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