src/programs/mdrun/md.c

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