src/gromacs/mdrun/md.cpp

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