src/gromacs/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 /*! \internal \file
  38  *
  39  * \brief Implements the integrator for normal molecular dynamics simulations
  40  *
  41  * \author David van der Spoel <david.vanderspoel@icm.uu.se>
  42  * \ingroup module_mdrun
  43  */
  44 #include "gmxpre.h"
  45
  46 #include <cinttypes>
  47 #include <cmath>
  48 #include <cstdio>
  49 #include <cstdlib>
  50
  51 #include <algorithm>
  52 #include <memory>
  53
  54 #include "gromacs/awh/awh.h"
  55 #include "gromacs/commandline/filenm.h"
  56 #include "gromacs/compat/make_unique.h"
  57 #include "gromacs/domdec/collect.h"
  58 #include "gromacs/domdec/domdec.h"
  59 #include "gromacs/domdec/domdec_network.h"
  60 #include "gromacs/domdec/domdec_struct.h"
  61 #include "gromacs/domdec/partition.h"
  62 #include "gromacs/essentialdynamics/edsam.h"
  63 #include "gromacs/ewald/pme.h"
  64 #include "gromacs/ewald/pme-load-balancing.h"
  65 #include "gromacs/fileio/trxio.h"
  66 #include "gromacs/gmxlib/network.h"
  67 #include "gromacs/gmxlib/nrnb.h"
  68 #include "gromacs/gpu_utils/gpu_utils.h"
  69 #include "gromacs/imd/imd.h"
  70 #include "gromacs/listed-forces/manage-threading.h"
  71 #include "gromacs/math/functions.h"
  72 #include "gromacs/math/utilities.h"
  73 #include "gromacs/math/vec.h"
  74 #include "gromacs/math/vectypes.h"
  75 #include "gromacs/mdlib/checkpointhandler.h"
  76 #include "gromacs/mdlib/compute_io.h"
  77 #include "gromacs/mdlib/constr.h"
  78 #include "gromacs/mdlib/ebin.h"
  79 #include "gromacs/mdlib/expanded.h"
  80 #include "gromacs/mdlib/force.h"
  81 #include "gromacs/mdlib/force_flags.h"
  82 #include "gromacs/mdlib/forcerec.h"
  83 #include "gromacs/mdlib/md_support.h"
  84 #include "gromacs/mdlib/mdatoms.h"
  85 #include "gromacs/mdlib/mdebin.h"
  86 #include "gromacs/mdlib/mdoutf.h"
  87 #include "gromacs/mdlib/mdrun.h"
  88 #include "gromacs/mdlib/mdsetup.h"
  89 #include "gromacs/mdlib/membed.h"
  90 #include "gromacs/mdlib/nb_verlet.h"
  91 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
  92 #include "gromacs/mdlib/ns.h"
  93 #include "gromacs/mdlib/resethandler.h"
  94 #include "gromacs/mdlib/shellfc.h"
  95 #include "gromacs/mdlib/sighandler.h"
  96 #include "gromacs/mdlib/sim_util.h"
  97 #include "gromacs/mdlib/simulationsignal.h"
  98 #include "gromacs/mdlib/stophandler.h"
  99 #include "gromacs/mdlib/tgroup.h"
 100 #include "gromacs/mdlib/trajectory_writing.h"
 101 #include "gromacs/mdlib/update.h"
 102 #include "gromacs/mdlib/vcm.h"
 103 #include "gromacs/mdlib/vsite.h"
 104 #include "gromacs/mdtypes/awh-history.h"
 105 #include "gromacs/mdtypes/awh-params.h"
 106 #include "gromacs/mdtypes/commrec.h"
 107 #include "gromacs/mdtypes/df_history.h"
 108 #include "gromacs/mdtypes/energyhistory.h"
 109 #include "gromacs/mdtypes/fcdata.h"
 110 #include "gromacs/mdtypes/forcerec.h"
 111 #include "gromacs/mdtypes/group.h"
 112 #include "gromacs/mdtypes/inputrec.h"
 113 #include "gromacs/mdtypes/interaction_const.h"
 114 #include "gromacs/mdtypes/md_enums.h"
 115 #include "gromacs/mdtypes/mdatom.h"
 116 #include "gromacs/mdtypes/observableshistory.h"
 117 #include "gromacs/mdtypes/pullhistory.h"
 118 #include "gromacs/mdtypes/state.h"
 119 #include "gromacs/pbcutil/mshift.h"
 120 #include "gromacs/pbcutil/pbc.h"
 121 #include "gromacs/pulling/output.h"
 122 #include "gromacs/pulling/pull.h"
 123 #include "gromacs/swap/swapcoords.h"
 124 #include "gromacs/timing/wallcycle.h"
 125 #include "gromacs/timing/walltime_accounting.h"
 126 #include "gromacs/topology/atoms.h"
 127 #include "gromacs/topology/idef.h"
 128 #include "gromacs/topology/mtop_util.h"
 129 #include "gromacs/topology/topology.h"
 130 #include "gromacs/trajectory/trajectoryframe.h"
 131 #include "gromacs/utility/basedefinitions.h"
 132 #include "gromacs/utility/cstringutil.h"
 133 #include "gromacs/utility/fatalerror.h"
 134 #include "gromacs/utility/logger.h"
 135 #include "gromacs/utility/real.h"
 136 #include "gromacs/utility/smalloc.h"
 137
 138 #include "integrator.h"
 139 #include "replicaexchange.h"
 140
 141 #if GMX_FAHCORE
 142 #include "corewrap.h"
 143 #endif
 144
 145 using gmx::SimulationSignaller;
 146
 147 void gmx::Integrator::do_md()
 148 {
 149     // TODO Historically, the EM and MD "integrators" used different
 150     // names for the t_inputrec *parameter, but these must have the
 151     // same name, now that it's a member of a struct. We use this ir
 152     // alias to avoid a large ripple of nearly useless changes.
 153     // t_inputrec is being replaced by IMdpOptionsProvider, so this
 154     // will go away eventually.
 155     t_inputrec             *ir   = inputrec;
 156     gmx_mdoutf             *outf = nullptr;
 157     int64_t                 step, step_rel;
 158     double                  t, t0, lam0[efptNR];
 159     gmx_bool                bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
 160     gmx_bool                bNS, bNStList, bSimAnn, bStopCM,
 161                             bFirstStep, bInitStep, bLastStep = FALSE;
 162     gmx_bool                bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
 163     gmx_bool                do_ene, do_log, do_verbose;
 164     gmx_bool                bMasterState;
 165     int                     force_flags, cglo_flags;
 166     tensor                  force_vir, shake_vir, total_vir, tmp_vir, pres;
 167     int                     i, m;
 168     rvec                    mu_tot;
 169     t_vcm                  *vcm;
 170     matrix                  parrinellorahmanMu, M;
 171     gmx_repl_ex_t           repl_ex = nullptr;
 172     gmx_localtop_t         *top;
 173     t_mdebin               *mdebin   = nullptr;
 174     gmx_enerdata_t         *enerd;
 175     PaddedVector<gmx::RVec> f {};
 176     gmx_global_stat_t       gstat;
 177     gmx_update_t           *upd   = nullptr;
 178     t_graph                *graph = nullptr;
 179     gmx_groups_t           *groups;
 180     gmx_ekindata_t         *ekind;
 181     gmx_shellfc_t          *shellfc;
 182     gmx_bool                bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
 183     gmx_bool                bTemp, bPres, bTrotter;
 184     real                    dvdl_constr;
 185     rvec                   *cbuf        = nullptr;
 186     int                     cbuf_nalloc = 0;
 187     matrix                  lastbox;
 188     int                     lamnew  = 0;
 189     /* for FEP */
 190     int                     nstfep = 0;
 191     double                  cycles;
 192     real                    saved_conserved_quantity = 0;
 193     real                    last_ekin                = 0;
 194     t_extmass               MassQ;
 195     int                   **trotter_seq;
 196     char                    sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
 197
 198     /* PME load balancing data for GPU kernels */
 199     pme_load_balancing_t *pme_loadbal      = nullptr;
 200     gmx_bool              bPMETune         = FALSE;
 201     gmx_bool              bPMETunePrinting = FALSE;
 202
 203     /* Interactive MD */
 204     gmx_bool          bIMDstep = FALSE;
 205
 206     /* Domain decomposition could incorrectly miss a bonded
 207        interaction, but checking for that requires a global
 208        communication stage, which does not otherwise happen in DD
 209        code. So we do that alongside the first global energy reduction
 210        after a new DD is made. These variables handle whether the
 211        check happens, and the result it returns. */
 212     bool              shouldCheckNumberOfBondedInteractions = false;
 213     int               totalNumberOfBondedInteractions       = -1;
 214
 215     SimulationSignals signals;
 216     // Most global communnication stages don't propagate mdrun
 217     // signals, and will use this object to achieve that.
 218     SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
 219
 220     if (!mdrunOptions.writeConfout)
 221     {
 222         // This is on by default, and the main known use case for
 223         // turning it off is for convenience in benchmarking, which is
 224         // something that should not show up in the general user
 225         // interface.
 226         GMX_LOG(mdlog.info).asParagraph().
 227             appendText("The -noconfout functionality is deprecated, and may be removed in a future version.");
 228     }
 229
 230     /* md-vv uses averaged full step velocities for T-control
 231        md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
 232        md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
 233     bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
 234
 235     const bool bRerunMD      = false;
 236     int        nstglobalcomm = mdrunOptions.globalCommunicationInterval;
 237
 238     nstglobalcomm   = check_nstglobalcomm(mdlog, nstglobalcomm, ir);
 239     bGStatEveryStep = (nstglobalcomm == 1);
 240
 241     groups = &top_global->groups;
 242
 243     std::unique_ptr<EssentialDynamics> ed = nullptr;
 244     if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
 245     {
 246         /* Initialize essential dynamics sampling */
 247         ed = init_edsam(mdlog,
 248                         opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
 249                         top_global,
 250                         ir, cr, constr,
 251                         state_global, observablesHistory,
 252                         oenv, mdrunOptions.continuationOptions.appendFiles);
 253     }
 254
 255     /* Initial values */
 256     init_md(fplog, cr, outputProvider, ir, oenv, mdrunOptions,
 257             &t, &t0, state_global, lam0,
 258             nrnb, top_global, &upd, deform,
 259             nfile, fnm, &outf, &mdebin,
 260             force_vir, shake_vir, total_vir, pres, mu_tot, &bSimAnn, &vcm, wcycle);
 261
 262     /* Energy terms and groups */
 263     snew(enerd, 1);
 264     init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
 265                   enerd);
 266
 267     /* Kinetic energy data */
 268     snew(ekind, 1);
 269     init_ekindata(fplog, top_global, &(ir->opts), ekind);
 270     /* Copy the cos acceleration to the groups struct */
 271     ekind->cosacc.cos_accel = ir->cos_accel;
 272
 273     gstat = global_stat_init(ir);
 274
 275     /* Check for polarizable models and flexible constraints */
 276     shellfc = init_shell_flexcon(fplog,
 277                                  top_global, constr ? constr->numFlexibleConstraints() : 0,
 278                                  ir->nstcalcenergy, DOMAINDECOMP(cr));
 279
 280     {
 281         double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1);
 282         if ((io > 2000) && MASTER(cr))
 283         {
 284             fprintf(stderr,
 285                     "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
 286                     io);
 287         }
 288     }
 289
 290     /* Set up interactive MD (IMD) */
 291     init_IMD(ir, cr, ms, top_global, fplog, ir->nstcalcenergy,
 292              MASTER(cr) ? state_global->x.rvec_array() : nullptr,
 293              nfile, fnm, oenv, mdrunOptions);
 294
 295     // Local state only becomes valid now.
 296     std::unique_ptr<t_state> stateInstance;
 297     t_state *                state;
 298
 299     if (DOMAINDECOMP(cr))
 300     {
 301         top = dd_init_local_top(top_global);
 302
 303         stateInstance = compat::make_unique<t_state>();
 304         state         = stateInstance.get();
 305         dd_init_local_state(cr->dd, state_global, state);
 306
 307         /* Distribute the charge groups over the nodes from the master node */
 308         dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
 309                             state_global, top_global, ir,
 310                             state, &f, mdAtoms, top, fr,
 311                             vsite, constr,
 312                             nrnb, nullptr, FALSE);
 313         shouldCheckNumberOfBondedInteractions = true;
 314         update_realloc(upd, state->natoms);
 315     }
 316     else
 317     {
 318         state_change_natoms(state_global, state_global->natoms);
 319         f.resizeWithPadding(state_global->natoms);
 320         /* Copy the pointer to the global state */
 321         state = state_global;
 322
 323         snew(top, 1);
 324         mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
 325                                   &graph, mdAtoms, constr, vsite, shellfc);
 326
 327         update_realloc(upd, state->natoms);
 328     }
 329
 330     auto mdatoms = mdAtoms->mdatoms();
 331
 332     // NOTE: The global state is no longer used at this point.
 333     // But state_global is still used as temporary storage space for writing
 334     // the global state to file and potentially for replica exchange.
 335     // (Global topology should persist.)
 336
 337     update_mdatoms(mdatoms, state->lambda[efptMASS]);
 338
 339     const ContinuationOptions &continuationOptions    = mdrunOptions.continuationOptions;
 340     bool                       startingFromCheckpoint = continuationOptions.startedFromCheckpoint;
 341
 342     if (ir->bExpanded)
 343     {
 344         /* Check nstexpanded here, because the grompp check was broken */
 345         if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
 346         {
 347             gmx_fatal(FARGS, "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
 348         }
 349         init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
 350     }
 351
 352     if (MASTER(cr))
 353     {
 354         if (startingFromCheckpoint)
 355         {
 356             /* Update mdebin with energy history if appending to output files */
 357             if (continuationOptions.appendFiles)
 358             {
 359                 /* If no history is available (because a checkpoint is from before
 360                  * it was written) make a new one later, otherwise restore it.
 361                  */
 362                 if (observablesHistory->energyHistory)
 363                 {
 364                     restore_energyhistory_from_state(mdebin, observablesHistory->energyHistory.get());
 365                 }
 366             }
 367             else if (observablesHistory->energyHistory)
 368             {
 369                 /* We might have read an energy history from checkpoint.
 370                  * As we are not appending, we want to restart the statistics.
 371                  * Free the allocated memory and reset the counts.
 372                  */
 373                 observablesHistory->energyHistory = {};
 374                 /* We might have read a pull history from checkpoint.
 375                  * We will still want to keep the statistics, so that the files
 376                  * can be joined and still be meaningful.
 377                  * This means that observablesHistory->pullHistory
 378                  * should not be reset.
 379                  */
 380             }
 381         }
 382         if (!observablesHistory->energyHistory)
 383         {
 384             observablesHistory->energyHistory = compat::make_unique<energyhistory_t>();
 385         }
 386         if (!observablesHistory->pullHistory)
 387         {
 388             observablesHistory->pullHistory = compat::make_unique<PullHistory>();
 389         }
 390         /* Set the initial energy history in state by updating once */
 391         update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
 392     }
 393
 394     preparePrevStepPullCom(ir, mdatoms, state, state_global, cr, startingFromCheckpoint);
 395
 396     // TODO: Remove this by converting AWH into a ForceProvider
 397     auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingFromCheckpoint,
 398                                 shellfc != nullptr,
 399                                 opt2fn("-awh", nfile, fnm), ir->pull_work);
 400
 401     const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
 402     if (useReplicaExchange && MASTER(cr))
 403     {
 404         repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
 405                                         replExParams);
 406     }
 407     /* PME tuning is only supported in the Verlet scheme, with PME for
 408      * Coulomb. It is not supported with only LJ PME. */
 409     bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) &&
 410                 !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
 411     if (bPMETune)
 412     {
 413         pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box,
 414                          *fr->ic, *fr->nbv->listParams, fr->pmedata, use_GPU(fr->nbv),
 415                          &bPMETunePrinting);
 416     }
 417
 418     if (!ir->bContinuation)
 419     {
 420         if (state->flags & (1 << estV))
 421         {
 422             auto v = makeArrayRef(state->v);
 423             /* Set the velocities of vsites, shells and frozen atoms to zero */
 424             for (i = 0; i < mdatoms->homenr; i++)
 425             {
 426                 if (mdatoms->ptype[i] == eptVSite ||
 427                     mdatoms->ptype[i] == eptShell)
 428                 {
 429                     clear_rvec(v[i]);
 430                 }
 431                 else if (mdatoms->cFREEZE)
 432                 {
 433                     for (m = 0; m < DIM; m++)
 434                     {
 435                         if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
 436                         {
 437                             v[i][m] = 0;
 438                         }
 439                     }
 440                 }
 441             }
 442         }
 443
 444         if (constr)
 445         {
 446             /* Constrain the initial coordinates and velocities */
 447             do_constrain_first(fplog, constr, ir, mdatoms, state);
 448         }
 449         if (vsite)
 450         {
 451             /* Construct the virtual sites for the initial configuration */
 452             construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
 453                              top->idef.iparams, top->idef.il,
 454                              fr->ePBC, fr->bMolPBC, cr, state->box);
 455         }
 456     }
 457
 458     if (ir->efep != efepNO)
 459     {
 460         /* Set free energy calculation frequency as the greatest common
 461          * denominator of nstdhdl and repl_ex_nst. */
 462         nstfep = ir->fepvals->nstdhdl;
 463         if (ir->bExpanded)
 464         {
 465             nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
 466         }
 467         if (useReplicaExchange)
 468         {
 469             nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
 470         }
 471     }
 472
 473     /* Be REALLY careful about what flags you set here. You CANNOT assume
 474      * this is the first step, since we might be restarting from a checkpoint,
 475      * and in that case we should not do any modifications to the state.
 476      */
 477     bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
 478
 479     if (continuationOptions.haveReadEkin)
 480     {
 481         restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
 482     }
 483
 484     cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT
 485                   | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
 486                   | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
 487                   | (continuationOptions.haveReadEkin ? CGLO_READEKIN : 0));
 488
 489     bSumEkinhOld = FALSE;
 490     /* To minimize communication, compute_globals computes the COM velocity
 491      * and the kinetic energy for the velocities without COM motion removed.
 492      * Thus to get the kinetic energy without the COM contribution, we need
 493      * to call compute_globals twice.
 494      */
 495     for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
 496     {
 497         int cglo_flags_iteration = cglo_flags;
 498         if (bStopCM && cgloIteration == 0)
 499         {
 500             cglo_flags_iteration |= CGLO_STOPCM;
 501             cglo_flags_iteration &= ~CGLO_TEMPERATURE;
 502         }
 503         compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 504                         nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 505                         constr, &nullSignaller, state->box,
 506                         &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
 507                         | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
 508     }
 509     checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
 510                                     top_global, top, state,
 511                                     &shouldCheckNumberOfBondedInteractions);
 512     if (ir->eI == eiVVAK)
 513     {
 514         /* a second call to get the half step temperature initialized as well */
 515         /* we do the same call as above, but turn the pressure off -- internally to
 516            compute_globals, this is recognized as a velocity verlet half-step
 517            kinetic energy calculation.  This minimized excess variables, but
 518            perhaps loses some logic?*/
 519
 520         compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 521                         nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 522                         constr, &nullSignaller, state->box,
 523                         nullptr, &bSumEkinhOld,
 524                         cglo_flags & ~CGLO_PRESSURE);
 525     }
 526
 527     /* Calculate the initial half step temperature, and save the ekinh_old */
 528     if (!continuationOptions.startedFromCheckpoint)
 529     {
 530         for (i = 0; (i < ir->opts.ngtc); i++)
 531         {
 532             copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
 533         }
 534     }
 535
 536     /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
 537        temperature control */
 538     trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
 539
 540     if (MASTER(cr))
 541     {
 542         if (!ir->bContinuation)
 543         {
 544             if (constr && ir->eConstrAlg == econtLINCS)
 545             {
 546                 fprintf(fplog,
 547                         "RMS relative constraint deviation after constraining: %.2e\n",
 548                         constr->rmsd());
 549             }
 550             if (EI_STATE_VELOCITY(ir->eI))
 551             {
 552                 real temp = enerd->term[F_TEMP];
 553                 if (ir->eI != eiVV)
 554                 {
 555                     /* Result of Ekin averaged over velocities of -half
 556                      * and +half step, while we only have -half step here.
 557                      */
 558                     temp *= 2;
 559                 }
 560                 fprintf(fplog, "Initial temperature: %g K\n", temp);
 561             }
 562         }
 563
 564         char tbuf[20];
 565         fprintf(stderr, "starting mdrun '%s'\n",
 566                 *(top_global->name));
 567         if (ir->nsteps >= 0)
 568         {
 569             sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
 570         }
 571         else
 572         {
 573             sprintf(tbuf, "%s", "infinite");
 574         }
 575         if (ir->init_step > 0)
 576         {
 577             fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
 578                     gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
 579                     gmx_step_str(ir->init_step, sbuf2),
 580                     ir->init_step*ir->delta_t);
 581         }
 582         else
 583         {
 584             fprintf(stderr, "%s steps, %s ps.\n",
 585                     gmx_step_str(ir->nsteps, sbuf), tbuf);
 586         }
 587         fprintf(fplog, "\n");
 588     }
 589
 590     walltime_accounting_start_time(walltime_accounting);
 591     wallcycle_start(wcycle, ewcRUN);
 592     print_start(fplog, cr, walltime_accounting, "mdrun");
 593
 594 #if GMX_FAHCORE
 595     /* safest point to do file checkpointing is here.  More general point would be immediately before integrator call */
 596     int chkpt_ret = fcCheckPointParallel( cr->nodeid,
 597                                           NULL, 0);
 598     if (chkpt_ret == 0)
 599     {
 600         gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
 601     }
 602 #endif
 603
 604     /***********************************************************
 605      *
 606      *             Loop over MD steps
 607      *
 608      ************************************************************/
 609
 610     bFirstStep       = TRUE;
 611     /* Skip the first Nose-Hoover integration when we get the state from tpx */
 612     bInitStep        = !startingFromCheckpoint || EI_VV(ir->eI);
 613     bSumEkinhOld     = FALSE;
 614     bExchanged       = FALSE;
 615     bNeedRepartition = FALSE;
 616
 617     bool simulationsShareState = false;
 618     int  nstSignalComm         = nstglobalcomm;
 619     {
 620         // TODO This implementation of ensemble orientation restraints is nasty because
 621         // a user can't just do multi-sim with single-sim orientation restraints.
 622         bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
 623         bool awhUsesMultiSim         = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
 624
 625         // Replica exchange, ensemble restraints and AWH need all
 626         // simulations to remain synchronized, so they need
 627         // checkpoints and stop conditions to act on the same step, so
 628         // the propagation of such signals must take place between
 629         // simulations, not just within simulations.
 630         // TODO: Make algorithm initializers set these flags.
 631         simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
 632
 633         if (simulationsShareState)
 634         {
 635             // Inter-simulation signal communication does not need to happen
 636             // often, so we use a minimum of 200 steps to reduce overhead.
 637             const int c_minimumInterSimulationSignallingInterval = 200;
 638             nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
 639         }
 640     }
 641
 642     auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
 643                 compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
 644                 MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
 645                 mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
 646
 647     auto checkpointHandler = compat::make_unique<CheckpointHandler>(
 648                 compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]),
 649                 simulationsShareState, ir->nstlist == 0, MASTER(cr),
 650                 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
 651
 652     const bool resetCountersIsLocal = true;
 653     auto       resetHandler         = compat::make_unique<ResetHandler>(
 654                 compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]), !resetCountersIsLocal,
 655                 ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
 656                 mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
 657
 658     DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion   = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
 659     DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion  = (DOMAINDECOMP(cr) ? DdCloseBalanceRegionAfterForceComputation::yes : DdCloseBalanceRegionAfterForceComputation::no);
 660
 661     step     = ir->init_step;
 662     step_rel = 0;
 663
 664     // TODO extract this to new multi-simulation module
 665     if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
 666     {
 667         if (!multisim_int_all_are_equal(ms, ir->nsteps))
 668         {
 669             GMX_LOG(mdlog.warning).appendText(
 670                     "Note: The number of steps is not consistent across multi simulations,\n"
 671                     "but we are proceeding anyway!");
 672         }
 673         if (!multisim_int_all_are_equal(ms, ir->init_step))
 674         {
 675             GMX_LOG(mdlog.warning).appendText(
 676                     "Note: The initial step is not consistent across multi simulations,\n"
 677                     "but we are proceeding anyway!");
 678         }
 679     }
 680
 681     /* and stop now if we should */
 682     bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
 683     while (!bLastStep)
 684     {
 685
 686         /* Determine if this is a neighbor search step */
 687         bNStList = (ir->nstlist > 0  && step % ir->nstlist == 0);
 688
 689         if (bPMETune && bNStList)
 690         {
 691             /* PME grid + cut-off optimization with GPUs or PME nodes */
 692             pme_loadbal_do(pme_loadbal, cr,
 693                            (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
 694                            fplog, mdlog,
 695                            *ir, fr, *state,
 696                            wcycle,
 697                            step, step_rel,
 698                            &bPMETunePrinting);
 699         }
 700
 701         wallcycle_start(wcycle, ewcSTEP);
 702
 703         bLastStep = (step_rel == ir->nsteps);
 704         t         = t0 + step*ir->delta_t;
 705
 706         // TODO Refactor this, so that nstfep does not need a default value of zero
 707         if (ir->efep != efepNO || ir->bSimTemp)
 708         {
 709             /* find and set the current lambdas */
 710             setCurrentLambdasLocal(step, ir->fepvals, lam0, state);
 711
 712             bDoDHDL      = do_per_step(step, ir->fepvals->nstdhdl);
 713             bDoFEP       = ((ir->efep != efepNO) && do_per_step(step, nstfep));
 714             bDoExpanded  = (do_per_step(step, ir->expandedvals->nstexpanded)
 715                             && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
 716         }
 717
 718         bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
 719                      do_per_step(step, replExParams.exchangeInterval));
 720
 721         if (bSimAnn)
 722         {
 723             update_annealing_target_temp(ir, t, upd);
 724         }
 725
 726         /* Stop Center of Mass motion */
 727         bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
 728
 729         /* Determine whether or not to do Neighbour Searching */
 730         bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
 731
 732         bLastStep = bLastStep || stopHandler->stoppingAfterCurrentStep(bNS);
 733
 734         /* do_log triggers energy and virial calculation. Because this leads
 735          * to different code paths, forces can be different. Thus for exact
 736          * continuation we should avoid extra log output.
 737          * Note that the || bLastStep can result in non-exact continuation
 738          * beyond the last step. But we don't consider that to be an issue.
 739          */
 740         do_log     = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep;
 741         do_verbose = mdrunOptions.verbose &&
 742             (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
 743
 744         if (bNS && !(bFirstStep && ir->bContinuation))
 745         {
 746             bMasterState = FALSE;
 747             /* Correct the new box if it is too skewed */
 748             if (inputrecDynamicBox(ir))
 749             {
 750                 if (correct_box(fplog, step, state->box, graph))
 751                 {
 752                     bMasterState = TRUE;
 753                 }
 754             }
 755             if (DOMAINDECOMP(cr) && bMasterState)
 756             {
 757                 dd_collect_state(cr->dd, state, state_global);
 758             }
 759
 760             if (DOMAINDECOMP(cr))
 761             {
 762                 /* Repartition the domain decomposition */
 763                 dd_partition_system(fplog, mdlog, step, cr,
 764                                     bMasterState, nstglobalcomm,
 765                                     state_global, top_global, ir,
 766                                     state, &f, mdAtoms, top, fr,
 767                                     vsite, constr,
 768                                     nrnb, wcycle,
 769                                     do_verbose && !bPMETunePrinting);
 770                 shouldCheckNumberOfBondedInteractions = true;
 771                 update_realloc(upd, state->natoms);
 772             }
 773         }
 774
 775         if (MASTER(cr) && do_log)
 776         {
 777             print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
 778         }
 779
 780         if (ir->efep != efepNO)
 781         {
 782             update_mdatoms(mdatoms, state->lambda[efptMASS]);
 783         }
 784
 785         if (bExchanged)
 786         {
 787
 788             /* We need the kinetic energy at minus the half step for determining
 789              * the full step kinetic energy and possibly for T-coupling.*/
 790             /* This may not be quite working correctly yet . . . . */
 791             compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 792                             wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
 793                             constr, &nullSignaller, state->box,
 794                             &totalNumberOfBondedInteractions, &bSumEkinhOld,
 795                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
 796             checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
 797                                             top_global, top, state,
 798                                             &shouldCheckNumberOfBondedInteractions);
 799         }
 800         clear_mat(force_vir);
 801
 802         checkpointHandler->decideIfCheckpointingThisStep(bNS, bFirstStep, bLastStep);
 803
 804         /* Determine the energy and pressure:
 805          * at nstcalcenergy steps and at energy output steps (set below).
 806          */
 807         if (EI_VV(ir->eI) && (!bInitStep))
 808         {
 809             bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
 810             bCalcVir      = bCalcEnerStep ||
 811                 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
 812         }
 813         else
 814         {
 815             bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
 816             bCalcVir      = bCalcEnerStep ||
 817                 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
 818         }
 819         bCalcEner = bCalcEnerStep;
 820
 821         do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
 822
 823         if (do_ene || do_log || bDoReplEx)
 824         {
 825             bCalcVir  = TRUE;
 826             bCalcEner = TRUE;
 827         }
 828
 829         /* Do we need global communication ? */
 830         bGStat = (bCalcVir || bCalcEner || bStopCM ||
 831                   do_per_step(step, nstglobalcomm) ||
 832                   (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
 833
 834         force_flags = (GMX_FORCE_STATECHANGED |
 835                        ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0) |
 836                        GMX_FORCE_ALLFORCES |
 837                        (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
 838                        (bCalcEner ? GMX_FORCE_ENERGY : 0) |
 839                        (bDoFEP ? GMX_FORCE_DHDL : 0)
 840                        );
 841
 842         if (shellfc)
 843         {
 844             /* Now is the time to relax the shells */
 845             relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
 846                                 enforcedRotation, step,
 847                                 ir, bNS, force_flags, top,
 848                                 constr, enerd, fcd,
 849                                 state, f.arrayRefWithPadding(), force_vir, mdatoms,
 850                                 nrnb, wcycle, graph, groups,
 851                                 shellfc, fr, ppForceWorkload, t, mu_tot,
 852                                 vsite,
 853                                 ddOpenBalanceRegion, ddCloseBalanceRegion);
 854         }
 855         else
 856         {
 857             /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
 858                (or the AWH update will be performed twice for one step when continuing). It would be best to
 859                call this update function from do_md_trajectory_writing but that would occur after do_force.
 860                One would have to divide the update_awh function into one function applying the AWH force
 861                and one doing the AWH bias update. The update AWH bias function could then be called after
 862                do_md_trajectory_writing (then containing update_awh_history).
 863                The checkpointing will in the future probably moved to the start of the md loop which will
 864                rid of this issue. */
 865             if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
 866             {
 867                 awh->updateHistory(state_global->awhHistory.get());
 868             }
 869
 870             /* The coordinates (x) are shifted (to get whole molecules)
 871              * in do_force.
 872              * This is parallellized as well, and does communication too.
 873              * Check comments in sim_util.c
 874              */
 875             do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation,
 876                      step, nrnb, wcycle, top, groups,
 877                      state->box, state->x.arrayRefWithPadding(), &state->hist,
 878                      f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
 879                      state->lambda, graph,
 880                      fr, ppForceWorkload, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
 881                      (bNS ? GMX_FORCE_NS : 0) | force_flags,
 882                      ddOpenBalanceRegion, ddCloseBalanceRegion);
 883         }
 884
 885         if (EI_VV(ir->eI) && !startingFromCheckpoint)
 886         /*  ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
 887         {
 888             rvec *vbuf = nullptr;
 889
 890             wallcycle_start(wcycle, ewcUPDATE);
 891             if (ir->eI == eiVV && bInitStep)
 892             {
 893                 /* if using velocity verlet with full time step Ekin,
 894                  * take the first half step only to compute the
 895                  * virial for the first step. From there,
 896                  * revert back to the initial coordinates
 897                  * so that the input is actually the initial step.
 898                  */
 899                 snew(vbuf, state->natoms);
 900                 copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
 901             }
 902             else
 903             {
 904                 /* this is for NHC in the Ekin(t+dt/2) version of vv */
 905                 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
 906             }
 907
 908             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
 909                           ekind, M, upd, etrtVELOCITY1,
 910                           cr, constr);
 911
 912             wallcycle_stop(wcycle, ewcUPDATE);
 913             constrain_velocities(step, nullptr,
 914                                  state,
 915                                  shake_vir,
 916                                  constr,
 917                                  bCalcVir, do_log, do_ene);
 918             wallcycle_start(wcycle, ewcUPDATE);
 919             /* if VV, compute the pressure and constraints */
 920             /* For VV2, we strictly only need this if using pressure
 921              * control, but we really would like to have accurate pressures
 922              * printed out.
 923              * Think about ways around this in the future?
 924              * For now, keep this choice in comments.
 925              */
 926             /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
 927             /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
 928             bPres = TRUE;
 929             bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
 930             if (bCalcEner && ir->eI == eiVVAK)
 931             {
 932                 bSumEkinhOld = TRUE;
 933             }
 934             /* for vv, the first half of the integration actually corresponds to the previous step.
 935                So we need information from the last step in the first half of the integration */
 936             if (bGStat || do_per_step(step-1, nstglobalcomm))
 937             {
 938                 wallcycle_stop(wcycle, ewcUPDATE);
 939                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 940                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
 941                                 constr, &nullSignaller, state->box,
 942                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
 943                                 (bGStat ? CGLO_GSTAT : 0)
 944                                 | (bCalcEner ? CGLO_ENERGY : 0)
 945                                 | (bTemp ? CGLO_TEMPERATURE : 0)
 946                                 | (bPres ? CGLO_PRESSURE : 0)
 947                                 | (bPres ? CGLO_CONSTRAINT : 0)
 948                                 | (bStopCM ? CGLO_STOPCM : 0)
 949                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
 950                                 | CGLO_SCALEEKIN
 951                                 );
 952                 /* explanation of above:
 953                    a) We compute Ekin at the full time step
 954                    if 1) we are using the AveVel Ekin, and it's not the
 955                    initial step, or 2) if we are using AveEkin, but need the full
 956                    time step kinetic energy for the pressure (always true now, since we want accurate statistics).
 957                    b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
 958                    EkinAveVel because it's needed for the pressure */
 959                 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
 960                                                 top_global, top, state,
 961                                                 &shouldCheckNumberOfBondedInteractions);
 962                 wallcycle_start(wcycle, ewcUPDATE);
 963             }
 964             /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
 965             if (!bInitStep)
 966             {
 967                 if (bTrotter)
 968                 {
 969                     m_add(force_vir, shake_vir, total_vir);     /* we need the un-dispersion corrected total vir here */
 970                     trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
 971
 972                     /* TODO This is only needed when we're about to write
 973                      * a checkpoint, because we use it after the restart
 974                      * (in a kludge?). But what should we be doing if
 975                      * startingFromCheckpoint or bInitStep are true? */
 976                     if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
 977                     {
 978                         copy_mat(shake_vir, state->svir_prev);
 979                         copy_mat(force_vir, state->fvir_prev);
 980                     }
 981                     if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
 982                     {
 983                         /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
 984                         enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
 985                         enerd->term[F_EKIN] = trace(ekind->ekin);
 986                     }
 987                 }
 988                 else if (bExchanged)
 989                 {
 990                     wallcycle_stop(wcycle, ewcUPDATE);
 991                     /* We need the kinetic energy at minus the half step for determining
 992                      * the full step kinetic energy and possibly for T-coupling.*/
 993                     /* This may not be quite working correctly yet . . . . */
 994                     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
 995                                     wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
 996                                     constr, &nullSignaller, state->box,
 997                                     nullptr, &bSumEkinhOld,
 998                                     CGLO_GSTAT | CGLO_TEMPERATURE);
 999                     wallcycle_start(wcycle, ewcUPDATE);
1000                 }
1001             }
1002             /* if it's the initial step, we performed this first step just to get the constraint virial */
1003             if (ir->eI == eiVV && bInitStep)
1004             {
1005                 copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
1006                 sfree(vbuf);
1007             }
1008             wallcycle_stop(wcycle, ewcUPDATE);
1009         }
1010
1011         /* compute the conserved quantity */
1012         if (EI_VV(ir->eI))
1013         {
1014             saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1015             if (ir->eI == eiVV)
1016             {
1017                 last_ekin = enerd->term[F_EKIN];
1018             }
1019             if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1020             {
1021                 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1022             }
1023             /* sum up the foreign energy and dhdl terms for vv.  currently done every step so that dhdl is correct in the .edr */
1024             if (ir->efep != efepNO)
1025             {
1026                 sum_dhdl(enerd, state->lambda, ir->fepvals);
1027             }
1028         }
1029
1030         /* ########  END FIRST UPDATE STEP  ############## */
1031         /* ########  If doing VV, we now have v(dt) ###### */
1032         if (bDoExpanded)
1033         {
1034             /* perform extended ensemble sampling in lambda - we don't
1035                actually move to the new state before outputting
1036                statistics, but if performing simulated tempering, we
1037                do update the velocities and the tau_t. */
1038
1039             lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v.rvec_array(), mdatoms);
1040             /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1041             if (MASTER(cr))
1042             {
1043                 copy_df_history(state_global->dfhist, state->dfhist);
1044             }
1045         }
1046
1047         /* Now we have the energies and forces corresponding to the
1048          * coordinates at time t. We must output all of this before
1049          * the update.
1050          */
1051         do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1052                                  ir, state, state_global, observablesHistory,
1053                                  top_global, fr,
1054                                  outf, mdebin, ekind, f,
1055                                  checkpointHandler->isCheckpointingStep(),
1056                                  bRerunMD, bLastStep,
1057                                  mdrunOptions.writeConfout,
1058                                  bSumEkinhOld);
1059         /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1060         bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x.rvec_array(), ir, t, wcycle);
1061
1062         /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1063         if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1064         {
1065             copy_mat(state->svir_prev, shake_vir);
1066             copy_mat(state->fvir_prev, force_vir);
1067         }
1068
1069         stopHandler->setSignal();
1070         resetHandler->setSignal(walltime_accounting);
1071
1072         if (bGStat || !PAR(cr))
1073         {
1074             /* In parallel we only have to check for checkpointing in steps
1075              * where we do global communication,
1076              *  otherwise the other nodes don't know.
1077              */
1078             checkpointHandler->setSignal(walltime_accounting);
1079         }
1080
1081         /* #########   START SECOND UPDATE STEP ################# */
1082
1083         /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1084            in preprocessing */
1085
1086         if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1087         {
1088             gmx_bool bIfRandomize;
1089             bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, upd, constr);
1090             /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1091             if (constr && bIfRandomize)
1092             {
1093                 constrain_velocities(step, nullptr,
1094                                      state,
1095                                      tmp_vir,
1096                                      constr,
1097                                      bCalcVir, do_log, do_ene);
1098             }
1099         }
1100         /* Box is changed in update() when we do pressure coupling,
1101          * but we should still use the old box for energy corrections and when
1102          * writing it to the energy file, so it matches the trajectory files for
1103          * the same timestep above. Make a copy in a separate array.
1104          */
1105         copy_mat(state->box, lastbox);
1106
1107         dvdl_constr = 0;
1108
1109         wallcycle_start(wcycle, ewcUPDATE);
1110         /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1111         if (bTrotter)
1112         {
1113             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1114             /* We can only do Berendsen coupling after we have summed
1115              * the kinetic energy or virial. Since the happens
1116              * in global_state after update, we should only do it at
1117              * step % nstlist = 1 with bGStatEveryStep=FALSE.
1118              */
1119         }
1120         else
1121         {
1122             update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1123             update_pcouple_before_coordinates(fplog, step, ir, state,
1124                                               parrinellorahmanMu, M,
1125                                               bInitStep);
1126         }
1127
1128         if (EI_VV(ir->eI))
1129         {
1130             /* velocity half-step update */
1131             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1132                           ekind, M, upd, etrtVELOCITY2,
1133                           cr, constr);
1134         }
1135
1136         /* Above, initialize just copies ekinh into ekin,
1137          * it doesn't copy position (for VV),
1138          * and entire integrator for MD.
1139          */
1140
1141         if (ir->eI == eiVVAK)
1142         {
1143             /* We probably only need md->homenr, not state->natoms */
1144             if (state->natoms > cbuf_nalloc)
1145             {
1146                 cbuf_nalloc = state->natoms;
1147                 srenew(cbuf, cbuf_nalloc);
1148             }
1149             copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1150         }
1151
1152         update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1153                       ekind, M, upd, etrtPOSITION, cr, constr);
1154         wallcycle_stop(wcycle, ewcUPDATE);
1155
1156         constrain_coordinates(step, &dvdl_constr, state,
1157                               shake_vir,
1158                               upd, constr,
1159                               bCalcVir, do_log, do_ene);
1160         update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
1161                               cr, nrnb, wcycle, upd, constr, do_log, do_ene);
1162         finish_update(ir, mdatoms,
1163                       state, graph,
1164                       nrnb, wcycle, upd, constr);
1165
1166         if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1167         {
1168             updatePrevStepPullCom(ir->pull_work, state);
1169         }
1170
1171         if (ir->eI == eiVVAK)
1172         {
1173             /* erase F_EKIN and F_TEMP here? */
1174             /* just compute the kinetic energy at the half step to perform a trotter step */
1175             compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1176                             wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1177                             constr, &nullSignaller, lastbox,
1178                             nullptr, &bSumEkinhOld,
1179                             (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1180                             );
1181             wallcycle_start(wcycle, ewcUPDATE);
1182             trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1183             /* now we know the scaling, we can compute the positions again again */
1184             copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1185
1186             update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1187                           ekind, M, upd, etrtPOSITION, cr, constr);
1188             wallcycle_stop(wcycle, ewcUPDATE);
1189
1190             /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1191             /* are the small terms in the shake_vir here due
1192              * to numerical errors, or are they important
1193              * physically? I'm thinking they are just errors, but not completely sure.
1194              * For now, will call without actually constraining, constr=NULL*/
1195             finish_update(ir, mdatoms,
1196                           state, graph,
1197                           nrnb, wcycle, upd, nullptr);
1198         }
1199         if (EI_VV(ir->eI))
1200         {
1201             /* this factor or 2 correction is necessary
1202                because half of the constraint force is removed
1203                in the vv step, so we have to double it.  See
1204                the Redmine issue #1255.  It is not yet clear
1205                if the factor of 2 is exact, or just a very
1206                good approximation, and this will be
1207                investigated.  The next step is to see if this
1208                can be done adding a dhdl contribution from the
1209                rattle step, but this is somewhat more
1210                complicated with the current code. Will be
1211                investigated, hopefully for 4.6.3. However,
1212                this current solution is much better than
1213                having it completely wrong.
1214              */
1215             enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1216         }
1217         else
1218         {
1219             enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1220         }
1221
1222         if (vsite != nullptr)
1223         {
1224             wallcycle_start(wcycle, ewcVSITECONSTR);
1225             if (graph != nullptr)
1226             {
1227                 shift_self(graph, state->box, state->x.rvec_array());
1228             }
1229             construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
1230                              top->idef.iparams, top->idef.il,
1231                              fr->ePBC, fr->bMolPBC, cr, state->box);
1232
1233             if (graph != nullptr)
1234             {
1235                 unshift_self(graph, state->box, state->x.rvec_array());
1236             }
1237             wallcycle_stop(wcycle, ewcVSITECONSTR);
1238         }
1239
1240         /* ############## IF NOT VV, Calculate globals HERE  ############ */
1241         /* With Leap-Frog we can skip compute_globals at
1242          * non-communication steps, but we need to calculate
1243          * the kinetic energy one step before communication.
1244          */
1245         {
1246             // Organize to do inter-simulation signalling on steps if
1247             // and when algorithms require it.
1248             bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1249
1250             if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1251             {
1252                 // Since we're already communicating at this step, we
1253                 // can propagate intra-simulation signals. Note that
1254                 // check_nstglobalcomm has the responsibility for
1255                 // choosing the value of nstglobalcomm that is one way
1256                 // bGStat becomes true, so we can't get into a
1257                 // situation where e.g. checkpointing can't be
1258                 // signalled.
1259                 bool                doIntraSimSignal = true;
1260                 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1261
1262                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1263                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1264                                 constr, &signaller,
1265                                 lastbox,
1266                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1267                                 (bGStat ? CGLO_GSTAT : 0)
1268                                 | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1269                                 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1270                                 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1271                                 | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
1272                                 | CGLO_CONSTRAINT
1273                                 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1274                                 );
1275                 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1276                                                 top_global, top, state,
1277                                                 &shouldCheckNumberOfBondedInteractions);
1278             }
1279         }
1280
1281         /* #############  END CALC EKIN AND PRESSURE ################# */
1282
1283         /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1284            the virial that should probably be addressed eventually. state->veta has better properies,
1285            but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1286            generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
1287
1288         if (ir->efep != efepNO && !EI_VV(ir->eI))
1289         {
1290             /* Sum up the foreign energy and dhdl terms for md and sd.
1291                Currently done every step so that dhdl is correct in the .edr */
1292             sum_dhdl(enerd, state->lambda, ir->fepvals);
1293         }
1294
1295         update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1296                                          pres, force_vir, shake_vir,
1297                                          parrinellorahmanMu,
1298                                          state, nrnb, upd);
1299
1300         /* ################# END UPDATE STEP 2 ################# */
1301         /* #### We now have r(t+dt) and v(t+dt/2)  ############# */
1302
1303         /* The coordinates (x) were unshifted in update */
1304         if (!bGStat)
1305         {
1306             /* We will not sum ekinh_old,
1307              * so signal that we still have to do it.
1308              */
1309             bSumEkinhOld = TRUE;
1310         }
1311
1312         if (bCalcEner)
1313         {
1314             /* #########  BEGIN PREPARING EDR OUTPUT  ###########  */
1315
1316             /* use the directly determined last velocity, not actually the averaged half steps */
1317             if (bTrotter && ir->eI == eiVV)
1318             {
1319                 enerd->term[F_EKIN] = last_ekin;
1320             }
1321             enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1322
1323             if (integratorHasConservedEnergyQuantity(ir))
1324             {
1325                 if (EI_VV(ir->eI))
1326                 {
1327                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1328                 }
1329                 else
1330                 {
1331                     enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1332                 }
1333             }
1334             /* #########  END PREPARING EDR OUTPUT  ###########  */
1335         }
1336
1337         /* Output stuff */
1338         if (MASTER(cr))
1339         {
1340             if (fplog && do_log && bDoExpanded)
1341             {
1342                 /* only needed if doing expanded ensemble */
1343                 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1344                                           state_global->dfhist, state->fep_state, ir->nstlog, step);
1345             }
1346             if (bCalcEner)
1347             {
1348                 upd_mdebin(mdebin, bDoDHDL, bCalcEnerStep,
1349                            t, mdatoms->tmass, enerd, state,
1350                            ir->fepvals, ir->expandedvals, lastbox,
1351                            shake_vir, force_vir, total_vir, pres,
1352                            ekind, mu_tot, constr);
1353             }
1354             else
1355             {
1356                 upd_mdebin_step(mdebin);
1357             }
1358
1359             gmx_bool do_dr  = do_per_step(step, ir->nstdisreout);
1360             gmx_bool do_or  = do_per_step(step, ir->nstorireout);
1361
1362             print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
1363                        step, t,
1364                        eprNORMAL, mdebin, fcd, groups, &(ir->opts), awh.get());
1365
1366             if (ir->bPull)
1367             {
1368                 pull_print_output(ir->pull_work, step, t);
1369             }
1370
1371             if (do_per_step(step, ir->nstlog))
1372             {
1373                 if (fflush(fplog) != 0)
1374                 {
1375                     gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1376                 }
1377             }
1378         }
1379         if (bDoExpanded)
1380         {
1381             /* Have to do this part _after_ outputting the logfile and the edr file */
1382             /* Gets written into the state at the beginning of next loop*/
1383             state->fep_state = lamnew;
1384         }
1385         /* Print the remaining wall clock time for the run */
1386         if (isMasterSimMasterRank(ms, cr) &&
1387             (do_verbose || gmx_got_usr_signal()) &&
1388             !bPMETunePrinting)
1389         {
1390             if (shellfc)
1391             {
1392                 fprintf(stderr, "\n");
1393             }
1394             print_time(stderr, walltime_accounting, step, ir, cr);
1395         }
1396
1397         /* Ion/water position swapping.
1398          * Not done in last step since trajectory writing happens before this call
1399          * in the MD loop and exchanges would be lost anyway. */
1400         bNeedRepartition = FALSE;
1401         if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1402             do_per_step(step, ir->swap->nstswap))
1403         {
1404             bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
1405                                              as_rvec_array(state->x.data()),
1406                                              state->box,
1407                                              MASTER(cr) && mdrunOptions.verbose,
1408                                              bRerunMD);
1409
1410             if (bNeedRepartition && DOMAINDECOMP(cr))
1411             {
1412                 dd_collect_state(cr->dd, state, state_global);
1413             }
1414         }
1415
1416         /* Replica exchange */
1417         bExchanged = FALSE;
1418         if (bDoReplEx)
1419         {
1420             bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
1421                                           state_global, enerd,
1422                                           state, step, t);
1423         }
1424
1425         if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1426         {
1427             dd_partition_system(fplog, mdlog, step, cr, TRUE, 1,
1428                                 state_global, top_global, ir,
1429                                 state, &f, mdAtoms, top, fr,
1430                                 vsite, constr,
1431                                 nrnb, wcycle, FALSE);
1432             shouldCheckNumberOfBondedInteractions = true;
1433             update_realloc(upd, state->natoms);
1434         }
1435
1436         bFirstStep             = FALSE;
1437         bInitStep              = FALSE;
1438         startingFromCheckpoint = false;
1439
1440         /* #######  SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1441         /* With all integrators, except VV, we need to retain the pressure
1442          * at the current step for coupling at the next step.
1443          */
1444         if ((state->flags & (1<<estPRES_PREV)) &&
1445             (bGStatEveryStep ||
1446              (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1447         {
1448             /* Store the pressure in t_state for pressure coupling
1449              * at the next MD step.
1450              */
1451             copy_mat(pres, state->pres_prev);
1452         }
1453
1454         /* #######  END SET VARIABLES FOR NEXT ITERATION ###### */
1455
1456         if ( (membed != nullptr) && (!bLastStep) )
1457         {
1458             rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1459         }
1460
1461         cycles = wallcycle_stop(wcycle, ewcSTEP);
1462         if (DOMAINDECOMP(cr) && wcycle)
1463         {
1464             dd_cycles_add(cr->dd, cycles, ddCyclStep);
1465         }
1466
1467         /* increase the MD step number */
1468         step++;
1469         step_rel++;
1470
1471         resetHandler->resetCounters(
1472                 step, step_rel, mdlog, fplog, cr, (use_GPU(fr->nbv) ? fr->nbv : nullptr),
1473                 nrnb, fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
1474
1475         /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1476         IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
1477
1478     }
1479     /* End of main MD loop */
1480
1481     /* Closing TNG files can include compressing data. Therefore it is good to do that
1482      * before stopping the time measurements. */
1483     mdoutf_tng_close(outf);
1484
1485     /* Stop measuring walltime */
1486     walltime_accounting_end_time(walltime_accounting);
1487
1488     if (!thisRankHasDuty(cr, DUTY_PME))
1489     {
1490         /* Tell the PME only node to finish */
1491         gmx_pme_send_finish(cr);
1492     }
1493
1494     if (MASTER(cr))
1495     {
1496         if (ir->nstcalcenergy > 0)
1497         {
1498             print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
1499                        eprAVER, mdebin, fcd, groups, &(ir->opts), awh.get());
1500         }
1501     }
1502     done_mdebin(mdebin);
1503     done_mdoutf(outf);
1504
1505     if (bPMETune)
1506     {
1507         pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
1508     }
1509
1510     done_shellfc(fplog, shellfc, step_rel);
1511
1512     if (useReplicaExchange && MASTER(cr))
1513     {
1514         print_replica_exchange_statistics(fplog, repl_ex);
1515     }
1516
1517     // Clean up swapcoords
1518     if (ir->eSwapCoords != eswapNO)
1519     {
1520         finish_swapcoords(ir->swap);
1521     }
1522
1523     /* IMD cleanup, if bIMD is TRUE. */
1524     IMD_finalize(ir->bIMD, ir->imd);
1525
1526     walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1527
1528     destroy_enerdata(enerd);
1529
1530     sfree(enerd);
1531
1532     /* Clean up topology. top->atomtypes has an allocated pointer if no domain decomposition*/
1533     if (!DOMAINDECOMP(cr))
1534     {
1535         done_atomtypes(&top->atomtypes);
1536     }
1537     sfree(top);
1538 }