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