src/programs/mdrun/md.cpp

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