2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
39 * \brief Implements the integrator for normal molecular dynamics simulations
41 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
42 * \ingroup module_mdrun
55 #include "gromacs/awh/awh.h"
56 #include "gromacs/commandline/filenm.h"
57 #include "gromacs/compat/make_unique.h"
58 #include "gromacs/domdec/collect.h"
59 #include "gromacs/domdec/domdec.h"
60 #include "gromacs/domdec/domdec_network.h"
61 #include "gromacs/domdec/domdec_struct.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/compute_io.h"
76 #include "gromacs/mdlib/constr.h"
77 #include "gromacs/mdlib/ebin.h"
78 #include "gromacs/mdlib/expanded.h"
79 #include "gromacs/mdlib/force.h"
80 #include "gromacs/mdlib/force_flags.h"
81 #include "gromacs/mdlib/forcerec.h"
82 #include "gromacs/mdlib/md_support.h"
83 #include "gromacs/mdlib/mdatoms.h"
84 #include "gromacs/mdlib/mdebin.h"
85 #include "gromacs/mdlib/mdoutf.h"
86 #include "gromacs/mdlib/mdrun.h"
87 #include "gromacs/mdlib/mdsetup.h"
88 #include "gromacs/mdlib/membed.h"
89 #include "gromacs/mdlib/nb_verlet.h"
90 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
91 #include "gromacs/mdlib/ns.h"
92 #include "gromacs/mdlib/shellfc.h"
93 #include "gromacs/mdlib/sighandler.h"
94 #include "gromacs/mdlib/sim_util.h"
95 #include "gromacs/mdlib/simulationsignal.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/mdtypes/awh-history.h"
102 #include "gromacs/mdtypes/awh-params.h"
103 #include "gromacs/mdtypes/commrec.h"
104 #include "gromacs/mdtypes/df_history.h"
105 #include "gromacs/mdtypes/energyhistory.h"
106 #include "gromacs/mdtypes/fcdata.h"
107 #include "gromacs/mdtypes/forcerec.h"
108 #include "gromacs/mdtypes/group.h"
109 #include "gromacs/mdtypes/inputrec.h"
110 #include "gromacs/mdtypes/interaction_const.h"
111 #include "gromacs/mdtypes/md_enums.h"
112 #include "gromacs/mdtypes/mdatom.h"
113 #include "gromacs/mdtypes/observableshistory.h"
114 #include "gromacs/mdtypes/state.h"
115 #include "gromacs/pbcutil/mshift.h"
116 #include "gromacs/pbcutil/pbc.h"
117 #include "gromacs/pulling/pull.h"
118 #include "gromacs/swap/swapcoords.h"
119 #include "gromacs/timing/wallcycle.h"
120 #include "gromacs/timing/walltime_accounting.h"
121 #include "gromacs/topology/atoms.h"
122 #include "gromacs/topology/idef.h"
123 #include "gromacs/topology/mtop_util.h"
124 #include "gromacs/topology/topology.h"
125 #include "gromacs/trajectory/trajectoryframe.h"
126 #include "gromacs/utility/basedefinitions.h"
127 #include "gromacs/utility/cstringutil.h"
128 #include "gromacs/utility/fatalerror.h"
129 #include "gromacs/utility/logger.h"
130 #include "gromacs/utility/real.h"
131 #include "gromacs/utility/smalloc.h"
133 #include "integrator.h"
134 #include "replicaexchange.h"
137 #include "corewrap.h"
140 using gmx::SimulationSignaller;
142 //! Resets all the counters.
143 static void reset_all_counters(FILE *fplog, const gmx::MDLogger &mdlog, t_commrec *cr,
145 gmx_wallcycle_t wcycle, t_nrnb *nrnb,
146 gmx_walltime_accounting_t walltime_accounting,
147 struct nonbonded_verlet_t *nbv,
148 struct gmx_pme_t *pme)
150 char sbuf[STEPSTRSIZE];
152 /* Reset all the counters related to performance over the run */
153 GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
154 "step %s: resetting all time and cycle counters",
155 gmx_step_str(step, sbuf));
159 nbnxn_gpu_reset_timings(nbv);
162 if (pme_gpu_task_enabled(pme))
164 pme_gpu_reset_timings(pme);
167 if (use_GPU(nbv) || pme_gpu_task_enabled(pme))
172 wallcycle_stop(wcycle, ewcRUN);
173 wallcycle_reset_all(wcycle);
174 if (DOMAINDECOMP(cr))
176 reset_dd_statistics_counters(cr->dd);
179 wallcycle_start(wcycle, ewcRUN);
180 walltime_accounting_reset_time(walltime_accounting, step);
181 print_date_and_time(fplog, cr->nodeid, "Restarted time", gmx_gettime());
184 /*! \brief Copy the state from \p rerunFrame to \p globalState and, if requested, construct vsites
186 * \param[in] rerunFrame The trajectory frame to compute energy/forces for
187 * \param[in,out] globalState The global state container
188 * \param[in] constructVsites When true, vsite coordinates are constructed
189 * \param[in] vsite Vsite setup, can be nullptr when \p constructVsites = false
190 * \param[in] idef Topology parameters, used for constructing vsites
191 * \param[in] timeStep Time step, used for constructing vsites
192 * \param[in] forceRec Force record, used for constructing vsites
193 * \param[in,out] graph The molecular graph, used for constructing vsites when != nullptr
194 * \param[in,out] warnWhenNoV When true, issue a warning when no velocities are present in \p rerunFrame; is set to false when a warning was issued
196 static void prepareRerunState(const t_trxframe &rerunFrame,
197 t_state *globalState,
198 bool constructVsites,
199 const gmx_vsite_t *vsite,
202 const t_forcerec &forceRec,
204 gmx_bool *warnWhenNoV)
206 for (int i = 0; i < globalState->natoms; i++)
208 copy_rvec(rerunFrame.x[i], globalState->x[i]);
212 for (int i = 0; i < globalState->natoms; i++)
214 copy_rvec(rerunFrame.v[i], globalState->v[i]);
219 for (int i = 0; i < globalState->natoms; i++)
221 clear_rvec(globalState->v[i]);
225 fprintf(stderr, "\nWARNING: Some frames do not contain velocities.\n"
226 " Ekin, temperature and pressure are incorrect,\n"
227 " the virial will be incorrect when constraints are present.\n"
229 *warnWhenNoV = FALSE;
232 copy_mat(rerunFrame.box, globalState->box);
236 GMX_ASSERT(vsite, "Need valid vsite for constructing vsites");
240 /* Following is necessary because the graph may get out of sync
241 * with the coordinates if we only have every N'th coordinate set
243 mk_mshift(nullptr, graph, forceRec.ePBC, globalState->box, as_rvec_array(globalState->x.data()));
244 shift_self(graph, globalState->box, as_rvec_array(globalState->x.data()));
246 construct_vsites(vsite, as_rvec_array(globalState->x.data()), timeStep, as_rvec_array(globalState->v.data()),
247 idef.iparams, idef.il,
248 forceRec.ePBC, forceRec.bMolPBC, nullptr, globalState->box);
251 unshift_self(graph, globalState->box, as_rvec_array(globalState->x.data()));
256 void gmx::Integrator::do_md()
258 // TODO Historically, the EM and MD "integrators" used different
259 // names for the t_inputrec *parameter, but these must have the
260 // same name, now that it's a member of a struct. We use this ir
261 // alias to avoid a large ripple of nearly useless changes.
262 // t_inputrec is being replaced by IMdpOptionsProvider, so this
263 // will go away eventually.
264 t_inputrec *ir = inputrec;
265 gmx_mdoutf *outf = nullptr;
266 int64_t step, step_rel;
267 double t, t0, lam0[efptNR];
268 gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
269 gmx_bool bNS, bNStList, bSimAnn, bStopCM,
270 bFirstStep, bInitStep, bLastStep = FALSE;
271 gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
272 gmx_bool do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE,
273 bForceUpdate = FALSE, bCPT;
274 gmx_bool bMasterState;
275 int force_flags, cglo_flags;
276 tensor force_vir, shake_vir, total_vir, tmp_vir, pres;
281 matrix parrinellorahmanMu, M;
283 gmx_repl_ex_t repl_ex = nullptr;
286 t_mdebin *mdebin = nullptr;
287 gmx_enerdata_t *enerd;
288 PaddedRVecVector f {};
289 gmx_global_stat_t gstat;
290 gmx_update_t *upd = nullptr;
291 t_graph *graph = nullptr;
292 gmx_groups_t *groups;
293 gmx_ekindata_t *ekind;
294 gmx_shellfc_t *shellfc;
295 gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
296 gmx_bool bResetCountersHalfMaxH = FALSE;
297 gmx_bool bTemp, bPres, bTrotter;
299 rvec *cbuf = nullptr;
306 real saved_conserved_quantity = 0;
310 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
311 int handled_stop_condition = gmx_stop_cond_none; /* compare to get_stop_condition*/
314 /* PME load balancing data for GPU kernels */
315 pme_load_balancing_t *pme_loadbal = nullptr;
316 gmx_bool bPMETune = FALSE;
317 gmx_bool bPMETunePrinting = FALSE;
320 gmx_bool bIMDstep = FALSE;
323 /* Temporary addition for FAHCORE checkpointing */
326 /* Domain decomposition could incorrectly miss a bonded
327 interaction, but checking for that requires a global
328 communication stage, which does not otherwise happen in DD
329 code. So we do that alongside the first global energy reduction
330 after a new DD is made. These variables handle whether the
331 check happens, and the result it returns. */
332 bool shouldCheckNumberOfBondedInteractions = false;
333 int totalNumberOfBondedInteractions = -1;
335 SimulationSignals signals;
336 // Most global communnication stages don't propagate mdrun
337 // signals, and will use this object to achieve that.
338 SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
340 if (!mdrunOptions.writeConfout)
342 // This is on by default, and the main known use case for
343 // turning it off is for convenience in benchmarking, which is
344 // something that should not show up in the general user
346 GMX_LOG(mdlog.info).asParagraph().
347 appendText("The -noconfout functionality is deprecated, and may be removed in a future version.");
349 if (mdrunOptions.timingOptions.resetHalfway)
351 GMX_LOG(mdlog.info).asParagraph().
352 appendText("The -resethway functionality is deprecated, and may be removed in a future version.");
355 /* Signal to reset the counters half the simulation steps. */
356 wcycle_set_reset_counters(wcycle, ir->nsteps/2);
358 /* Signal to reset the counters halfway the simulation time. */
359 bResetCountersHalfMaxH = (mdrunOptions.maximumHoursToRun > 0);
362 /* md-vv uses averaged full step velocities for T-control
363 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
364 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
365 bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
367 const gmx_bool bRerunMD = mdrunOptions.rerun;
368 int nstglobalcomm = mdrunOptions.globalCommunicationInterval;
372 /* Since we don't know if the frames read are related in any way,
373 * rebuild the neighborlist at every step.
376 ir->nstcalcenergy = 1;
380 nstglobalcomm = check_nstglobalcomm(mdlog, nstglobalcomm, ir);
381 bGStatEveryStep = (nstglobalcomm == 1);
385 ir->nstxout_compressed = 0;
387 groups = &top_global->groups;
389 std::unique_ptr<EssentialDynamics> ed = nullptr;
390 if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
392 /* Initialize essential dynamics sampling */
393 ed = init_edsam(opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
396 state_global, observablesHistory,
397 oenv, mdrunOptions.continuationOptions.appendFiles);
401 init_md(fplog, cr, outputProvider, ir, oenv, mdrunOptions,
402 &t, &t0, state_global, lam0,
403 nrnb, top_global, &upd, deform,
404 nfile, fnm, &outf, &mdebin,
405 force_vir, shake_vir, total_vir, pres, mu_tot, &bSimAnn, &vcm, wcycle);
407 /* Energy terms and groups */
409 init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
412 /* Kinetic energy data */
414 init_ekindata(fplog, top_global, &(ir->opts), ekind);
415 /* Copy the cos acceleration to the groups struct */
416 ekind->cosacc.cos_accel = ir->cos_accel;
418 gstat = global_stat_init(ir);
420 /* Check for polarizable models and flexible constraints */
421 shellfc = init_shell_flexcon(fplog,
422 top_global, constr ? constr->numFlexibleConstraints() : 0,
423 ir->nstcalcenergy, DOMAINDECOMP(cr));
426 double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1);
427 if ((io > 2000) && MASTER(cr))
430 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
435 /* Set up interactive MD (IMD) */
436 init_IMD(ir, cr, ms, top_global, fplog, ir->nstcalcenergy,
437 MASTER(cr) ? as_rvec_array(state_global->x.data()) : nullptr,
438 nfile, fnm, oenv, mdrunOptions);
440 // Local state only becomes valid now.
441 std::unique_ptr<t_state> stateInstance;
444 if (DOMAINDECOMP(cr))
446 top = dd_init_local_top(top_global);
448 stateInstance = compat::make_unique<t_state>();
449 state = stateInstance.get();
450 dd_init_local_state(cr->dd, state_global, state);
452 /* Distribute the charge groups over the nodes from the master node */
453 dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
454 state_global, top_global, ir,
455 state, &f, mdAtoms, top, fr,
457 nrnb, nullptr, FALSE);
458 shouldCheckNumberOfBondedInteractions = true;
459 update_realloc(upd, state->natoms);
463 state_change_natoms(state_global, state_global->natoms);
464 /* We need to allocate one element extra, since we might use
465 * (unaligned) 4-wide SIMD loads to access rvec entries.
467 f.resize(gmx::paddedRVecVectorSize(state_global->natoms));
468 /* Copy the pointer to the global state */
469 state = state_global;
472 mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
473 &graph, mdAtoms, constr, vsite, shellfc);
475 update_realloc(upd, state->natoms);
478 auto mdatoms = mdAtoms->mdatoms();
480 // NOTE: The global state is no longer used at this point.
481 // But state_global is still used as temporary storage space for writing
482 // the global state to file and potentially for replica exchange.
483 // (Global topology should persist.)
485 update_mdatoms(mdatoms, state->lambda[efptMASS]);
487 const ContinuationOptions &continuationOptions = mdrunOptions.continuationOptions;
488 bool startingFromCheckpoint = continuationOptions.startedFromCheckpoint;
492 init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
497 if (startingFromCheckpoint)
499 /* Update mdebin with energy history if appending to output files */
500 if (continuationOptions.appendFiles)
502 restore_energyhistory_from_state(mdebin, observablesHistory->energyHistory.get());
504 else if (observablesHistory->energyHistory != nullptr)
506 /* We might have read an energy history from checkpoint.
507 * As we are not appending, we want to restart the statistics.
508 * Free the allocated memory and reset the counts.
510 observablesHistory->energyHistory = {};
513 if (observablesHistory->energyHistory == nullptr)
515 observablesHistory->energyHistory = compat::make_unique<energyhistory_t>();
517 /* Set the initial energy history in state by updating once */
518 update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
521 // TODO: Remove this by converting AWH into a ForceProvider
522 auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingFromCheckpoint,
524 opt2fn("-awh", nfile, fnm), ir->pull_work);
526 const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
527 if (useReplicaExchange && MASTER(cr))
529 repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
532 /* PME tuning is only supported in the Verlet scheme, with PME for
533 * Coulomb. It is not supported with only LJ PME, or for
535 bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) && !bRerunMD &&
536 !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
539 pme_loadbal_init(&pme_loadbal, cr, mdlog, ir, state->box,
540 fr->ic, fr->nbv->listParams.get(), fr->pmedata, use_GPU(fr->nbv),
544 if (!ir->bContinuation && !bRerunMD)
546 if (state->flags & (1 << estV))
548 /* Set the velocities of vsites, shells and frozen atoms to zero */
549 for (i = 0; i < mdatoms->homenr; i++)
551 if (mdatoms->ptype[i] == eptVSite ||
552 mdatoms->ptype[i] == eptShell)
554 clear_rvec(state->v[i]);
556 else if (mdatoms->cFREEZE)
558 for (m = 0; m < DIM; m++)
560 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
571 /* Constrain the initial coordinates and velocities */
572 do_constrain_first(fplog, constr, ir, mdatoms, state);
576 /* Construct the virtual sites for the initial configuration */
577 construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, nullptr,
578 top->idef.iparams, top->idef.il,
579 fr->ePBC, fr->bMolPBC, cr, state->box);
583 if (ir->efep != efepNO)
585 /* Set free energy calculation frequency as the greatest common
586 * denominator of nstdhdl and repl_ex_nst. */
587 nstfep = ir->fepvals->nstdhdl;
590 nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
592 if (useReplicaExchange)
594 nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
598 /* Be REALLY careful about what flags you set here. You CANNOT assume
599 * this is the first step, since we might be restarting from a checkpoint,
600 * and in that case we should not do any modifications to the state.
602 bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
604 if (continuationOptions.haveReadEkin)
606 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
609 cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT
610 | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
611 | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
612 | (continuationOptions.haveReadEkin ? CGLO_READEKIN : 0));
614 bSumEkinhOld = FALSE;
615 /* To minimize communication, compute_globals computes the COM velocity
616 * and the kinetic energy for the velocities without COM motion removed.
617 * Thus to get the kinetic energy without the COM contribution, we need
618 * to call compute_globals twice.
620 for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
622 int cglo_flags_iteration = cglo_flags;
623 if (bStopCM && cgloIteration == 0)
625 cglo_flags_iteration |= CGLO_STOPCM;
626 cglo_flags_iteration &= ~CGLO_TEMPERATURE;
628 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
629 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
630 constr, &nullSignaller, state->box,
631 &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
632 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
634 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
635 top_global, top, state,
636 &shouldCheckNumberOfBondedInteractions);
637 if (ir->eI == eiVVAK)
639 /* a second call to get the half step temperature initialized as well */
640 /* we do the same call as above, but turn the pressure off -- internally to
641 compute_globals, this is recognized as a velocity verlet half-step
642 kinetic energy calculation. This minimized excess variables, but
643 perhaps loses some logic?*/
645 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
646 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
647 constr, &nullSignaller, state->box,
648 nullptr, &bSumEkinhOld,
649 cglo_flags & ~CGLO_PRESSURE);
652 /* Calculate the initial half step temperature, and save the ekinh_old */
653 if (!continuationOptions.startedFromCheckpoint)
655 for (i = 0; (i < ir->opts.ngtc); i++)
657 copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
661 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
662 temperature control */
663 trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
667 if (!ir->bContinuation)
669 if (constr && ir->eConstrAlg == econtLINCS)
672 "RMS relative constraint deviation after constraining: %.2e\n",
675 if (EI_STATE_VELOCITY(ir->eI))
677 real temp = enerd->term[F_TEMP];
680 /* Result of Ekin averaged over velocities of -half
681 * and +half step, while we only have -half step here.
685 fprintf(fplog, "Initial temperature: %g K\n", temp);
691 fprintf(stderr, "starting md rerun '%s', reading coordinates from"
692 " input trajectory '%s'\n\n",
693 *(top_global->name), opt2fn("-rerun", nfile, fnm));
694 if (mdrunOptions.verbose)
696 fprintf(stderr, "Calculated time to finish depends on nsteps from "
697 "run input file,\nwhich may not correspond to the time "
698 "needed to process input trajectory.\n\n");
704 fprintf(stderr, "starting mdrun '%s'\n",
705 *(top_global->name));
708 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
712 sprintf(tbuf, "%s", "infinite");
714 if (ir->init_step > 0)
716 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
717 gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
718 gmx_step_str(ir->init_step, sbuf2),
719 ir->init_step*ir->delta_t);
723 fprintf(stderr, "%s steps, %s ps.\n",
724 gmx_step_str(ir->nsteps, sbuf), tbuf);
727 fprintf(fplog, "\n");
730 walltime_accounting_start_time(walltime_accounting);
731 wallcycle_start(wcycle, ewcRUN);
732 print_start(fplog, cr, walltime_accounting, "mdrun");
734 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
736 chkpt_ret = fcCheckPointParallel( cr->nodeid,
740 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
744 /***********************************************************
748 ************************************************************/
750 /* if rerunMD then read coordinates and velocities from input trajectory */
753 if (getenv("GMX_FORCE_UPDATE"))
761 bLastStep = !read_first_frame(oenv, &status,
762 opt2fn("-rerun", nfile, fnm),
763 &rerun_fr, TRX_NEED_X | TRX_READ_V);
764 if (rerun_fr.natoms != top_global->natoms)
767 "Number of atoms in trajectory (%d) does not match the "
768 "run input file (%d)\n",
769 rerun_fr.natoms, top_global->natoms);
771 if (ir->ePBC != epbcNONE)
775 gmx_fatal(FARGS, "Rerun trajectory frame step %ld time %f does not contain a box, while pbc is used", rerun_fr.step, rerun_fr.time);
777 if (max_cutoff2(ir->ePBC, rerun_fr.box) < gmx::square(fr->rlist))
779 gmx_fatal(FARGS, "Rerun trajectory frame step %ld time %f has too small box dimensions", rerun_fr.step, rerun_fr.time);
786 rerun_parallel_comm(cr, &rerun_fr, &bLastStep);
789 if (ir->ePBC != epbcNONE)
791 /* Set the shift vectors.
792 * Necessary here when have a static box different from the tpr box.
794 calc_shifts(rerun_fr.box, fr->shift_vec);
799 /* Skip the first Nose-Hoover integration when we get the state from tpx */
800 bInitStep = !startingFromCheckpoint || EI_VV(ir->eI);
801 bSumEkinhOld = FALSE;
803 bNeedRepartition = FALSE;
805 bool simulationsShareState = false;
806 int nstSignalComm = nstglobalcomm;
808 // TODO This implementation of ensemble orientation restraints is nasty because
809 // a user can't just do multi-sim with single-sim orientation restraints.
810 bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
811 bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
813 // Replica exchange, ensemble restraints and AWH need all
814 // simulations to remain synchronized, so they need
815 // checkpoints and stop conditions to act on the same step, so
816 // the propagation of such signals must take place between
817 // simulations, not just within simulations.
818 // TODO: Make algorithm initializers set these flags.
819 simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
820 bool resetCountersIsLocal = true;
821 signals[eglsCHKPT] = SimulationSignal(!simulationsShareState);
822 signals[eglsSTOPCOND] = SimulationSignal(!simulationsShareState);
823 signals[eglsRESETCOUNTERS] = SimulationSignal(resetCountersIsLocal);
825 if (simulationsShareState)
827 // Inter-simulation signal communication does not need to happen
828 // often, so we use a minimum of 200 steps to reduce overhead.
829 const int c_minimumInterSimulationSignallingInterval = 200;
830 nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
834 DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
835 DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion = (DOMAINDECOMP(cr) ? DdCloseBalanceRegionAfterForceComputation::yes : DdCloseBalanceRegionAfterForceComputation::no);
837 step = ir->init_step;
840 // TODO extract this to new multi-simulation module
841 if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
843 if (!multisim_int_all_are_equal(ms, ir->nsteps))
845 GMX_LOG(mdlog.warning).appendText(
846 "Note: The number of steps is not consistent across multi simulations,\n"
847 "but we are proceeding anyway!");
849 if (!multisim_int_all_are_equal(ms, ir->init_step))
851 GMX_LOG(mdlog.warning).appendText(
852 "Note: The initial step is not consistent across multi simulations,\n"
853 "but we are proceeding anyway!");
857 /* and stop now if we should */
858 bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
862 /* Determine if this is a neighbor search step */
863 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
865 if (bPMETune && bNStList)
867 /* PME grid + cut-off optimization with GPUs or PME nodes */
868 pme_loadbal_do(pme_loadbal, cr,
869 (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
877 wallcycle_start(wcycle, ewcSTEP);
883 step = rerun_fr.step;
884 step_rel = step - ir->init_step;
897 bLastStep = (step_rel == ir->nsteps);
898 t = t0 + step*ir->delta_t;
901 // TODO Refactor this, so that nstfep does not need a default value of zero
902 if (ir->efep != efepNO || ir->bSimTemp)
904 /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
905 requiring different logic. */
910 setCurrentLambdasRerun(step, ir->fepvals, &rerun_fr, lam0, state_global);
915 setCurrentLambdasLocal(step, ir->fepvals, lam0, state);
917 bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
918 bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
919 bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
920 && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
923 bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
924 do_per_step(step, replExParams.exchangeInterval));
928 update_annealing_target_temp(ir, t, upd);
931 if (bRerunMD && MASTER(cr))
933 const bool constructVsites = ((vsite != nullptr) && mdrunOptions.rerunConstructVsites);
934 if (constructVsites && DOMAINDECOMP(cr))
936 gmx_fatal(FARGS, "Vsite recalculation with -rerun is not implemented with domain decomposition, use a single rank");
938 prepareRerunState(rerun_fr, state_global, constructVsites, vsite, top->idef, ir->delta_t, *fr, graph, &bRerunWarnNoV);
941 /* Stop Center of Mass motion */
942 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
946 /* for rerun MD always do Neighbour Searching */
947 bNS = (bFirstStep || ir->nstlist != 0);
951 /* Determine whether or not to do Neighbour Searching */
952 bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
955 /* < 0 means stop at next step, > 0 means stop at next NS step */
956 if ( (signals[eglsSTOPCOND].set < 0) ||
957 ( (signals[eglsSTOPCOND].set > 0 ) && ( bNS || ir->nstlist == 0)))
962 /* do_log triggers energy and virial calculation. Because this leads
963 * to different code paths, forces can be different. Thus for exact
964 * continuation we should avoid extra log output.
965 * Note that the || bLastStep can result in non-exact continuation
966 * beyond the last step. But we don't consider that to be an issue.
968 do_log = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep || bRerunMD;
969 do_verbose = mdrunOptions.verbose &&
970 (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep || bRerunMD);
972 if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD))
980 bMasterState = FALSE;
981 /* Correct the new box if it is too skewed */
982 if (inputrecDynamicBox(ir))
984 if (correct_box(fplog, step, state->box, graph))
989 if (DOMAINDECOMP(cr) && bMasterState)
991 dd_collect_state(cr->dd, state, state_global);
995 if (DOMAINDECOMP(cr))
997 /* Repartition the domain decomposition */
998 dd_partition_system(fplog, mdlog, step, cr,
999 bMasterState, nstglobalcomm,
1000 state_global, top_global, ir,
1001 state, &f, mdAtoms, top, fr,
1004 do_verbose && !bPMETunePrinting);
1005 shouldCheckNumberOfBondedInteractions = true;
1006 update_realloc(upd, state->natoms);
1010 if (MASTER(cr) && do_log)
1012 print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
1015 if (ir->efep != efepNO)
1017 update_mdatoms(mdatoms, state->lambda[efptMASS]);
1020 if ((bRerunMD && rerun_fr.bV) || bExchanged)
1023 /* We need the kinetic energy at minus the half step for determining
1024 * the full step kinetic energy and possibly for T-coupling.*/
1025 /* This may not be quite working correctly yet . . . . */
1026 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1027 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1028 constr, &nullSignaller, state->box,
1029 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1030 CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
1031 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1032 top_global, top, state,
1033 &shouldCheckNumberOfBondedInteractions);
1035 clear_mat(force_vir);
1037 /* We write a checkpoint at this MD step when:
1038 * either at an NS step when we signalled through gs,
1039 * or at the last step (but not when we do not want confout),
1040 * but never at the first step or with rerun.
1042 bCPT = ((((signals[eglsCHKPT].set != 0) && (bNS || ir->nstlist == 0)) ||
1043 (bLastStep && mdrunOptions.writeConfout)) &&
1044 step > ir->init_step && !bRerunMD);
1047 signals[eglsCHKPT].set = 0;
1050 /* Determine the energy and pressure:
1051 * at nstcalcenergy steps and at energy output steps (set below).
1053 if (EI_VV(ir->eI) && (!bInitStep))
1055 /* for vv, the first half of the integration actually corresponds
1056 to the previous step. bCalcEner is only required to be evaluated on the 'next' step,
1057 but the virial needs to be calculated on both the current step and the 'next' step. Future
1058 reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
1060 /* TODO: This is probably not what we want, we will write to energy file one step after nstcalcenergy steps. */
1061 bCalcEnerStep = do_per_step(step - 1, ir->nstcalcenergy);
1062 bCalcVir = bCalcEnerStep ||
1063 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
1067 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
1068 bCalcVir = bCalcEnerStep ||
1069 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
1071 bCalcEner = bCalcEnerStep;
1073 do_ene = (do_per_step(step, ir->nstenergy) || bLastStep || bRerunMD);
1075 if (do_ene || do_log || bDoReplEx)
1081 /* Do we need global communication ? */
1082 bGStat = (bCalcVir || bCalcEner || bStopCM ||
1083 do_per_step(step, nstglobalcomm) ||
1084 (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
1086 force_flags = (GMX_FORCE_STATECHANGED |
1087 ((inputrecDynamicBox(ir) || bRerunMD) ? GMX_FORCE_DYNAMICBOX : 0) |
1088 GMX_FORCE_ALLFORCES |
1089 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
1090 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
1091 (bDoFEP ? GMX_FORCE_DHDL : 0)
1096 /* Now is the time to relax the shells */
1097 relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
1098 enforcedRotation, step,
1099 ir, bNS, force_flags, top,
1101 state, f, force_vir, mdatoms,
1102 nrnb, wcycle, graph, groups,
1103 shellfc, fr, t, mu_tot,
1105 ddOpenBalanceRegion, ddCloseBalanceRegion);
1109 /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
1110 (or the AWH update will be performed twice for one step when continuing). It would be best to
1111 call this update function from do_md_trajectory_writing but that would occur after do_force.
1112 One would have to divide the update_awh function into one function applying the AWH force
1113 and one doing the AWH bias update. The update AWH bias function could then be called after
1114 do_md_trajectory_writing (then containing update_awh_history).
1115 The checkpointing will in the future probably moved to the start of the md loop which will
1116 rid of this issue. */
1117 if (awh && bCPT && MASTER(cr))
1119 awh->updateHistory(state_global->awhHistory.get());
1122 /* The coordinates (x) are shifted (to get whole molecules)
1124 * This is parallellized as well, and does communication too.
1125 * Check comments in sim_util.c
1127 do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation,
1128 step, nrnb, wcycle, top, groups,
1129 state->box, state->x, &state->hist,
1130 f, force_vir, mdatoms, enerd, fcd,
1131 state->lambda, graph,
1132 fr, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
1133 (bNS ? GMX_FORCE_NS : 0) | force_flags,
1134 ddOpenBalanceRegion, ddCloseBalanceRegion);
1137 if (EI_VV(ir->eI) && !startingFromCheckpoint && !bRerunMD)
1138 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1140 rvec *vbuf = nullptr;
1142 wallcycle_start(wcycle, ewcUPDATE);
1143 if (ir->eI == eiVV && bInitStep)
1145 /* if using velocity verlet with full time step Ekin,
1146 * take the first half step only to compute the
1147 * virial for the first step. From there,
1148 * revert back to the initial coordinates
1149 * so that the input is actually the initial step.
1151 snew(vbuf, state->natoms);
1152 copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
1156 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1157 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
1160 update_coords(step, ir, mdatoms, state, f, fcd,
1161 ekind, M, upd, etrtVELOCITY1,
1164 if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
1166 wallcycle_stop(wcycle, ewcUPDATE);
1167 constrain_velocities(step, nullptr,
1171 bCalcVir, do_log, do_ene);
1172 wallcycle_start(wcycle, ewcUPDATE);
1176 /* Need to unshift here if a do_force has been
1177 called in the previous step */
1178 unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1180 /* if VV, compute the pressure and constraints */
1181 /* For VV2, we strictly only need this if using pressure
1182 * control, but we really would like to have accurate pressures
1184 * Think about ways around this in the future?
1185 * For now, keep this choice in comments.
1187 /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
1188 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
1190 bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
1191 if (bCalcEner && ir->eI == eiVVAK)
1193 bSumEkinhOld = TRUE;
1195 /* for vv, the first half of the integration actually corresponds to the previous step.
1196 So we need information from the last step in the first half of the integration */
1197 if (bGStat || do_per_step(step-1, nstglobalcomm))
1199 wallcycle_stop(wcycle, ewcUPDATE);
1200 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1201 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1202 constr, &nullSignaller, state->box,
1203 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1204 (bGStat ? CGLO_GSTAT : 0)
1206 | (bTemp ? CGLO_TEMPERATURE : 0)
1207 | (bPres ? CGLO_PRESSURE : 0)
1208 | (bPres ? CGLO_CONSTRAINT : 0)
1209 | (bStopCM ? CGLO_STOPCM : 0)
1210 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1213 /* explanation of above:
1214 a) We compute Ekin at the full time step
1215 if 1) we are using the AveVel Ekin, and it's not the
1216 initial step, or 2) if we are using AveEkin, but need the full
1217 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1218 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1219 EkinAveVel because it's needed for the pressure */
1220 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1221 top_global, top, state,
1222 &shouldCheckNumberOfBondedInteractions);
1223 wallcycle_start(wcycle, ewcUPDATE);
1225 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1230 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
1231 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
1233 /* TODO This is only needed when we're about to write
1234 * a checkpoint, because we use it after the restart
1235 * (in a kludge?). But what should we be doing if
1236 * startingFromCheckpoint or bInitStep are true? */
1237 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1239 copy_mat(shake_vir, state->svir_prev);
1240 copy_mat(force_vir, state->fvir_prev);
1242 if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
1244 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1245 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
1246 enerd->term[F_EKIN] = trace(ekind->ekin);
1249 else if (bExchanged)
1251 wallcycle_stop(wcycle, ewcUPDATE);
1252 /* We need the kinetic energy at minus the half step for determining
1253 * the full step kinetic energy and possibly for T-coupling.*/
1254 /* This may not be quite working correctly yet . . . . */
1255 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1256 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1257 constr, &nullSignaller, state->box,
1258 nullptr, &bSumEkinhOld,
1259 CGLO_GSTAT | CGLO_TEMPERATURE);
1260 wallcycle_start(wcycle, ewcUPDATE);
1263 /* if it's the initial step, we performed this first step just to get the constraint virial */
1264 if (ir->eI == eiVV && bInitStep)
1266 copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
1269 wallcycle_stop(wcycle, ewcUPDATE);
1272 /* compute the conserved quantity */
1275 saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1278 last_ekin = enerd->term[F_EKIN];
1280 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1282 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1284 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1285 if (ir->efep != efepNO && !bRerunMD)
1287 sum_dhdl(enerd, state->lambda, ir->fepvals);
1291 /* ######## END FIRST UPDATE STEP ############## */
1292 /* ######## If doing VV, we now have v(dt) ###### */
1295 /* perform extended ensemble sampling in lambda - we don't
1296 actually move to the new state before outputting
1297 statistics, but if performing simulated tempering, we
1298 do update the velocities and the tau_t. */
1300 lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
1301 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1304 copy_df_history(state_global->dfhist, state->dfhist);
1308 /* Now we have the energies and forces corresponding to the
1309 * coordinates at time t. We must output all of this before
1312 do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1313 ir, state, state_global, observablesHistory,
1315 outf, mdebin, ekind, f,
1317 bCPT, bRerunMD, bLastStep,
1318 mdrunOptions.writeConfout,
1320 /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1321 bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
1323 /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1324 if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1326 copy_mat(state->svir_prev, shake_vir);
1327 copy_mat(state->fvir_prev, force_vir);
1330 double secondsSinceStart = walltime_accounting_get_time_since_start(walltime_accounting);
1332 /* Check whether everything is still allright */
1333 if ((static_cast<int>(gmx_get_stop_condition()) > handled_stop_condition)
1339 int nsteps_stop = -1;
1341 /* this just makes signals[].sig compatible with the hack
1342 of sending signals around by MPI_Reduce together with
1344 if ((gmx_get_stop_condition() == gmx_stop_cond_next_ns) ||
1345 (mdrunOptions.reproducible &&
1346 gmx_get_stop_condition() == gmx_stop_cond_next))
1348 /* We need at least two global communication steps to pass
1349 * around the signal. We stop at a pair-list creation step
1350 * to allow for exact continuation, when possible.
1352 signals[eglsSTOPCOND].sig = 1;
1353 nsteps_stop = std::max(ir->nstlist, 2*nstSignalComm);
1355 else if (gmx_get_stop_condition() == gmx_stop_cond_next)
1357 /* Stop directly after the next global communication step.
1358 * This breaks exact continuation.
1360 signals[eglsSTOPCOND].sig = -1;
1361 nsteps_stop = nstSignalComm + 1;
1366 "\n\nReceived the %s signal, stopping within %d steps\n\n",
1367 gmx_get_signal_name(), nsteps_stop);
1371 "\n\nReceived the %s signal, stopping within %d steps\n\n",
1372 gmx_get_signal_name(), nsteps_stop);
1374 handled_stop_condition = static_cast<int>(gmx_get_stop_condition());
1376 else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
1377 (mdrunOptions.maximumHoursToRun > 0 &&
1378 secondsSinceStart > mdrunOptions.maximumHoursToRun*60.0*60.0*0.99) &&
1379 signals[eglsSTOPCOND].sig == 0 && signals[eglsSTOPCOND].set == 0)
1381 /* Signal to terminate the run */
1382 signals[eglsSTOPCOND].sig = 1;
1385 fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",
1386 gmx_step_str(step, sbuf), mdrunOptions.maximumHoursToRun*0.99);
1388 fprintf(stderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n",
1389 gmx_step_str(step, sbuf), mdrunOptions.maximumHoursToRun*0.99);
1392 if (bResetCountersHalfMaxH && MASTER(cr) &&
1393 secondsSinceStart > mdrunOptions.maximumHoursToRun*60.0*60.0*0.495)
1395 /* Set flag that will communicate the signal to all ranks in the simulation */
1396 signals[eglsRESETCOUNTERS].sig = 1;
1399 /* In parallel we only have to check for checkpointing in steps
1400 * where we do global communication,
1401 * otherwise the other nodes don't know.
1403 const real cpt_period = mdrunOptions.checkpointOptions.period;
1404 if (MASTER(cr) && ((bGStat || !PAR(cr)) &&
1407 secondsSinceStart >= nchkpt*cpt_period*60.0)) &&
1408 signals[eglsCHKPT].set == 0)
1410 signals[eglsCHKPT].sig = 1;
1413 /* ######### START SECOND UPDATE STEP ################# */
1415 /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1418 if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1420 gmx_bool bIfRandomize;
1421 bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr);
1422 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1423 if (constr && bIfRandomize)
1425 constrain_velocities(step, nullptr,
1429 bCalcVir, do_log, do_ene);
1432 /* Box is changed in update() when we do pressure coupling,
1433 * but we should still use the old box for energy corrections and when
1434 * writing it to the energy file, so it matches the trajectory files for
1435 * the same timestep above. Make a copy in a separate array.
1437 copy_mat(state->box, lastbox);
1441 if (!bRerunMD || rerun_fr.bV || bForceUpdate)
1443 wallcycle_start(wcycle, ewcUPDATE);
1444 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1447 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1448 /* We can only do Berendsen coupling after we have summed
1449 * the kinetic energy or virial. Since the happens
1450 * in global_state after update, we should only do it at
1451 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1456 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1457 update_pcouple_before_coordinates(fplog, step, ir, state,
1458 parrinellorahmanMu, M,
1464 /* velocity half-step update */
1465 update_coords(step, ir, mdatoms, state, f, fcd,
1466 ekind, M, upd, etrtVELOCITY2,
1470 /* Above, initialize just copies ekinh into ekin,
1471 * it doesn't copy position (for VV),
1472 * and entire integrator for MD.
1475 if (ir->eI == eiVVAK)
1477 /* We probably only need md->homenr, not state->natoms */
1478 if (state->natoms > cbuf_nalloc)
1480 cbuf_nalloc = state->natoms;
1481 srenew(cbuf, cbuf_nalloc);
1483 copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1486 update_coords(step, ir, mdatoms, state, f, fcd,
1487 ekind, M, upd, etrtPOSITION, cr, constr);
1488 wallcycle_stop(wcycle, ewcUPDATE);
1490 constrain_coordinates(step, &dvdl_constr, state,
1492 wcycle, upd, constr,
1493 bCalcVir, do_log, do_ene);
1494 update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
1495 cr, nrnb, wcycle, upd, constr, do_log, do_ene);
1496 finish_update(ir, mdatoms,
1498 nrnb, wcycle, upd, constr);
1500 if (ir->eI == eiVVAK)
1502 /* erase F_EKIN and F_TEMP here? */
1503 /* just compute the kinetic energy at the half step to perform a trotter step */
1504 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1505 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1506 constr, &nullSignaller, lastbox,
1507 nullptr, &bSumEkinhOld,
1508 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1510 wallcycle_start(wcycle, ewcUPDATE);
1511 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1512 /* now we know the scaling, we can compute the positions again again */
1513 copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1515 update_coords(step, ir, mdatoms, state, f, fcd,
1516 ekind, M, upd, etrtPOSITION, cr, constr);
1517 wallcycle_stop(wcycle, ewcUPDATE);
1519 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1520 /* are the small terms in the shake_vir here due
1521 * to numerical errors, or are they important
1522 * physically? I'm thinking they are just errors, but not completely sure.
1523 * For now, will call without actually constraining, constr=NULL*/
1524 finish_update(ir, mdatoms,
1526 nrnb, wcycle, upd, nullptr);
1530 /* this factor or 2 correction is necessary
1531 because half of the constraint force is removed
1532 in the vv step, so we have to double it. See
1533 the Redmine issue #1255. It is not yet clear
1534 if the factor of 2 is exact, or just a very
1535 good approximation, and this will be
1536 investigated. The next step is to see if this
1537 can be done adding a dhdl contribution from the
1538 rattle step, but this is somewhat more
1539 complicated with the current code. Will be
1540 investigated, hopefully for 4.6.3. However,
1541 this current solution is much better than
1542 having it completely wrong.
1544 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1548 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1553 /* Need to unshift here */
1554 unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1557 if (vsite != nullptr)
1559 wallcycle_start(wcycle, ewcVSITECONSTR);
1560 if (graph != nullptr)
1562 shift_self(graph, state->box, as_rvec_array(state->x.data()));
1564 construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
1565 top->idef.iparams, top->idef.il,
1566 fr->ePBC, fr->bMolPBC, cr, state->box);
1568 if (graph != nullptr)
1570 unshift_self(graph, state->box, as_rvec_array(state->x.data()));
1572 wallcycle_stop(wcycle, ewcVSITECONSTR);
1575 /* ############## IF NOT VV, Calculate globals HERE ############ */
1576 /* With Leap-Frog we can skip compute_globals at
1577 * non-communication steps, but we need to calculate
1578 * the kinetic energy one step before communication.
1581 // Organize to do inter-simulation signalling on steps if
1582 // and when algorithms require it.
1583 bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1585 if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1587 // Since we're already communicating at this step, we
1588 // can propagate intra-simulation signals. Note that
1589 // check_nstglobalcomm has the responsibility for
1590 // choosing the value of nstglobalcomm that is one way
1591 // bGStat becomes true, so we can't get into a
1592 // situation where e.g. checkpointing can't be
1594 bool doIntraSimSignal = true;
1595 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1597 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
1598 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1601 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1602 (bGStat ? CGLO_GSTAT : 0)
1603 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
1604 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1605 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1606 | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
1608 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1610 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1611 top_global, top, state,
1612 &shouldCheckNumberOfBondedInteractions);
1616 /* ############# END CALC EKIN AND PRESSURE ################# */
1618 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1619 the virial that should probably be addressed eventually. state->veta has better properies,
1620 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1621 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1623 if (ir->efep != efepNO && (!EI_VV(ir->eI) || bRerunMD))
1625 /* Sum up the foreign energy and dhdl terms for md and sd.
1626 Currently done every step so that dhdl is correct in the .edr */
1627 sum_dhdl(enerd, state->lambda, ir->fepvals);
1630 update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1631 pres, force_vir, shake_vir,
1635 /* ################# END UPDATE STEP 2 ################# */
1636 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1638 /* The coordinates (x) were unshifted in update */
1641 /* We will not sum ekinh_old,
1642 * so signal that we still have to do it.
1644 bSumEkinhOld = TRUE;
1649 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1651 /* use the directly determined last velocity, not actually the averaged half steps */
1652 if (bTrotter && ir->eI == eiVV)
1654 enerd->term[F_EKIN] = last_ekin;
1656 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1658 if (integratorHasConservedEnergyQuantity(ir))
1662 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1666 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1669 /* ######### END PREPARING EDR OUTPUT ########### */
1675 if (fplog && do_log && bDoExpanded)
1677 /* only needed if doing expanded ensemble */
1678 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1679 state_global->dfhist, state->fep_state, ir->nstlog, step);
1683 upd_mdebin(mdebin, bDoDHDL, bCalcEnerStep,
1684 t, mdatoms->tmass, enerd, state,
1685 ir->fepvals, ir->expandedvals, lastbox,
1686 shake_vir, force_vir, total_vir, pres,
1687 ekind, mu_tot, constr);
1691 upd_mdebin_step(mdebin);
1694 gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
1695 gmx_bool do_or = do_per_step(step, ir->nstorireout);
1697 print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
1699 eprNORMAL, mdebin, fcd, groups, &(ir->opts), awh.get());
1703 pull_print_output(ir->pull_work, step, t);
1706 if (do_per_step(step, ir->nstlog))
1708 if (fflush(fplog) != 0)
1710 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1716 /* Have to do this part _after_ outputting the logfile and the edr file */
1717 /* Gets written into the state at the beginning of next loop*/
1718 state->fep_state = lamnew;
1720 /* Print the remaining wall clock time for the run */
1721 if (isMasterSimMasterRank(ms, cr) &&
1722 (do_verbose || gmx_got_usr_signal()) &&
1727 fprintf(stderr, "\n");
1729 print_time(stderr, walltime_accounting, step, ir, cr);
1732 /* Ion/water position swapping.
1733 * Not done in last step since trajectory writing happens before this call
1734 * in the MD loop and exchanges would be lost anyway. */
1735 bNeedRepartition = FALSE;
1736 if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1737 do_per_step(step, ir->swap->nstswap))
1739 bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
1740 bRerunMD ? rerun_fr.x : as_rvec_array(state->x.data()),
1741 bRerunMD ? rerun_fr.box : state->box,
1742 MASTER(cr) && mdrunOptions.verbose,
1745 if (bNeedRepartition && DOMAINDECOMP(cr))
1747 dd_collect_state(cr->dd, state, state_global);
1751 /* Replica exchange */
1755 bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
1756 state_global, enerd,
1760 if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1762 dd_partition_system(fplog, mdlog, step, cr, TRUE, 1,
1763 state_global, top_global, ir,
1764 state, &f, mdAtoms, top, fr,
1766 nrnb, wcycle, FALSE);
1767 shouldCheckNumberOfBondedInteractions = true;
1768 update_realloc(upd, state->natoms);
1773 startingFromCheckpoint = false;
1775 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1776 /* With all integrators, except VV, we need to retain the pressure
1777 * at the current step for coupling at the next step.
1779 if ((state->flags & (1<<estPRES_PREV)) &&
1781 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1783 /* Store the pressure in t_state for pressure coupling
1784 * at the next MD step.
1786 copy_mat(pres, state->pres_prev);
1789 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1791 if ( (membed != nullptr) && (!bLastStep) )
1793 rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1800 /* read next frame from input trajectory */
1801 bLastStep = !read_next_frame(oenv, status, &rerun_fr);
1806 rerun_parallel_comm(cr, &rerun_fr, &bLastStep);
1810 cycles = wallcycle_stop(wcycle, ewcSTEP);
1811 if (DOMAINDECOMP(cr) && wcycle)
1813 dd_cycles_add(cr->dd, cycles, ddCyclStep);
1816 if (!bRerunMD || !rerun_fr.bStep)
1818 /* increase the MD step number */
1823 /* TODO make a counter-reset module */
1824 /* If it is time to reset counters, set a flag that remains
1825 true until counters actually get reset */
1826 if (step_rel == wcycle_get_reset_counters(wcycle) ||
1827 signals[eglsRESETCOUNTERS].set != 0)
1829 if (pme_loadbal_is_active(pme_loadbal))
1831 /* Do not permit counter reset while PME load
1832 * balancing is active. The only purpose for resetting
1833 * counters is to measure reliable performance data,
1834 * and that can't be done before balancing
1837 * TODO consider fixing this by delaying the reset
1838 * until after load balancing completes,
1839 * e.g. https://gerrit.gromacs.org/#/c/4964/2 */
1840 gmx_fatal(FARGS, "PME tuning was still active when attempting to "
1841 "reset mdrun counters at step %" PRId64 ". Try "
1842 "resetting counters later in the run, e.g. with gmx "
1843 "mdrun -resetstep.", step);
1845 reset_all_counters(fplog, mdlog, cr, step, wcycle, nrnb, walltime_accounting,
1846 use_GPU(fr->nbv) ? fr->nbv : nullptr, fr->pmedata);
1847 wcycle_set_reset_counters(wcycle, -1);
1848 if (!thisRankHasDuty(cr, DUTY_PME))
1850 /* Tell our PME node to reset its counters */
1851 gmx_pme_send_resetcounters(cr, step);
1853 /* If mdrun -maxh -resethway was active, it can only trigger once */
1854 bResetCountersHalfMaxH = FALSE; /* TODO move this to where signals[eglsRESETCOUNTERS].sig is set */
1855 /* Reset can only happen once, so clear the triggering flag. */
1856 signals[eglsRESETCOUNTERS].set = 0;
1859 /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1860 IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
1863 /* End of main MD loop */
1865 /* Closing TNG files can include compressing data. Therefore it is good to do that
1866 * before stopping the time measurements. */
1867 mdoutf_tng_close(outf);
1869 /* Stop measuring walltime */
1870 walltime_accounting_end_time(walltime_accounting);
1872 if (bRerunMD && MASTER(cr))
1877 if (!thisRankHasDuty(cr, DUTY_PME))
1879 /* Tell the PME only node to finish */
1880 gmx_pme_send_finish(cr);
1885 if (ir->nstcalcenergy > 0 && !bRerunMD)
1887 print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
1888 eprAVER, mdebin, fcd, groups, &(ir->opts), awh.get());
1891 done_mdebin(mdebin);
1896 pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
1899 done_shellfc(fplog, shellfc, step_rel);
1901 if (useReplicaExchange && MASTER(cr))
1903 print_replica_exchange_statistics(fplog, repl_ex);
1906 // Clean up swapcoords
1907 if (ir->eSwapCoords != eswapNO)
1909 finish_swapcoords(ir->swap);
1912 /* IMD cleanup, if bIMD is TRUE. */
1913 IMD_finalize(ir->bIMD, ir->imd);
1915 walltime_accounting_set_nsteps_done(walltime_accounting, step);
1916 if (step_rel >= wcycle_get_reset_counters(wcycle) &&
1917 signals[eglsRESETCOUNTERS].set == 0 &&
1918 !bResetCountersHalfMaxH)
1920 walltime_accounting_set_valid_finish(walltime_accounting);
1923 destroy_enerdata(enerd);