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/mdrunutility/printtime.h"
103 #include "gromacs/mdtypes/awh_history.h"
104 #include "gromacs/mdtypes/awh_params.h"
105 #include "gromacs/mdtypes/commrec.h"
106 #include "gromacs/mdtypes/df_history.h"
107 #include "gromacs/mdtypes/energyhistory.h"
108 #include "gromacs/mdtypes/fcdata.h"
109 #include "gromacs/mdtypes/forcerec.h"
110 #include "gromacs/mdtypes/group.h"
111 #include "gromacs/mdtypes/inputrec.h"
112 #include "gromacs/mdtypes/interaction_const.h"
113 #include "gromacs/mdtypes/md_enums.h"
114 #include "gromacs/mdtypes/mdatom.h"
115 #include "gromacs/mdtypes/observableshistory.h"
116 #include "gromacs/mdtypes/pullhistory.h"
117 #include "gromacs/mdtypes/state.h"
118 #include "gromacs/nbnxm/nbnxm.h"
119 #include "gromacs/pbcutil/mshift.h"
120 #include "gromacs/pbcutil/pbc.h"
121 #include "gromacs/pulling/output.h"
122 #include "gromacs/pulling/pull.h"
123 #include "gromacs/swap/swapcoords.h"
124 #include "gromacs/timing/wallcycle.h"
125 #include "gromacs/timing/walltime_accounting.h"
126 #include "gromacs/topology/atoms.h"
127 #include "gromacs/topology/idef.h"
128 #include "gromacs/topology/mtop_util.h"
129 #include "gromacs/topology/topology.h"
130 #include "gromacs/trajectory/trajectoryframe.h"
131 #include "gromacs/utility/basedefinitions.h"
132 #include "gromacs/utility/cstringutil.h"
133 #include "gromacs/utility/fatalerror.h"
134 #include "gromacs/utility/logger.h"
135 #include "gromacs/utility/real.h"
136 #include "gromacs/utility/smalloc.h"
138 #include "integrator.h"
139 #include "replicaexchange.h"
142 #include "corewrap.h"
145 using gmx::SimulationSignaller;
147 void gmx::Integrator::do_md()
149 // TODO Historically, the EM and MD "integrators" used different
150 // names for the t_inputrec *parameter, but these must have the
151 // same name, now that it's a member of a struct. We use this ir
152 // alias to avoid a large ripple of nearly useless changes.
153 // t_inputrec is being replaced by IMdpOptionsProvider, so this
154 // will go away eventually.
155 t_inputrec *ir = inputrec;
156 int64_t step, step_rel;
157 double t = ir->init_t, t0 = ir->init_t, lam0[efptNR];
158 gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
159 gmx_bool bNS, bNStList, bStopCM,
160 bFirstStep, bInitStep, bLastStep = FALSE;
161 gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
162 gmx_bool do_ene, do_log, do_verbose;
163 gmx_bool bMasterState;
164 int force_flags, cglo_flags;
165 tensor force_vir = {{0}}, shake_vir = {{0}}, total_vir = {{0}},
166 tmp_vir = {{0}}, pres = {{0}};
169 matrix parrinellorahmanMu, M;
170 gmx_repl_ex_t repl_ex = nullptr;
172 gmx_enerdata_t *enerd;
173 PaddedVector<gmx::RVec> f {};
174 gmx_global_stat_t gstat;
175 t_graph *graph = nullptr;
176 gmx_groups_t *groups;
177 gmx_shellfc_t *shellfc;
178 gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
179 gmx_bool bTemp, bPres, bTrotter;
181 rvec *cbuf = nullptr;
188 real saved_conserved_quantity = 0;
191 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
193 /* PME load balancing data for GPU kernels */
194 gmx_bool bPMETune = FALSE;
195 gmx_bool bPMETunePrinting = FALSE;
198 gmx_bool bIMDstep = FALSE;
200 /* Domain decomposition could incorrectly miss a bonded
201 interaction, but checking for that requires a global
202 communication stage, which does not otherwise happen in DD
203 code. So we do that alongside the first global energy reduction
204 after a new DD is made. These variables handle whether the
205 check happens, and the result it returns. */
206 bool shouldCheckNumberOfBondedInteractions = false;
207 int totalNumberOfBondedInteractions = -1;
209 SimulationSignals signals;
210 // Most global communnication stages don't propagate mdrun
211 // signals, and will use this object to achieve that.
212 SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
214 if (!mdrunOptions.writeConfout)
216 // This is on by default, and the main known use case for
217 // turning it off is for convenience in benchmarking, which is
218 // something that should not show up in the general user
220 GMX_LOG(mdlog.info).asParagraph().
221 appendText("The -noconfout functionality is deprecated, and may be removed in a future version.");
224 /* md-vv uses averaged full step velocities for T-control
225 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
226 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
227 bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
229 const bool bRerunMD = false;
230 int nstglobalcomm = mdrunOptions.globalCommunicationInterval;
232 nstglobalcomm = check_nstglobalcomm(mdlog, nstglobalcomm, ir, cr);
233 bGStatEveryStep = (nstglobalcomm == 1);
235 groups = &top_global->groups;
237 std::unique_ptr<EssentialDynamics> ed = nullptr;
238 if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
240 /* Initialize essential dynamics sampling */
241 ed = init_edsam(mdlog,
242 opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
245 state_global, observablesHistory,
246 oenv, mdrunOptions.continuationOptions.appendFiles);
249 initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda, lam0);
250 Update upd(ir, deform);
251 bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
252 if (!mdrunOptions.continuationOptions.appendFiles)
254 pleaseCiteCouplingAlgorithms(fplog, *ir);
257 gmx_mdoutf *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, ir, top_global, oenv, wcycle);
258 gmx::EnergyOutput energyOutput;
259 energyOutput.prepare(mdoutf_get_fp_ene(outf), top_global, ir, mdoutf_get_fp_dhdl(outf));
261 /* Energy terms and groups */
263 init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
266 /* Kinetic energy data */
267 std::unique_ptr<gmx_ekindata_t> eKinData = std::make_unique<gmx_ekindata_t>();
268 gmx_ekindata_t *ekind = eKinData.get();
269 init_ekindata(fplog, top_global, &(ir->opts), ekind);
270 /* Copy the cos acceleration to the groups struct */
271 ekind->cosacc.cos_accel = ir->cos_accel;
273 gstat = global_stat_init(ir);
275 /* Check for polarizable models and flexible constraints */
276 shellfc = init_shell_flexcon(fplog,
277 top_global, constr ? constr->numFlexibleConstraints() : 0,
278 ir->nstcalcenergy, DOMAINDECOMP(cr));
281 double io = compute_io(ir, top_global->natoms, groups, energyOutput.numEnergyTerms(), 1);
282 if ((io > 2000) && MASTER(cr))
285 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
290 /* Set up interactive MD (IMD) */
291 init_IMD(ir, cr, ms, top_global, fplog, ir->nstcalcenergy,
292 MASTER(cr) ? state_global->x.rvec_array() : nullptr,
293 nfile, fnm, oenv, mdrunOptions);
295 // Local state only becomes valid now.
296 std::unique_ptr<t_state> stateInstance;
299 if (DOMAINDECOMP(cr))
301 dd_init_local_top(*top_global, &top);
303 stateInstance = std::make_unique<t_state>();
304 state = stateInstance.get();
305 dd_init_local_state(cr->dd, state_global, state);
307 /* Distribute the charge groups over the nodes from the master node */
308 dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
309 state_global, *top_global, ir,
310 state, &f, mdAtoms, &top, fr,
312 nrnb, nullptr, FALSE);
313 shouldCheckNumberOfBondedInteractions = true;
314 upd.setNumAtoms(state->natoms);
318 state_change_natoms(state_global, state_global->natoms);
319 f.resizeWithPadding(state_global->natoms);
320 /* Copy the pointer to the global state */
321 state = state_global;
323 /* Generate and initialize new topology */
324 mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr,
325 &graph, mdAtoms, constr, vsite, shellfc);
327 upd.setNumAtoms(state->natoms);
330 auto mdatoms = mdAtoms->mdatoms();
332 // NOTE: The global state is no longer used at this point.
333 // But state_global is still used as temporary storage space for writing
334 // the global state to file and potentially for replica exchange.
335 // (Global topology should persist.)
337 update_mdatoms(mdatoms, state->lambda[efptMASS]);
339 const ContinuationOptions &continuationOptions = mdrunOptions.continuationOptions;
340 bool startingFromCheckpoint = continuationOptions.startedFromCheckpoint;
344 /* Check nstexpanded here, because the grompp check was broken */
345 if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
347 gmx_fatal(FARGS, "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
349 init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
354 if (startingFromCheckpoint)
356 /* Restore from energy history if appending to output files */
357 if (continuationOptions.appendFiles)
359 /* If no history is available (because a checkpoint is from before
360 * it was written) make a new one later, otherwise restore it.
362 if (observablesHistory->energyHistory)
364 energyOutput.restoreFromEnergyHistory(*observablesHistory->energyHistory);
367 else if (observablesHistory->energyHistory)
369 /* We might have read an energy history from checkpoint.
370 * As we are not appending, we want to restart the statistics.
371 * Free the allocated memory and reset the counts.
373 observablesHistory->energyHistory = {};
374 /* We might have read a pull history from checkpoint.
375 * We will still want to keep the statistics, so that the files
376 * can be joined and still be meaningful.
377 * This means that observablesHistory->pullHistory
378 * should not be reset.
382 if (!observablesHistory->energyHistory)
384 observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
386 if (!observablesHistory->pullHistory)
388 observablesHistory->pullHistory = std::make_unique<PullHistory>();
390 /* Set the initial energy history */
391 energyOutput.fillEnergyHistory(observablesHistory->energyHistory.get());
394 preparePrevStepPullCom(ir, mdatoms, state, state_global, cr, startingFromCheckpoint);
396 // TODO: Remove this by converting AWH into a ForceProvider
397 auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingFromCheckpoint,
399 opt2fn("-awh", nfile, fnm), ir->pull_work);
401 const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
402 if (useReplicaExchange && MASTER(cr))
404 repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
407 /* PME tuning is only supported in the Verlet scheme, with PME for
408 * Coulomb. It is not supported with only LJ PME. */
409 bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) &&
410 !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
412 pme_load_balancing_t *pme_loadbal = nullptr;
415 pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box,
416 *fr->ic, fr->nbv->pairlistSets().params(), fr->pmedata, use_GPU(fr->nbv.get()),
420 if (!ir->bContinuation)
422 if (state->flags & (1 << estV))
424 auto v = makeArrayRef(state->v);
425 /* Set the velocities of vsites, shells and frozen atoms to zero */
426 for (i = 0; i < mdatoms->homenr; i++)
428 if (mdatoms->ptype[i] == eptVSite ||
429 mdatoms->ptype[i] == eptShell)
433 else if (mdatoms->cFREEZE)
435 for (m = 0; m < DIM; m++)
437 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
448 /* Constrain the initial coordinates and velocities */
449 do_constrain_first(fplog, constr, ir, mdatoms, state);
453 /* Construct the virtual sites for the initial configuration */
454 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
455 top.idef.iparams, top.idef.il,
456 fr->ePBC, fr->bMolPBC, cr, state->box);
460 if (ir->efep != efepNO)
462 /* Set free energy calculation frequency as the greatest common
463 * denominator of nstdhdl and repl_ex_nst. */
464 nstfep = ir->fepvals->nstdhdl;
467 nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
469 if (useReplicaExchange)
471 nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
475 /* Be REALLY careful about what flags you set here. You CANNOT assume
476 * this is the first step, since we might be restarting from a checkpoint,
477 * and in that case we should not do any modifications to the state.
479 bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
481 if (continuationOptions.haveReadEkin)
483 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
486 cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT
487 | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
488 | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
489 | (continuationOptions.haveReadEkin ? CGLO_READEKIN : 0));
491 bSumEkinhOld = FALSE;
493 t_vcm vcm(top_global->groups, *ir);
494 reportComRemovalInfo(fplog, vcm);
496 /* To minimize communication, compute_globals computes the COM velocity
497 * and the kinetic energy for the velocities without COM motion removed.
498 * Thus to get the kinetic energy without the COM contribution, we need
499 * to call compute_globals twice.
501 for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
503 int cglo_flags_iteration = cglo_flags;
504 if (bStopCM && cgloIteration == 0)
506 cglo_flags_iteration |= CGLO_STOPCM;
507 cglo_flags_iteration &= ~CGLO_TEMPERATURE;
509 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
510 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
511 constr, &nullSignaller, state->box,
512 &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
513 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
515 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
516 top_global, &top, state,
517 &shouldCheckNumberOfBondedInteractions);
518 if (ir->eI == eiVVAK)
520 /* a second call to get the half step temperature initialized as well */
521 /* we do the same call as above, but turn the pressure off -- internally to
522 compute_globals, this is recognized as a velocity verlet half-step
523 kinetic energy calculation. This minimized excess variables, but
524 perhaps loses some logic?*/
526 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
527 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
528 constr, &nullSignaller, state->box,
529 nullptr, &bSumEkinhOld,
530 cglo_flags & ~CGLO_PRESSURE);
533 /* Calculate the initial half step temperature, and save the ekinh_old */
534 if (!continuationOptions.startedFromCheckpoint)
536 for (i = 0; (i < ir->opts.ngtc); i++)
538 copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
542 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
543 temperature control */
544 auto trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
548 if (!ir->bContinuation)
550 if (constr && ir->eConstrAlg == econtLINCS)
553 "RMS relative constraint deviation after constraining: %.2e\n",
556 if (EI_STATE_VELOCITY(ir->eI))
558 real temp = enerd->term[F_TEMP];
561 /* Result of Ekin averaged over velocities of -half
562 * and +half step, while we only have -half step here.
566 fprintf(fplog, "Initial temperature: %g K\n", temp);
571 fprintf(stderr, "starting mdrun '%s'\n",
572 *(top_global->name));
575 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
579 sprintf(tbuf, "%s", "infinite");
581 if (ir->init_step > 0)
583 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
584 gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
585 gmx_step_str(ir->init_step, sbuf2),
586 ir->init_step*ir->delta_t);
590 fprintf(stderr, "%s steps, %s ps.\n",
591 gmx_step_str(ir->nsteps, sbuf), tbuf);
593 fprintf(fplog, "\n");
596 walltime_accounting_start_time(walltime_accounting);
597 wallcycle_start(wcycle, ewcRUN);
598 print_start(fplog, cr, walltime_accounting, "mdrun");
601 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
602 int chkpt_ret = fcCheckPointParallel( cr->nodeid,
606 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
610 /***********************************************************
614 ************************************************************/
617 /* Skip the first Nose-Hoover integration when we get the state from tpx */
618 bInitStep = !startingFromCheckpoint || EI_VV(ir->eI);
619 bSumEkinhOld = FALSE;
621 bNeedRepartition = FALSE;
623 bool simulationsShareState = false;
624 int nstSignalComm = nstglobalcomm;
626 // TODO This implementation of ensemble orientation restraints is nasty because
627 // a user can't just do multi-sim with single-sim orientation restraints.
628 bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
629 bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
631 // Replica exchange, ensemble restraints and AWH need all
632 // simulations to remain synchronized, so they need
633 // checkpoints and stop conditions to act on the same step, so
634 // the propagation of such signals must take place between
635 // simulations, not just within simulations.
636 // TODO: Make algorithm initializers set these flags.
637 simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
639 if (simulationsShareState)
641 // Inter-simulation signal communication does not need to happen
642 // often, so we use a minimum of 200 steps to reduce overhead.
643 const int c_minimumInterSimulationSignallingInterval = 200;
644 nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
648 auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
649 compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
650 MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
651 mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
653 auto checkpointHandler = std::make_unique<CheckpointHandler>(
654 compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]),
655 simulationsShareState, ir->nstlist == 0, MASTER(cr),
656 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
658 const bool resetCountersIsLocal = true;
659 auto resetHandler = std::make_unique<ResetHandler>(
660 compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]), !resetCountersIsLocal,
661 ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
662 mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
664 const DDBalanceRegionHandler ddBalanceRegionHandler(cr);
666 step = ir->init_step;
669 // TODO extract this to new multi-simulation module
670 if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
672 if (!multisim_int_all_are_equal(ms, ir->nsteps))
674 GMX_LOG(mdlog.warning).appendText(
675 "Note: The number of steps is not consistent across multi simulations,\n"
676 "but we are proceeding anyway!");
678 if (!multisim_int_all_are_equal(ms, ir->init_step))
680 GMX_LOG(mdlog.warning).appendText(
681 "Note: The initial step is not consistent across multi simulations,\n"
682 "but we are proceeding anyway!");
686 /* and stop now if we should */
687 bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
691 /* Determine if this is a neighbor search step */
692 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
694 if (bPMETune && bNStList)
696 /* PME grid + cut-off optimization with GPUs or PME nodes */
697 pme_loadbal_do(pme_loadbal, cr,
698 (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
706 wallcycle_start(wcycle, ewcSTEP);
708 bLastStep = (step_rel == ir->nsteps);
709 t = t0 + step*ir->delta_t;
711 // TODO Refactor this, so that nstfep does not need a default value of zero
712 if (ir->efep != efepNO || ir->bSimTemp)
714 /* find and set the current lambdas */
715 setCurrentLambdasLocal(step, ir->fepvals, lam0, state);
717 bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
718 bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
719 bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
720 && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
723 bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
724 do_per_step(step, replExParams.exchangeInterval));
726 if (doSimulatedAnnealing)
728 update_annealing_target_temp(ir, t, &upd);
731 /* Stop Center of Mass motion */
732 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
734 /* Determine whether or not to do Neighbour Searching */
735 bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
737 bLastStep = bLastStep || stopHandler->stoppingAfterCurrentStep(bNS);
739 /* do_log triggers energy and virial calculation. Because this leads
740 * to different code paths, forces can be different. Thus for exact
741 * continuation we should avoid extra log output.
742 * Note that the || bLastStep can result in non-exact continuation
743 * beyond the last step. But we don't consider that to be an issue.
745 do_log = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep;
746 do_verbose = mdrunOptions.verbose &&
747 (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
749 if (bNS && !(bFirstStep && ir->bContinuation))
751 bMasterState = FALSE;
752 /* Correct the new box if it is too skewed */
753 if (inputrecDynamicBox(ir))
755 if (correct_box(fplog, step, state->box, graph))
760 if (DOMAINDECOMP(cr) && bMasterState)
762 dd_collect_state(cr->dd, state, state_global);
765 if (DOMAINDECOMP(cr))
767 /* Repartition the domain decomposition */
768 dd_partition_system(fplog, mdlog, step, cr,
769 bMasterState, nstglobalcomm,
770 state_global, *top_global, ir,
771 state, &f, mdAtoms, &top, fr,
774 do_verbose && !bPMETunePrinting);
775 shouldCheckNumberOfBondedInteractions = true;
776 upd.setNumAtoms(state->natoms);
780 if (MASTER(cr) && do_log)
782 print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
785 if (ir->efep != efepNO)
787 update_mdatoms(mdatoms, state->lambda[efptMASS]);
793 /* We need the kinetic energy at minus the half step for determining
794 * the full step kinetic energy and possibly for T-coupling.*/
795 /* This may not be quite working correctly yet . . . . */
796 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
797 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
798 constr, &nullSignaller, state->box,
799 &totalNumberOfBondedInteractions, &bSumEkinhOld,
800 CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
801 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
802 top_global, &top, state,
803 &shouldCheckNumberOfBondedInteractions);
805 clear_mat(force_vir);
807 checkpointHandler->decideIfCheckpointingThisStep(bNS, bFirstStep, bLastStep);
809 /* Determine the energy and pressure:
810 * at nstcalcenergy steps and at energy output steps (set below).
812 if (EI_VV(ir->eI) && (!bInitStep))
814 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
815 bCalcVir = bCalcEnerStep ||
816 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
820 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
821 bCalcVir = bCalcEnerStep ||
822 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
824 bCalcEner = bCalcEnerStep;
826 do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
828 if (do_ene || do_log || bDoReplEx)
834 /* Do we need global communication ? */
835 bGStat = (bCalcVir || bCalcEner || bStopCM ||
836 do_per_step(step, nstglobalcomm) ||
837 (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
839 force_flags = (GMX_FORCE_STATECHANGED |
840 ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0) |
841 GMX_FORCE_ALLFORCES |
842 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
843 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
844 (bDoFEP ? GMX_FORCE_DHDL : 0)
849 /* Now is the time to relax the shells */
850 relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
851 enforcedRotation, step,
852 ir, bNS, force_flags, &top,
854 state, f.arrayRefWithPadding(), force_vir, mdatoms,
855 nrnb, wcycle, graph, groups,
856 shellfc, fr, ppForceWorkload, t, mu_tot,
858 ddBalanceRegionHandler);
862 /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
863 (or the AWH update will be performed twice for one step when continuing). It would be best to
864 call this update function from do_md_trajectory_writing but that would occur after do_force.
865 One would have to divide the update_awh function into one function applying the AWH force
866 and one doing the AWH bias update. The update AWH bias function could then be called after
867 do_md_trajectory_writing (then containing update_awh_history).
868 The checkpointing will in the future probably moved to the start of the md loop which will
869 rid of this issue. */
870 if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
872 awh->updateHistory(state_global->awhHistory.get());
875 /* The coordinates (x) are shifted (to get whole molecules)
877 * This is parallellized as well, and does communication too.
878 * Check comments in sim_util.c
880 do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation,
881 step, nrnb, wcycle, &top, groups,
882 state->box, state->x.arrayRefWithPadding(), &state->hist,
883 f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
884 state->lambda, graph,
885 fr, ppForceWorkload, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
886 (bNS ? GMX_FORCE_NS : 0) | force_flags,
887 ddBalanceRegionHandler);
890 if (EI_VV(ir->eI) && !startingFromCheckpoint)
891 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
893 rvec *vbuf = nullptr;
895 wallcycle_start(wcycle, ewcUPDATE);
896 if (ir->eI == eiVV && bInitStep)
898 /* if using velocity verlet with full time step Ekin,
899 * take the first half step only to compute the
900 * virial for the first step. From there,
901 * revert back to the initial coordinates
902 * so that the input is actually the initial step.
904 snew(vbuf, state->natoms);
905 copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
909 /* this is for NHC in the Ekin(t+dt/2) version of vv */
910 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
913 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
914 ekind, M, &upd, etrtVELOCITY1,
917 wallcycle_stop(wcycle, ewcUPDATE);
918 constrain_velocities(step, nullptr,
922 bCalcVir, do_log, do_ene);
923 wallcycle_start(wcycle, ewcUPDATE);
924 /* if VV, compute the pressure and constraints */
925 /* For VV2, we strictly only need this if using pressure
926 * control, but we really would like to have accurate pressures
928 * Think about ways around this in the future?
929 * For now, keep this choice in comments.
931 /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
932 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
934 bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
935 if (bCalcEner && ir->eI == eiVVAK)
939 /* for vv, the first half of the integration actually corresponds to the previous step.
940 So we need information from the last step in the first half of the integration */
941 if (bGStat || do_per_step(step-1, nstglobalcomm))
943 wallcycle_stop(wcycle, ewcUPDATE);
944 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
945 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
946 constr, &nullSignaller, state->box,
947 &totalNumberOfBondedInteractions, &bSumEkinhOld,
948 (bGStat ? CGLO_GSTAT : 0)
949 | (bCalcEner ? CGLO_ENERGY : 0)
950 | (bTemp ? CGLO_TEMPERATURE : 0)
951 | (bPres ? CGLO_PRESSURE : 0)
952 | (bPres ? CGLO_CONSTRAINT : 0)
953 | (bStopCM ? CGLO_STOPCM : 0)
954 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
957 /* explanation of above:
958 a) We compute Ekin at the full time step
959 if 1) we are using the AveVel Ekin, and it's not the
960 initial step, or 2) if we are using AveEkin, but need the full
961 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
962 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
963 EkinAveVel because it's needed for the pressure */
964 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
965 top_global, &top, state,
966 &shouldCheckNumberOfBondedInteractions);
967 wallcycle_start(wcycle, ewcUPDATE);
969 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
974 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
975 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
977 /* TODO This is only needed when we're about to write
978 * a checkpoint, because we use it after the restart
979 * (in a kludge?). But what should we be doing if
980 * startingFromCheckpoint or bInitStep are true? */
981 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
983 copy_mat(shake_vir, state->svir_prev);
984 copy_mat(force_vir, state->fvir_prev);
986 if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
988 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
989 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
990 enerd->term[F_EKIN] = trace(ekind->ekin);
995 wallcycle_stop(wcycle, ewcUPDATE);
996 /* We need the kinetic energy at minus the half step for determining
997 * the full step kinetic energy and possibly for T-coupling.*/
998 /* This may not be quite working correctly yet . . . . */
999 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1000 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1001 constr, &nullSignaller, state->box,
1002 nullptr, &bSumEkinhOld,
1003 CGLO_GSTAT | CGLO_TEMPERATURE);
1004 wallcycle_start(wcycle, ewcUPDATE);
1007 /* if it's the initial step, we performed this first step just to get the constraint virial */
1008 if (ir->eI == eiVV && bInitStep)
1010 copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
1013 wallcycle_stop(wcycle, ewcUPDATE);
1016 /* compute the conserved quantity */
1019 saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1022 last_ekin = enerd->term[F_EKIN];
1024 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1026 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1028 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1029 if (ir->efep != efepNO)
1031 sum_dhdl(enerd, state->lambda, ir->fepvals);
1035 /* ######## END FIRST UPDATE STEP ############## */
1036 /* ######## If doing VV, we now have v(dt) ###### */
1039 /* perform extended ensemble sampling in lambda - we don't
1040 actually move to the new state before outputting
1041 statistics, but if performing simulated tempering, we
1042 do update the velocities and the tau_t. */
1044 lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v.rvec_array(), mdatoms);
1045 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1048 copy_df_history(state_global->dfhist, state->dfhist);
1052 /* Now we have the energies and forces corresponding to the
1053 * coordinates at time t. We must output all of this before
1056 do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1057 ir, state, state_global, observablesHistory,
1059 outf, energyOutput, ekind, f,
1060 checkpointHandler->isCheckpointingStep(),
1061 bRerunMD, bLastStep,
1062 mdrunOptions.writeConfout,
1064 /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1065 bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x.rvec_array(), ir, t, wcycle);
1067 /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1068 if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1070 copy_mat(state->svir_prev, shake_vir);
1071 copy_mat(state->fvir_prev, force_vir);
1074 stopHandler->setSignal();
1075 resetHandler->setSignal(walltime_accounting);
1077 if (bGStat || !PAR(cr))
1079 /* In parallel we only have to check for checkpointing in steps
1080 * where we do global communication,
1081 * otherwise the other nodes don't know.
1083 checkpointHandler->setSignal(walltime_accounting);
1086 /* ######### START SECOND UPDATE STEP ################# */
1088 /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1091 if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1093 gmx_bool bIfRandomize;
1094 bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
1095 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1096 if (constr && bIfRandomize)
1098 constrain_velocities(step, nullptr,
1102 bCalcVir, do_log, do_ene);
1105 /* Box is changed in update() when we do pressure coupling,
1106 * but we should still use the old box for energy corrections and when
1107 * writing it to the energy file, so it matches the trajectory files for
1108 * the same timestep above. Make a copy in a separate array.
1110 copy_mat(state->box, lastbox);
1114 wallcycle_start(wcycle, ewcUPDATE);
1115 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1118 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1119 /* We can only do Berendsen coupling after we have summed
1120 * the kinetic energy or virial. Since the happens
1121 * in global_state after update, we should only do it at
1122 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1127 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1128 update_pcouple_before_coordinates(fplog, step, ir, state,
1129 parrinellorahmanMu, M,
1135 /* velocity half-step update */
1136 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1137 ekind, M, &upd, etrtVELOCITY2,
1141 /* Above, initialize just copies ekinh into ekin,
1142 * it doesn't copy position (for VV),
1143 * and entire integrator for MD.
1146 if (ir->eI == eiVVAK)
1148 /* We probably only need md->homenr, not state->natoms */
1149 if (state->natoms > cbuf_nalloc)
1151 cbuf_nalloc = state->natoms;
1152 srenew(cbuf, cbuf_nalloc);
1154 copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1157 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1158 ekind, M, &upd, etrtPOSITION, cr, constr);
1159 wallcycle_stop(wcycle, ewcUPDATE);
1161 constrain_coordinates(step, &dvdl_constr, state,
1164 bCalcVir, do_log, do_ene);
1165 update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
1166 cr, nrnb, wcycle, &upd, constr, do_log, do_ene);
1167 finish_update(ir, mdatoms,
1169 nrnb, wcycle, &upd, constr);
1171 if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1173 updatePrevStepPullCom(ir->pull_work, state);
1176 if (ir->eI == eiVVAK)
1178 /* erase F_EKIN and F_TEMP here? */
1179 /* just compute the kinetic energy at the half step to perform a trotter step */
1180 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1181 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1182 constr, &nullSignaller, lastbox,
1183 nullptr, &bSumEkinhOld,
1184 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1186 wallcycle_start(wcycle, ewcUPDATE);
1187 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1188 /* now we know the scaling, we can compute the positions again again */
1189 copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1191 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1192 ekind, M, &upd, etrtPOSITION, cr, constr);
1193 wallcycle_stop(wcycle, ewcUPDATE);
1195 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1196 /* are the small terms in the shake_vir here due
1197 * to numerical errors, or are they important
1198 * physically? I'm thinking they are just errors, but not completely sure.
1199 * For now, will call without actually constraining, constr=NULL*/
1200 finish_update(ir, mdatoms,
1202 nrnb, wcycle, &upd, nullptr);
1206 /* this factor or 2 correction is necessary
1207 because half of the constraint force is removed
1208 in the vv step, so we have to double it. See
1209 the Redmine issue #1255. It is not yet clear
1210 if the factor of 2 is exact, or just a very
1211 good approximation, and this will be
1212 investigated. The next step is to see if this
1213 can be done adding a dhdl contribution from the
1214 rattle step, but this is somewhat more
1215 complicated with the current code. Will be
1216 investigated, hopefully for 4.6.3. However,
1217 this current solution is much better than
1218 having it completely wrong.
1220 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1224 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1227 if (vsite != nullptr)
1229 wallcycle_start(wcycle, ewcVSITECONSTR);
1230 if (graph != nullptr)
1232 shift_self(graph, state->box, state->x.rvec_array());
1234 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
1235 top.idef.iparams, top.idef.il,
1236 fr->ePBC, fr->bMolPBC, cr, state->box);
1238 if (graph != nullptr)
1240 unshift_self(graph, state->box, state->x.rvec_array());
1242 wallcycle_stop(wcycle, ewcVSITECONSTR);
1245 /* ############## IF NOT VV, Calculate globals HERE ############ */
1246 /* With Leap-Frog we can skip compute_globals at
1247 * non-communication steps, but we need to calculate
1248 * the kinetic energy one step before communication.
1251 // Organize to do inter-simulation signalling on steps if
1252 // and when algorithms require it.
1253 bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1255 if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1257 // Since we're already communicating at this step, we
1258 // can propagate intra-simulation signals. Note that
1259 // check_nstglobalcomm has the responsibility for
1260 // choosing the value of nstglobalcomm that is one way
1261 // bGStat becomes true, so we can't get into a
1262 // situation where e.g. checkpointing can't be
1264 bool doIntraSimSignal = true;
1265 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1267 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1268 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1271 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1272 (bGStat ? CGLO_GSTAT : 0)
1273 | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1274 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1275 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1276 | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
1278 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1280 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1281 top_global, &top, state,
1282 &shouldCheckNumberOfBondedInteractions);
1286 /* ############# END CALC EKIN AND PRESSURE ################# */
1288 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1289 the virial that should probably be addressed eventually. state->veta has better properies,
1290 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1291 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1293 if (ir->efep != efepNO && !EI_VV(ir->eI))
1295 /* Sum up the foreign energy and dhdl terms for md and sd.
1296 Currently done every step so that dhdl is correct in the .edr */
1297 sum_dhdl(enerd, state->lambda, ir->fepvals);
1300 update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1301 pres, force_vir, shake_vir,
1305 /* ################# END UPDATE STEP 2 ################# */
1306 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1308 /* The coordinates (x) were unshifted in update */
1311 /* We will not sum ekinh_old,
1312 * so signal that we still have to do it.
1314 bSumEkinhOld = TRUE;
1319 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1321 /* use the directly determined last velocity, not actually the averaged half steps */
1322 if (bTrotter && ir->eI == eiVV)
1324 enerd->term[F_EKIN] = last_ekin;
1326 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1328 if (integratorHasConservedEnergyQuantity(ir))
1332 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1336 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1339 /* ######### END PREPARING EDR OUTPUT ########### */
1345 if (fplog && do_log && bDoExpanded)
1347 /* only needed if doing expanded ensemble */
1348 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1349 state_global->dfhist, state->fep_state, ir->nstlog, step);
1353 energyOutput.addDataAtEnergyStep(bDoDHDL, bCalcEnerStep,
1354 t, mdatoms->tmass, enerd, state,
1355 ir->fepvals, ir->expandedvals, lastbox,
1356 shake_vir, force_vir, total_vir, pres,
1357 ekind, mu_tot, constr);
1361 energyOutput.recordNonEnergyStep();
1364 gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
1365 gmx_bool do_or = do_per_step(step, ir->nstorireout);
1367 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
1368 do_log ? fplog : nullptr,
1370 eprNORMAL, fcd, groups, &(ir->opts), awh.get());
1374 pull_print_output(ir->pull_work, step, t);
1377 if (do_per_step(step, ir->nstlog))
1379 if (fflush(fplog) != 0)
1381 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1387 /* Have to do this part _after_ outputting the logfile and the edr file */
1388 /* Gets written into the state at the beginning of next loop*/
1389 state->fep_state = lamnew;
1391 /* Print the remaining wall clock time for the run */
1392 if (isMasterSimMasterRank(ms, cr) &&
1393 (do_verbose || gmx_got_usr_signal()) &&
1398 fprintf(stderr, "\n");
1400 print_time(stderr, walltime_accounting, step, ir, cr);
1403 /* Ion/water position swapping.
1404 * Not done in last step since trajectory writing happens before this call
1405 * in the MD loop and exchanges would be lost anyway. */
1406 bNeedRepartition = FALSE;
1407 if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1408 do_per_step(step, ir->swap->nstswap))
1410 bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
1411 as_rvec_array(state->x.data()),
1413 MASTER(cr) && mdrunOptions.verbose,
1416 if (bNeedRepartition && DOMAINDECOMP(cr))
1418 dd_collect_state(cr->dd, state, state_global);
1422 /* Replica exchange */
1426 bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
1427 state_global, enerd,
1431 if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1433 dd_partition_system(fplog, mdlog, step, cr, TRUE, 1,
1434 state_global, *top_global, ir,
1435 state, &f, mdAtoms, &top, fr,
1437 nrnb, wcycle, FALSE);
1438 shouldCheckNumberOfBondedInteractions = true;
1439 upd.setNumAtoms(state->natoms);
1444 startingFromCheckpoint = false;
1446 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1447 /* With all integrators, except VV, we need to retain the pressure
1448 * at the current step for coupling at the next step.
1450 if ((state->flags & (1<<estPRES_PREV)) &&
1452 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1454 /* Store the pressure in t_state for pressure coupling
1455 * at the next MD step.
1457 copy_mat(pres, state->pres_prev);
1460 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1462 if ( (membed != nullptr) && (!bLastStep) )
1464 rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1467 cycles = wallcycle_stop(wcycle, ewcSTEP);
1468 if (DOMAINDECOMP(cr) && wcycle)
1470 dd_cycles_add(cr->dd, cycles, ddCyclStep);
1473 /* increase the MD step number */
1477 resetHandler->resetCounters(
1478 step, step_rel, mdlog, fplog, cr, (use_GPU(fr->nbv.get()) ? fr->nbv.get() : nullptr),
1479 nrnb, fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
1481 /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1482 IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
1485 /* End of main MD loop */
1487 /* Closing TNG files can include compressing data. Therefore it is good to do that
1488 * before stopping the time measurements. */
1489 mdoutf_tng_close(outf);
1491 /* Stop measuring walltime */
1492 walltime_accounting_end_time(walltime_accounting);
1494 if (!thisRankHasDuty(cr, DUTY_PME))
1496 /* Tell the PME only node to finish */
1497 gmx_pme_send_finish(cr);
1502 if (ir->nstcalcenergy > 0)
1504 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE,
1506 eprAVER, fcd, groups, &(ir->opts), awh.get());
1513 pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv.get()));
1516 done_shellfc(fplog, shellfc, step_rel);
1518 if (useReplicaExchange && MASTER(cr))
1520 print_replica_exchange_statistics(fplog, repl_ex);
1523 // Clean up swapcoords
1524 if (ir->eSwapCoords != eswapNO)
1526 finish_swapcoords(ir->swap);
1529 /* IMD cleanup, if bIMD is TRUE. */
1530 IMD_finalize(ir->bIMD, ir->imd);
1532 walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1534 destroy_enerdata(enerd);
1537 global_stat_destroy(gstat);