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