src/programs/mdrun/md.cpp

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