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