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