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