src/programs/mdrun/md.cpp

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