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