src/programs/mdrun/md.cpp

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