src/programs/mdrun/md.cpp

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