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