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,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
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
54 #include "gromacs/awh/awh.h"
55 #include "gromacs/commandline/filenm.h"
56 #include "gromacs/domdec/collect.h"
57 #include "gromacs/domdec/dlbtiming.h"
58 #include "gromacs/domdec/domdec.h"
59 #include "gromacs/domdec/domdec_network.h"
60 #include "gromacs/domdec/domdec_struct.h"
61 #include "gromacs/domdec/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/energyoutput.h"
80 #include "gromacs/mdlib/expanded.h"
81 #include "gromacs/mdlib/force.h"
82 #include "gromacs/mdlib/force_flags.h"
83 #include "gromacs/mdlib/forcerec.h"
84 #include "gromacs/mdlib/md_support.h"
85 #include "gromacs/mdlib/mdatoms.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/ns.h"
91 #include "gromacs/mdlib/resethandler.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/stophandler.h"
97 #include "gromacs/mdlib/tgroup.h"
98 #include "gromacs/mdlib/trajectory_writing.h"
99 #include "gromacs/mdlib/update.h"
100 #include "gromacs/mdlib/vcm.h"
101 #include "gromacs/mdlib/vsite.h"
102 #include "gromacs/mdtypes/awh_history.h"
103 #include "gromacs/mdtypes/awh_params.h"
104 #include "gromacs/mdtypes/commrec.h"
105 #include "gromacs/mdtypes/df_history.h"
106 #include "gromacs/mdtypes/energyhistory.h"
107 #include "gromacs/mdtypes/fcdata.h"
108 #include "gromacs/mdtypes/forcerec.h"
109 #include "gromacs/mdtypes/group.h"
110 #include "gromacs/mdtypes/inputrec.h"
111 #include "gromacs/mdtypes/interaction_const.h"
112 #include "gromacs/mdtypes/md_enums.h"
113 #include "gromacs/mdtypes/mdatom.h"
114 #include "gromacs/mdtypes/observableshistory.h"
115 #include "gromacs/mdtypes/pullhistory.h"
116 #include "gromacs/mdtypes/state.h"
117 #include "gromacs/nbnxm/nbnxm.h"
118 #include "gromacs/pbcutil/mshift.h"
119 #include "gromacs/pbcutil/pbc.h"
120 #include "gromacs/pulling/output.h"
121 #include "gromacs/pulling/pull.h"
122 #include "gromacs/swap/swapcoords.h"
123 #include "gromacs/timing/wallcycle.h"
124 #include "gromacs/timing/walltime_accounting.h"
125 #include "gromacs/topology/atoms.h"
126 #include "gromacs/topology/idef.h"
127 #include "gromacs/topology/mtop_util.h"
128 #include "gromacs/topology/topology.h"
129 #include "gromacs/trajectory/trajectoryframe.h"
130 #include "gromacs/utility/basedefinitions.h"
131 #include "gromacs/utility/cstringutil.h"
132 #include "gromacs/utility/fatalerror.h"
133 #include "gromacs/utility/logger.h"
134 #include "gromacs/utility/real.h"
135 #include "gromacs/utility/smalloc.h"
137 #include "integrator.h"
138 #include "replicaexchange.h"
141 #include "corewrap.h"
144 using gmx::SimulationSignaller;
146 void gmx::Integrator::do_md()
148 // TODO Historically, the EM and MD "integrators" used different
149 // names for the t_inputrec *parameter, but these must have the
150 // same name, now that it's a member of a struct. We use this ir
151 // alias to avoid a large ripple of nearly useless changes.
152 // t_inputrec is being replaced by IMdpOptionsProvider, so this
153 // will go away eventually.
154 t_inputrec *ir = inputrec;
155 int64_t step, step_rel;
156 double t = ir->init_t, t0 = ir->init_t, lam0[efptNR];
157 gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
158 gmx_bool bNS, bNStList, bStopCM,
159 bFirstStep, bInitStep, bLastStep = FALSE;
160 gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
161 gmx_bool do_ene, do_log, do_verbose;
162 gmx_bool bMasterState;
163 int force_flags, cglo_flags;
164 tensor force_vir = {{0}}, shake_vir = {{0}}, total_vir = {{0}},
165 tmp_vir = {{0}}, pres = {{0}};
168 matrix parrinellorahmanMu, M;
169 gmx_repl_ex_t repl_ex = nullptr;
171 gmx_enerdata_t *enerd;
172 PaddedVector<gmx::RVec> f {};
173 gmx_global_stat_t gstat;
174 t_graph *graph = nullptr;
175 gmx_groups_t *groups;
176 gmx_shellfc_t *shellfc;
177 gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
178 gmx_bool bTemp, bPres, bTrotter;
180 rvec *cbuf = nullptr;
187 real saved_conserved_quantity = 0;
190 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
192 /* PME load balancing data for GPU kernels */
193 gmx_bool bPMETune = FALSE;
194 gmx_bool bPMETunePrinting = FALSE;
197 gmx_bool bIMDstep = FALSE;
199 /* Domain decomposition could incorrectly miss a bonded
200 interaction, but checking for that requires a global
201 communication stage, which does not otherwise happen in DD
202 code. So we do that alongside the first global energy reduction
203 after a new DD is made. These variables handle whether the
204 check happens, and the result it returns. */
205 bool shouldCheckNumberOfBondedInteractions = false;
206 int totalNumberOfBondedInteractions = -1;
208 SimulationSignals signals;
209 // Most global communnication stages don't propagate mdrun
210 // signals, and will use this object to achieve that.
211 SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
213 if (!mdrunOptions.writeConfout)
215 // This is on by default, and the main known use case for
216 // turning it off is for convenience in benchmarking, which is
217 // something that should not show up in the general user
219 GMX_LOG(mdlog.info).asParagraph().
220 appendText("The -noconfout functionality is deprecated, and may be removed in a future version.");
223 /* md-vv uses averaged full step velocities for T-control
224 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
225 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
226 bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
228 const bool bRerunMD = false;
229 int nstglobalcomm = mdrunOptions.globalCommunicationInterval;
231 nstglobalcomm = check_nstglobalcomm(mdlog, nstglobalcomm, ir, cr);
232 bGStatEveryStep = (nstglobalcomm == 1);
234 groups = &top_global->groups;
236 std::unique_ptr<EssentialDynamics> ed = nullptr;
237 if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
239 /* Initialize essential dynamics sampling */
240 ed = init_edsam(mdlog,
241 opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
244 state_global, observablesHistory,
245 oenv, mdrunOptions.continuationOptions.appendFiles);
248 initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda, lam0);
249 Update upd(ir, deform);
250 bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
251 if (!mdrunOptions.continuationOptions.appendFiles)
253 pleaseCiteCouplingAlgorithms(fplog, *ir);
256 gmx_mdoutf *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, ir, top_global, oenv, wcycle);
257 gmx::EnergyOutput energyOutput;
258 energyOutput.prepare(mdoutf_get_fp_ene(outf), top_global, ir, mdoutf_get_fp_dhdl(outf));
260 /* Energy terms and groups */
262 init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
265 /* Kinetic energy data */
266 std::unique_ptr<gmx_ekindata_t> eKinData = std::make_unique<gmx_ekindata_t>();
267 gmx_ekindata_t *ekind = eKinData.get();
268 init_ekindata(fplog, top_global, &(ir->opts), ekind);
269 /* Copy the cos acceleration to the groups struct */
270 ekind->cosacc.cos_accel = ir->cos_accel;
272 gstat = global_stat_init(ir);
274 /* Check for polarizable models and flexible constraints */
275 shellfc = init_shell_flexcon(fplog,
276 top_global, constr ? constr->numFlexibleConstraints() : 0,
277 ir->nstcalcenergy, DOMAINDECOMP(cr));
280 double io = compute_io(ir, top_global->natoms, groups, energyOutput.numEnergyTerms(), 1);
281 if ((io > 2000) && MASTER(cr))
284 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
289 /* Set up interactive MD (IMD) */
290 init_IMD(ir, cr, ms, top_global, fplog, ir->nstcalcenergy,
291 MASTER(cr) ? state_global->x.rvec_array() : nullptr,
292 nfile, fnm, oenv, mdrunOptions);
294 // Local state only becomes valid now.
295 std::unique_ptr<t_state> stateInstance;
298 if (DOMAINDECOMP(cr))
300 dd_init_local_top(*top_global, &top);
302 stateInstance = std::make_unique<t_state>();
303 state = stateInstance.get();
304 dd_init_local_state(cr->dd, state_global, state);
306 /* Distribute the charge groups over the nodes from the master node */
307 dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
308 state_global, *top_global, ir,
309 state, &f, mdAtoms, &top, fr,
311 nrnb, nullptr, FALSE);
312 shouldCheckNumberOfBondedInteractions = true;
313 upd.setNumAtoms(state->natoms);
317 state_change_natoms(state_global, state_global->natoms);
318 f.resizeWithPadding(state_global->natoms);
319 /* Copy the pointer to the global state */
320 state = state_global;
322 /* Generate and initialize new topology */
323 mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr,
324 &graph, mdAtoms, constr, vsite, shellfc);
326 upd.setNumAtoms(state->natoms);
329 auto mdatoms = mdAtoms->mdatoms();
331 // NOTE: The global state is no longer used at this point.
332 // But state_global is still used as temporary storage space for writing
333 // the global state to file and potentially for replica exchange.
334 // (Global topology should persist.)
336 update_mdatoms(mdatoms, state->lambda[efptMASS]);
338 const ContinuationOptions &continuationOptions = mdrunOptions.continuationOptions;
339 bool startingFromCheckpoint = continuationOptions.startedFromCheckpoint;
343 /* Check nstexpanded here, because the grompp check was broken */
344 if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
346 gmx_fatal(FARGS, "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
348 init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
353 if (startingFromCheckpoint)
355 /* Restore from energy history if appending to output files */
356 if (continuationOptions.appendFiles)
358 /* If no history is available (because a checkpoint is from before
359 * it was written) make a new one later, otherwise restore it.
361 if (observablesHistory->energyHistory)
363 energyOutput.restoreFromEnergyHistory(*observablesHistory->energyHistory);
366 else if (observablesHistory->energyHistory)
368 /* We might have read an energy history from checkpoint.
369 * As we are not appending, we want to restart the statistics.
370 * Free the allocated memory and reset the counts.
372 observablesHistory->energyHistory = {};
373 /* We might have read a pull history from checkpoint.
374 * We will still want to keep the statistics, so that the files
375 * can be joined and still be meaningful.
376 * This means that observablesHistory->pullHistory
377 * should not be reset.
381 if (!observablesHistory->energyHistory)
383 observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
385 if (!observablesHistory->pullHistory)
387 observablesHistory->pullHistory = std::make_unique<PullHistory>();
389 /* Set the initial energy history */
390 energyOutput.fillEnergyHistory(observablesHistory->energyHistory.get());
393 preparePrevStepPullCom(ir, mdatoms, state, state_global, cr, startingFromCheckpoint);
395 // TODO: Remove this by converting AWH into a ForceProvider
396 auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingFromCheckpoint,
398 opt2fn("-awh", nfile, fnm), ir->pull_work);
400 const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
401 if (useReplicaExchange && MASTER(cr))
403 repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
406 /* PME tuning is only supported in the Verlet scheme, with PME for
407 * Coulomb. It is not supported with only LJ PME. */
408 bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) &&
409 !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
411 pme_load_balancing_t *pme_loadbal = nullptr;
414 pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box,
415 *fr->ic, fr->nbv->pairlistSets().params(), fr->pmedata, use_GPU(fr->nbv),
419 if (!ir->bContinuation)
421 if (state->flags & (1 << estV))
423 auto v = makeArrayRef(state->v);
424 /* Set the velocities of vsites, shells and frozen atoms to zero */
425 for (i = 0; i < mdatoms->homenr; i++)
427 if (mdatoms->ptype[i] == eptVSite ||
428 mdatoms->ptype[i] == eptShell)
432 else if (mdatoms->cFREEZE)
434 for (m = 0; m < DIM; m++)
436 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
447 /* Constrain the initial coordinates and velocities */
448 do_constrain_first(fplog, constr, ir, mdatoms, state);
452 /* Construct the virtual sites for the initial configuration */
453 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
454 top.idef.iparams, top.idef.il,
455 fr->ePBC, fr->bMolPBC, cr, state->box);
459 if (ir->efep != efepNO)
461 /* Set free energy calculation frequency as the greatest common
462 * denominator of nstdhdl and repl_ex_nst. */
463 nstfep = ir->fepvals->nstdhdl;
466 nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
468 if (useReplicaExchange)
470 nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
474 /* Be REALLY careful about what flags you set here. You CANNOT assume
475 * this is the first step, since we might be restarting from a checkpoint,
476 * and in that case we should not do any modifications to the state.
478 bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
480 if (continuationOptions.haveReadEkin)
482 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
485 cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT
486 | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
487 | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
488 | (continuationOptions.haveReadEkin ? CGLO_READEKIN : 0));
490 bSumEkinhOld = FALSE;
492 t_vcm vcm(top_global->groups, *ir);
493 reportComRemovalInfo(fplog, vcm);
495 /* To minimize communication, compute_globals computes the COM velocity
496 * and the kinetic energy for the velocities without COM motion removed.
497 * Thus to get the kinetic energy without the COM contribution, we need
498 * to call compute_globals twice.
500 for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
502 int cglo_flags_iteration = cglo_flags;
503 if (bStopCM && cgloIteration == 0)
505 cglo_flags_iteration |= CGLO_STOPCM;
506 cglo_flags_iteration &= ~CGLO_TEMPERATURE;
508 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
509 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
510 constr, &nullSignaller, state->box,
511 &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
512 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
514 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
515 top_global, &top, state,
516 &shouldCheckNumberOfBondedInteractions);
517 if (ir->eI == eiVVAK)
519 /* a second call to get the half step temperature initialized as well */
520 /* we do the same call as above, but turn the pressure off -- internally to
521 compute_globals, this is recognized as a velocity verlet half-step
522 kinetic energy calculation. This minimized excess variables, but
523 perhaps loses some logic?*/
525 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
526 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
527 constr, &nullSignaller, state->box,
528 nullptr, &bSumEkinhOld,
529 cglo_flags & ~CGLO_PRESSURE);
532 /* Calculate the initial half step temperature, and save the ekinh_old */
533 if (!continuationOptions.startedFromCheckpoint)
535 for (i = 0; (i < ir->opts.ngtc); i++)
537 copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
541 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
542 temperature control */
543 auto trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
547 if (!ir->bContinuation)
549 if (constr && ir->eConstrAlg == econtLINCS)
552 "RMS relative constraint deviation after constraining: %.2e\n",
555 if (EI_STATE_VELOCITY(ir->eI))
557 real temp = enerd->term[F_TEMP];
560 /* Result of Ekin averaged over velocities of -half
561 * and +half step, while we only have -half step here.
565 fprintf(fplog, "Initial temperature: %g K\n", temp);
570 fprintf(stderr, "starting mdrun '%s'\n",
571 *(top_global->name));
574 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
578 sprintf(tbuf, "%s", "infinite");
580 if (ir->init_step > 0)
582 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
583 gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
584 gmx_step_str(ir->init_step, sbuf2),
585 ir->init_step*ir->delta_t);
589 fprintf(stderr, "%s steps, %s ps.\n",
590 gmx_step_str(ir->nsteps, sbuf), tbuf);
592 fprintf(fplog, "\n");
595 walltime_accounting_start_time(walltime_accounting);
596 wallcycle_start(wcycle, ewcRUN);
597 print_start(fplog, cr, walltime_accounting, "mdrun");
600 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
601 int chkpt_ret = fcCheckPointParallel( cr->nodeid,
605 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
609 /***********************************************************
613 ************************************************************/
616 /* Skip the first Nose-Hoover integration when we get the state from tpx */
617 bInitStep = !startingFromCheckpoint || EI_VV(ir->eI);
618 bSumEkinhOld = FALSE;
620 bNeedRepartition = FALSE;
622 bool simulationsShareState = false;
623 int nstSignalComm = nstglobalcomm;
625 // TODO This implementation of ensemble orientation restraints is nasty because
626 // a user can't just do multi-sim with single-sim orientation restraints.
627 bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
628 bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
630 // Replica exchange, ensemble restraints and AWH need all
631 // simulations to remain synchronized, so they need
632 // checkpoints and stop conditions to act on the same step, so
633 // the propagation of such signals must take place between
634 // simulations, not just within simulations.
635 // TODO: Make algorithm initializers set these flags.
636 simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
638 if (simulationsShareState)
640 // Inter-simulation signal communication does not need to happen
641 // often, so we use a minimum of 200 steps to reduce overhead.
642 const int c_minimumInterSimulationSignallingInterval = 200;
643 nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
647 auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
648 compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
649 MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
650 mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
652 auto checkpointHandler = std::make_unique<CheckpointHandler>(
653 compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]),
654 simulationsShareState, ir->nstlist == 0, MASTER(cr),
655 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
657 const bool resetCountersIsLocal = true;
658 auto resetHandler = std::make_unique<ResetHandler>(
659 compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]), !resetCountersIsLocal,
660 ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
661 mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
663 const DDBalanceRegionHandler ddBalanceRegionHandler(cr);
665 step = ir->init_step;
668 // TODO extract this to new multi-simulation module
669 if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
671 if (!multisim_int_all_are_equal(ms, ir->nsteps))
673 GMX_LOG(mdlog.warning).appendText(
674 "Note: The number of steps is not consistent across multi simulations,\n"
675 "but we are proceeding anyway!");
677 if (!multisim_int_all_are_equal(ms, ir->init_step))
679 GMX_LOG(mdlog.warning).appendText(
680 "Note: The initial step is not consistent across multi simulations,\n"
681 "but we are proceeding anyway!");
685 /* and stop now if we should */
686 bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
690 /* Determine if this is a neighbor search step */
691 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
693 if (bPMETune && bNStList)
695 /* PME grid + cut-off optimization with GPUs or PME nodes */
696 pme_loadbal_do(pme_loadbal, cr,
697 (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
705 wallcycle_start(wcycle, ewcSTEP);
707 bLastStep = (step_rel == ir->nsteps);
708 t = t0 + step*ir->delta_t;
710 // TODO Refactor this, so that nstfep does not need a default value of zero
711 if (ir->efep != efepNO || ir->bSimTemp)
713 /* find and set the current lambdas */
714 setCurrentLambdasLocal(step, ir->fepvals, lam0, state);
716 bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
717 bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
718 bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
719 && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
722 bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
723 do_per_step(step, replExParams.exchangeInterval));
725 if (doSimulatedAnnealing)
727 update_annealing_target_temp(ir, t, &upd);
730 /* Stop Center of Mass motion */
731 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
733 /* Determine whether or not to do Neighbour Searching */
734 bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
736 bLastStep = bLastStep || stopHandler->stoppingAfterCurrentStep(bNS);
738 /* do_log triggers energy and virial calculation. Because this leads
739 * to different code paths, forces can be different. Thus for exact
740 * continuation we should avoid extra log output.
741 * Note that the || bLastStep can result in non-exact continuation
742 * beyond the last step. But we don't consider that to be an issue.
744 do_log = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep;
745 do_verbose = mdrunOptions.verbose &&
746 (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
748 if (bNS && !(bFirstStep && ir->bContinuation))
750 bMasterState = FALSE;
751 /* Correct the new box if it is too skewed */
752 if (inputrecDynamicBox(ir))
754 if (correct_box(fplog, step, state->box, graph))
759 if (DOMAINDECOMP(cr) && bMasterState)
761 dd_collect_state(cr->dd, state, state_global);
764 if (DOMAINDECOMP(cr))
766 /* Repartition the domain decomposition */
767 dd_partition_system(fplog, mdlog, step, cr,
768 bMasterState, nstglobalcomm,
769 state_global, *top_global, ir,
770 state, &f, mdAtoms, &top, fr,
773 do_verbose && !bPMETunePrinting);
774 shouldCheckNumberOfBondedInteractions = true;
775 upd.setNumAtoms(state->natoms);
779 if (MASTER(cr) && do_log)
781 print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
784 if (ir->efep != efepNO)
786 update_mdatoms(mdatoms, state->lambda[efptMASS]);
792 /* We need the kinetic energy at minus the half step for determining
793 * the full step kinetic energy and possibly for T-coupling.*/
794 /* This may not be quite working correctly yet . . . . */
795 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
796 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
797 constr, &nullSignaller, state->box,
798 &totalNumberOfBondedInteractions, &bSumEkinhOld,
799 CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
800 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
801 top_global, &top, state,
802 &shouldCheckNumberOfBondedInteractions);
804 clear_mat(force_vir);
806 checkpointHandler->decideIfCheckpointingThisStep(bNS, bFirstStep, bLastStep);
808 /* Determine the energy and pressure:
809 * at nstcalcenergy steps and at energy output steps (set below).
811 if (EI_VV(ir->eI) && (!bInitStep))
813 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
814 bCalcVir = bCalcEnerStep ||
815 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
819 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
820 bCalcVir = bCalcEnerStep ||
821 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
823 bCalcEner = bCalcEnerStep;
825 do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
827 if (do_ene || do_log || bDoReplEx)
833 /* Do we need global communication ? */
834 bGStat = (bCalcVir || bCalcEner || bStopCM ||
835 do_per_step(step, nstglobalcomm) ||
836 (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
838 force_flags = (GMX_FORCE_STATECHANGED |
839 ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0) |
840 GMX_FORCE_ALLFORCES |
841 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
842 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
843 (bDoFEP ? GMX_FORCE_DHDL : 0)
848 /* Now is the time to relax the shells */
849 relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
850 enforcedRotation, step,
851 ir, bNS, force_flags, &top,
853 state, f.arrayRefWithPadding(), force_vir, mdatoms,
854 nrnb, wcycle, graph, groups,
855 shellfc, fr, ppForceWorkload, t, mu_tot,
857 ddBalanceRegionHandler);
861 /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
862 (or the AWH update will be performed twice for one step when continuing). It would be best to
863 call this update function from do_md_trajectory_writing but that would occur after do_force.
864 One would have to divide the update_awh function into one function applying the AWH force
865 and one doing the AWH bias update. The update AWH bias function could then be called after
866 do_md_trajectory_writing (then containing update_awh_history).
867 The checkpointing will in the future probably moved to the start of the md loop which will
868 rid of this issue. */
869 if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
871 awh->updateHistory(state_global->awhHistory.get());
874 /* The coordinates (x) are shifted (to get whole molecules)
876 * This is parallellized as well, and does communication too.
877 * Check comments in sim_util.c
879 do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation,
880 step, nrnb, wcycle, &top, groups,
881 state->box, state->x.arrayRefWithPadding(), &state->hist,
882 f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
883 state->lambda, graph,
884 fr, ppForceWorkload, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
885 (bNS ? GMX_FORCE_NS : 0) | force_flags,
886 ddBalanceRegionHandler);
889 if (EI_VV(ir->eI) && !startingFromCheckpoint)
890 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
892 rvec *vbuf = nullptr;
894 wallcycle_start(wcycle, ewcUPDATE);
895 if (ir->eI == eiVV && bInitStep)
897 /* if using velocity verlet with full time step Ekin,
898 * take the first half step only to compute the
899 * virial for the first step. From there,
900 * revert back to the initial coordinates
901 * so that the input is actually the initial step.
903 snew(vbuf, state->natoms);
904 copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
908 /* this is for NHC in the Ekin(t+dt/2) version of vv */
909 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
912 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
913 ekind, M, &upd, etrtVELOCITY1,
916 wallcycle_stop(wcycle, ewcUPDATE);
917 constrain_velocities(step, nullptr,
921 bCalcVir, do_log, do_ene);
922 wallcycle_start(wcycle, ewcUPDATE);
923 /* if VV, compute the pressure and constraints */
924 /* For VV2, we strictly only need this if using pressure
925 * control, but we really would like to have accurate pressures
927 * Think about ways around this in the future?
928 * For now, keep this choice in comments.
930 /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
931 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
933 bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
934 if (bCalcEner && ir->eI == eiVVAK)
938 /* for vv, the first half of the integration actually corresponds to the previous step.
939 So we need information from the last step in the first half of the integration */
940 if (bGStat || do_per_step(step-1, nstglobalcomm))
942 wallcycle_stop(wcycle, ewcUPDATE);
943 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
944 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
945 constr, &nullSignaller, state->box,
946 &totalNumberOfBondedInteractions, &bSumEkinhOld,
947 (bGStat ? CGLO_GSTAT : 0)
948 | (bCalcEner ? CGLO_ENERGY : 0)
949 | (bTemp ? CGLO_TEMPERATURE : 0)
950 | (bPres ? CGLO_PRESSURE : 0)
951 | (bPres ? CGLO_CONSTRAINT : 0)
952 | (bStopCM ? CGLO_STOPCM : 0)
953 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
956 /* explanation of above:
957 a) We compute Ekin at the full time step
958 if 1) we are using the AveVel Ekin, and it's not the
959 initial step, or 2) if we are using AveEkin, but need the full
960 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
961 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
962 EkinAveVel because it's needed for the pressure */
963 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
964 top_global, &top, state,
965 &shouldCheckNumberOfBondedInteractions);
966 wallcycle_start(wcycle, ewcUPDATE);
968 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
973 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
974 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
976 /* TODO This is only needed when we're about to write
977 * a checkpoint, because we use it after the restart
978 * (in a kludge?). But what should we be doing if
979 * startingFromCheckpoint or bInitStep are true? */
980 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
982 copy_mat(shake_vir, state->svir_prev);
983 copy_mat(force_vir, state->fvir_prev);
985 if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
987 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
988 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
989 enerd->term[F_EKIN] = trace(ekind->ekin);
994 wallcycle_stop(wcycle, ewcUPDATE);
995 /* We need the kinetic energy at minus the half step for determining
996 * the full step kinetic energy and possibly for T-coupling.*/
997 /* This may not be quite working correctly yet . . . . */
998 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
999 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1000 constr, &nullSignaller, state->box,
1001 nullptr, &bSumEkinhOld,
1002 CGLO_GSTAT | CGLO_TEMPERATURE);
1003 wallcycle_start(wcycle, ewcUPDATE);
1006 /* if it's the initial step, we performed this first step just to get the constraint virial */
1007 if (ir->eI == eiVV && bInitStep)
1009 copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
1012 wallcycle_stop(wcycle, ewcUPDATE);
1015 /* compute the conserved quantity */
1018 saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1021 last_ekin = enerd->term[F_EKIN];
1023 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1025 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1027 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1028 if (ir->efep != efepNO)
1030 sum_dhdl(enerd, state->lambda, ir->fepvals);
1034 /* ######## END FIRST UPDATE STEP ############## */
1035 /* ######## If doing VV, we now have v(dt) ###### */
1038 /* perform extended ensemble sampling in lambda - we don't
1039 actually move to the new state before outputting
1040 statistics, but if performing simulated tempering, we
1041 do update the velocities and the tau_t. */
1043 lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v.rvec_array(), mdatoms);
1044 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1047 copy_df_history(state_global->dfhist, state->dfhist);
1051 /* Now we have the energies and forces corresponding to the
1052 * coordinates at time t. We must output all of this before
1055 do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1056 ir, state, state_global, observablesHistory,
1058 outf, energyOutput, ekind, f,
1059 checkpointHandler->isCheckpointingStep(),
1060 bRerunMD, bLastStep,
1061 mdrunOptions.writeConfout,
1063 /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1064 bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x.rvec_array(), ir, t, wcycle);
1066 /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1067 if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1069 copy_mat(state->svir_prev, shake_vir);
1070 copy_mat(state->fvir_prev, force_vir);
1073 stopHandler->setSignal();
1074 resetHandler->setSignal(walltime_accounting);
1076 if (bGStat || !PAR(cr))
1078 /* In parallel we only have to check for checkpointing in steps
1079 * where we do global communication,
1080 * otherwise the other nodes don't know.
1082 checkpointHandler->setSignal(walltime_accounting);
1085 /* ######### START SECOND UPDATE STEP ################# */
1087 /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1090 if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1092 gmx_bool bIfRandomize;
1093 bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
1094 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1095 if (constr && bIfRandomize)
1097 constrain_velocities(step, nullptr,
1101 bCalcVir, do_log, do_ene);
1104 /* Box is changed in update() when we do pressure coupling,
1105 * but we should still use the old box for energy corrections and when
1106 * writing it to the energy file, so it matches the trajectory files for
1107 * the same timestep above. Make a copy in a separate array.
1109 copy_mat(state->box, lastbox);
1113 wallcycle_start(wcycle, ewcUPDATE);
1114 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1117 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1118 /* We can only do Berendsen coupling after we have summed
1119 * the kinetic energy or virial. Since the happens
1120 * in global_state after update, we should only do it at
1121 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1126 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1127 update_pcouple_before_coordinates(fplog, step, ir, state,
1128 parrinellorahmanMu, M,
1134 /* velocity half-step update */
1135 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1136 ekind, M, &upd, etrtVELOCITY2,
1140 /* Above, initialize just copies ekinh into ekin,
1141 * it doesn't copy position (for VV),
1142 * and entire integrator for MD.
1145 if (ir->eI == eiVVAK)
1147 /* We probably only need md->homenr, not state->natoms */
1148 if (state->natoms > cbuf_nalloc)
1150 cbuf_nalloc = state->natoms;
1151 srenew(cbuf, cbuf_nalloc);
1153 copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1156 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1157 ekind, M, &upd, etrtPOSITION, cr, constr);
1158 wallcycle_stop(wcycle, ewcUPDATE);
1160 constrain_coordinates(step, &dvdl_constr, state,
1163 bCalcVir, do_log, do_ene);
1164 update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
1165 cr, nrnb, wcycle, &upd, constr, do_log, do_ene);
1166 finish_update(ir, mdatoms,
1168 nrnb, wcycle, &upd, constr);
1170 if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1172 updatePrevStepPullCom(ir->pull_work, state);
1175 if (ir->eI == eiVVAK)
1177 /* erase F_EKIN and F_TEMP here? */
1178 /* just compute the kinetic energy at the half step to perform a trotter step */
1179 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1180 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1181 constr, &nullSignaller, lastbox,
1182 nullptr, &bSumEkinhOld,
1183 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1185 wallcycle_start(wcycle, ewcUPDATE);
1186 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1187 /* now we know the scaling, we can compute the positions again again */
1188 copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1190 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1191 ekind, M, &upd, etrtPOSITION, cr, constr);
1192 wallcycle_stop(wcycle, ewcUPDATE);
1194 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1195 /* are the small terms in the shake_vir here due
1196 * to numerical errors, or are they important
1197 * physically? I'm thinking they are just errors, but not completely sure.
1198 * For now, will call without actually constraining, constr=NULL*/
1199 finish_update(ir, mdatoms,
1201 nrnb, wcycle, &upd, nullptr);
1205 /* this factor or 2 correction is necessary
1206 because half of the constraint force is removed
1207 in the vv step, so we have to double it. See
1208 the Redmine issue #1255. It is not yet clear
1209 if the factor of 2 is exact, or just a very
1210 good approximation, and this will be
1211 investigated. The next step is to see if this
1212 can be done adding a dhdl contribution from the
1213 rattle step, but this is somewhat more
1214 complicated with the current code. Will be
1215 investigated, hopefully for 4.6.3. However,
1216 this current solution is much better than
1217 having it completely wrong.
1219 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1223 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1226 if (vsite != nullptr)
1228 wallcycle_start(wcycle, ewcVSITECONSTR);
1229 if (graph != nullptr)
1231 shift_self(graph, state->box, state->x.rvec_array());
1233 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
1234 top.idef.iparams, top.idef.il,
1235 fr->ePBC, fr->bMolPBC, cr, state->box);
1237 if (graph != nullptr)
1239 unshift_self(graph, state->box, state->x.rvec_array());
1241 wallcycle_stop(wcycle, ewcVSITECONSTR);
1244 /* ############## IF NOT VV, Calculate globals HERE ############ */
1245 /* With Leap-Frog we can skip compute_globals at
1246 * non-communication steps, but we need to calculate
1247 * the kinetic energy one step before communication.
1250 // Organize to do inter-simulation signalling on steps if
1251 // and when algorithms require it.
1252 bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1254 if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1256 // Since we're already communicating at this step, we
1257 // can propagate intra-simulation signals. Note that
1258 // check_nstglobalcomm has the responsibility for
1259 // choosing the value of nstglobalcomm that is one way
1260 // bGStat becomes true, so we can't get into a
1261 // situation where e.g. checkpointing can't be
1263 bool doIntraSimSignal = true;
1264 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1266 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1267 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1270 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1271 (bGStat ? CGLO_GSTAT : 0)
1272 | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1273 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1274 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1275 | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
1277 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1279 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1280 top_global, &top, state,
1281 &shouldCheckNumberOfBondedInteractions);
1285 /* ############# END CALC EKIN AND PRESSURE ################# */
1287 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1288 the virial that should probably be addressed eventually. state->veta has better properies,
1289 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1290 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1292 if (ir->efep != efepNO && !EI_VV(ir->eI))
1294 /* Sum up the foreign energy and dhdl terms for md and sd.
1295 Currently done every step so that dhdl is correct in the .edr */
1296 sum_dhdl(enerd, state->lambda, ir->fepvals);
1299 update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1300 pres, force_vir, shake_vir,
1304 /* ################# END UPDATE STEP 2 ################# */
1305 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1307 /* The coordinates (x) were unshifted in update */
1310 /* We will not sum ekinh_old,
1311 * so signal that we still have to do it.
1313 bSumEkinhOld = TRUE;
1318 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1320 /* use the directly determined last velocity, not actually the averaged half steps */
1321 if (bTrotter && ir->eI == eiVV)
1323 enerd->term[F_EKIN] = last_ekin;
1325 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1327 if (integratorHasConservedEnergyQuantity(ir))
1331 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1335 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1338 /* ######### END PREPARING EDR OUTPUT ########### */
1344 if (fplog && do_log && bDoExpanded)
1346 /* only needed if doing expanded ensemble */
1347 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1348 state_global->dfhist, state->fep_state, ir->nstlog, step);
1352 energyOutput.addDataAtEnergyStep(bDoDHDL, bCalcEnerStep,
1353 t, mdatoms->tmass, enerd, state,
1354 ir->fepvals, ir->expandedvals, lastbox,
1355 shake_vir, force_vir, total_vir, pres,
1356 ekind, mu_tot, constr);
1360 energyOutput.recordNonEnergyStep();
1363 gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
1364 gmx_bool do_or = do_per_step(step, ir->nstorireout);
1366 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
1367 do_log ? fplog : nullptr,
1369 eprNORMAL, fcd, groups, &(ir->opts), awh.get());
1373 pull_print_output(ir->pull_work, step, t);
1376 if (do_per_step(step, ir->nstlog))
1378 if (fflush(fplog) != 0)
1380 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1386 /* Have to do this part _after_ outputting the logfile and the edr file */
1387 /* Gets written into the state at the beginning of next loop*/
1388 state->fep_state = lamnew;
1390 /* Print the remaining wall clock time for the run */
1391 if (isMasterSimMasterRank(ms, cr) &&
1392 (do_verbose || gmx_got_usr_signal()) &&
1397 fprintf(stderr, "\n");
1399 print_time(stderr, walltime_accounting, step, ir, cr);
1402 /* Ion/water position swapping.
1403 * Not done in last step since trajectory writing happens before this call
1404 * in the MD loop and exchanges would be lost anyway. */
1405 bNeedRepartition = FALSE;
1406 if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1407 do_per_step(step, ir->swap->nstswap))
1409 bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
1410 as_rvec_array(state->x.data()),
1412 MASTER(cr) && mdrunOptions.verbose,
1415 if (bNeedRepartition && DOMAINDECOMP(cr))
1417 dd_collect_state(cr->dd, state, state_global);
1421 /* Replica exchange */
1425 bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
1426 state_global, enerd,
1430 if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1432 dd_partition_system(fplog, mdlog, step, cr, TRUE, 1,
1433 state_global, *top_global, ir,
1434 state, &f, mdAtoms, &top, fr,
1436 nrnb, wcycle, FALSE);
1437 shouldCheckNumberOfBondedInteractions = true;
1438 upd.setNumAtoms(state->natoms);
1443 startingFromCheckpoint = false;
1445 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1446 /* With all integrators, except VV, we need to retain the pressure
1447 * at the current step for coupling at the next step.
1449 if ((state->flags & (1<<estPRES_PREV)) &&
1451 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1453 /* Store the pressure in t_state for pressure coupling
1454 * at the next MD step.
1456 copy_mat(pres, state->pres_prev);
1459 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1461 if ( (membed != nullptr) && (!bLastStep) )
1463 rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1466 cycles = wallcycle_stop(wcycle, ewcSTEP);
1467 if (DOMAINDECOMP(cr) && wcycle)
1469 dd_cycles_add(cr->dd, cycles, ddCyclStep);
1472 /* increase the MD step number */
1476 resetHandler->resetCounters(
1477 step, step_rel, mdlog, fplog, cr, (use_GPU(fr->nbv) ? fr->nbv : nullptr),
1478 nrnb, fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
1480 /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1481 IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
1484 /* End of main MD loop */
1486 /* Closing TNG files can include compressing data. Therefore it is good to do that
1487 * before stopping the time measurements. */
1488 mdoutf_tng_close(outf);
1490 /* Stop measuring walltime */
1491 walltime_accounting_end_time(walltime_accounting);
1493 if (!thisRankHasDuty(cr, DUTY_PME))
1495 /* Tell the PME only node to finish */
1496 gmx_pme_send_finish(cr);
1501 if (ir->nstcalcenergy > 0)
1503 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE,
1505 eprAVER, fcd, groups, &(ir->opts), awh.get());
1512 pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
1515 done_shellfc(fplog, shellfc, step_rel);
1517 if (useReplicaExchange && MASTER(cr))
1519 print_replica_exchange_statistics(fplog, repl_ex);
1522 // Clean up swapcoords
1523 if (ir->eSwapCoords != eswapNO)
1525 finish_swapcoords(ir->swap);
1528 /* IMD cleanup, if bIMD is TRUE. */
1529 IMD_finalize(ir->bIMD, ir->imd);
1531 walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1533 destroy_enerdata(enerd);
1536 global_stat_destroy(gstat);