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