src/programs/mdrun/md.cpp

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