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