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