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/domdec.h"
58 #include "gromacs/domdec/domdec_network.h"
59 #include "gromacs/domdec/domdec_struct.h"
60 #include "gromacs/domdec/partition.h"
61 #include "gromacs/essentialdynamics/edsam.h"
62 #include "gromacs/ewald/pme.h"
63 #include "gromacs/ewald/pme-load-balancing.h"
64 #include "gromacs/fileio/trxio.h"
65 #include "gromacs/gmxlib/network.h"
66 #include "gromacs/gmxlib/nrnb.h"
67 #include "gromacs/gpu_utils/gpu_utils.h"
68 #include "gromacs/imd/imd.h"
69 #include "gromacs/listed_forces/manage-threading.h"
70 #include "gromacs/math/functions.h"
71 #include "gromacs/math/utilities.h"
72 #include "gromacs/math/vec.h"
73 #include "gromacs/math/vectypes.h"
74 #include "gromacs/mdlib/checkpointhandler.h"
75 #include "gromacs/mdlib/compute_io.h"
76 #include "gromacs/mdlib/constr.h"
77 #include "gromacs/mdlib/ebin.h"
78 #include "gromacs/mdlib/expanded.h"
79 #include "gromacs/mdlib/force.h"
80 #include "gromacs/mdlib/force_flags.h"
81 #include "gromacs/mdlib/forcerec.h"
82 #include "gromacs/mdlib/md_support.h"
83 #include "gromacs/mdlib/mdatoms.h"
84 #include "gromacs/mdlib/mdebin.h"
85 #include "gromacs/mdlib/mdoutf.h"
86 #include "gromacs/mdlib/mdrun.h"
87 #include "gromacs/mdlib/mdsetup.h"
88 #include "gromacs/mdlib/membed.h"
89 #include "gromacs/mdlib/nb_verlet.h"
90 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
91 #include "gromacs/mdlib/ns.h"
92 #include "gromacs/mdlib/resethandler.h"
93 #include "gromacs/mdlib/shellfc.h"
94 #include "gromacs/mdlib/sighandler.h"
95 #include "gromacs/mdlib/sim_util.h"
96 #include "gromacs/mdlib/simulationsignal.h"
97 #include "gromacs/mdlib/stophandler.h"
98 #include "gromacs/mdlib/tgroup.h"
99 #include "gromacs/mdlib/trajectory_writing.h"
100 #include "gromacs/mdlib/update.h"
101 #include "gromacs/mdlib/vcm.h"
102 #include "gromacs/mdlib/vsite.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/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,
257 outputProvider, ir, top_global, oenv, wcycle);
258 t_mdebin *mdebin = init_mdebin(mdrunOptions.continuationOptions.appendFiles ? nullptr : mdoutf_get_fp_ene(outf),
259 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, mdebin->ebin->nener, 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 /* Update mdebin with 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 restore_energyhistory_from_state(mdebin, observablesHistory->energyHistory.get());
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 in state by updating once */
391 update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
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->listParams, fr->pmedata, use_GPU(fr->nbv),
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 DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
665 DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion = (DOMAINDECOMP(cr) ? DdCloseBalanceRegionAfterForceComputation::yes : DdCloseBalanceRegionAfterForceComputation::no);
667 step = ir->init_step;
670 // TODO extract this to new multi-simulation module
671 if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
673 if (!multisim_int_all_are_equal(ms, ir->nsteps))
675 GMX_LOG(mdlog.warning).appendText(
676 "Note: The number of steps is not consistent across multi simulations,\n"
677 "but we are proceeding anyway!");
679 if (!multisim_int_all_are_equal(ms, ir->init_step))
681 GMX_LOG(mdlog.warning).appendText(
682 "Note: The initial step is not consistent across multi simulations,\n"
683 "but we are proceeding anyway!");
687 /* and stop now if we should */
688 bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
692 /* Determine if this is a neighbor search step */
693 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
695 if (bPMETune && bNStList)
697 /* PME grid + cut-off optimization with GPUs or PME nodes */
698 pme_loadbal_do(pme_loadbal, cr,
699 (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
707 wallcycle_start(wcycle, ewcSTEP);
709 bLastStep = (step_rel == ir->nsteps);
710 t = t0 + step*ir->delta_t;
712 // TODO Refactor this, so that nstfep does not need a default value of zero
713 if (ir->efep != efepNO || ir->bSimTemp)
715 /* find and set the current lambdas */
716 setCurrentLambdasLocal(step, ir->fepvals, lam0, state);
718 bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
719 bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
720 bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
721 && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
724 bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
725 do_per_step(step, replExParams.exchangeInterval));
727 if (doSimulatedAnnealing)
729 update_annealing_target_temp(ir, t, &upd);
732 /* Stop Center of Mass motion */
733 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
735 /* Determine whether or not to do Neighbour Searching */
736 bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
738 bLastStep = bLastStep || stopHandler->stoppingAfterCurrentStep(bNS);
740 /* do_log triggers energy and virial calculation. Because this leads
741 * to different code paths, forces can be different. Thus for exact
742 * continuation we should avoid extra log output.
743 * Note that the || bLastStep can result in non-exact continuation
744 * beyond the last step. But we don't consider that to be an issue.
746 do_log = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep;
747 do_verbose = mdrunOptions.verbose &&
748 (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
750 if (bNS && !(bFirstStep && ir->bContinuation))
752 bMasterState = FALSE;
753 /* Correct the new box if it is too skewed */
754 if (inputrecDynamicBox(ir))
756 if (correct_box(fplog, step, state->box, graph))
761 if (DOMAINDECOMP(cr) && bMasterState)
763 dd_collect_state(cr->dd, state, state_global);
766 if (DOMAINDECOMP(cr))
768 /* Repartition the domain decomposition */
769 dd_partition_system(fplog, mdlog, step, cr,
770 bMasterState, nstglobalcomm,
771 state_global, *top_global, ir,
772 state, &f, mdAtoms, &top, fr,
775 do_verbose && !bPMETunePrinting);
776 shouldCheckNumberOfBondedInteractions = true;
777 upd.setNumAtoms(state->natoms);
781 if (MASTER(cr) && do_log)
783 print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
786 if (ir->efep != efepNO)
788 update_mdatoms(mdatoms, state->lambda[efptMASS]);
794 /* We need the kinetic energy at minus the half step for determining
795 * the full step kinetic energy and possibly for T-coupling.*/
796 /* This may not be quite working correctly yet . . . . */
797 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
798 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
799 constr, &nullSignaller, state->box,
800 &totalNumberOfBondedInteractions, &bSumEkinhOld,
801 CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
802 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
803 top_global, &top, state,
804 &shouldCheckNumberOfBondedInteractions);
806 clear_mat(force_vir);
808 checkpointHandler->decideIfCheckpointingThisStep(bNS, bFirstStep, bLastStep);
810 /* Determine the energy and pressure:
811 * at nstcalcenergy steps and at energy output steps (set below).
813 if (EI_VV(ir->eI) && (!bInitStep))
815 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
816 bCalcVir = bCalcEnerStep ||
817 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
821 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
822 bCalcVir = bCalcEnerStep ||
823 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
825 bCalcEner = bCalcEnerStep;
827 do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
829 if (do_ene || do_log || bDoReplEx)
835 /* Do we need global communication ? */
836 bGStat = (bCalcVir || bCalcEner || bStopCM ||
837 do_per_step(step, nstglobalcomm) ||
838 (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
840 force_flags = (GMX_FORCE_STATECHANGED |
841 ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0) |
842 GMX_FORCE_ALLFORCES |
843 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
844 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
845 (bDoFEP ? GMX_FORCE_DHDL : 0)
850 /* Now is the time to relax the shells */
851 relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
852 enforcedRotation, step,
853 ir, bNS, force_flags, &top,
855 state, f.arrayRefWithPadding(), force_vir, mdatoms,
856 nrnb, wcycle, graph, groups,
857 shellfc, fr, ppForceWorkload, t, mu_tot,
859 ddOpenBalanceRegion, ddCloseBalanceRegion);
863 /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
864 (or the AWH update will be performed twice for one step when continuing). It would be best to
865 call this update function from do_md_trajectory_writing but that would occur after do_force.
866 One would have to divide the update_awh function into one function applying the AWH force
867 and one doing the AWH bias update. The update AWH bias function could then be called after
868 do_md_trajectory_writing (then containing update_awh_history).
869 The checkpointing will in the future probably moved to the start of the md loop which will
870 rid of this issue. */
871 if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
873 awh->updateHistory(state_global->awhHistory.get());
876 /* The coordinates (x) are shifted (to get whole molecules)
878 * This is parallellized as well, and does communication too.
879 * Check comments in sim_util.c
881 do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation,
882 step, nrnb, wcycle, &top, groups,
883 state->box, state->x.arrayRefWithPadding(), &state->hist,
884 f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
885 state->lambda, graph,
886 fr, ppForceWorkload, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
887 (bNS ? GMX_FORCE_NS : 0) | force_flags,
888 ddOpenBalanceRegion, ddCloseBalanceRegion);
891 if (EI_VV(ir->eI) && !startingFromCheckpoint)
892 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
894 rvec *vbuf = nullptr;
896 wallcycle_start(wcycle, ewcUPDATE);
897 if (ir->eI == eiVV && bInitStep)
899 /* if using velocity verlet with full time step Ekin,
900 * take the first half step only to compute the
901 * virial for the first step. From there,
902 * revert back to the initial coordinates
903 * so that the input is actually the initial step.
905 snew(vbuf, state->natoms);
906 copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
910 /* this is for NHC in the Ekin(t+dt/2) version of vv */
911 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
914 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
915 ekind, M, &upd, etrtVELOCITY1,
918 wallcycle_stop(wcycle, ewcUPDATE);
919 constrain_velocities(step, nullptr,
923 bCalcVir, do_log, do_ene);
924 wallcycle_start(wcycle, ewcUPDATE);
925 /* if VV, compute the pressure and constraints */
926 /* For VV2, we strictly only need this if using pressure
927 * control, but we really would like to have accurate pressures
929 * Think about ways around this in the future?
930 * For now, keep this choice in comments.
932 /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
933 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
935 bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
936 if (bCalcEner && ir->eI == eiVVAK)
940 /* for vv, the first half of the integration actually corresponds to the previous step.
941 So we need information from the last step in the first half of the integration */
942 if (bGStat || do_per_step(step-1, nstglobalcomm))
944 wallcycle_stop(wcycle, ewcUPDATE);
945 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
946 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
947 constr, &nullSignaller, state->box,
948 &totalNumberOfBondedInteractions, &bSumEkinhOld,
949 (bGStat ? CGLO_GSTAT : 0)
950 | (bCalcEner ? CGLO_ENERGY : 0)
951 | (bTemp ? CGLO_TEMPERATURE : 0)
952 | (bPres ? CGLO_PRESSURE : 0)
953 | (bPres ? CGLO_CONSTRAINT : 0)
954 | (bStopCM ? CGLO_STOPCM : 0)
955 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
958 /* explanation of above:
959 a) We compute Ekin at the full time step
960 if 1) we are using the AveVel Ekin, and it's not the
961 initial step, or 2) if we are using AveEkin, but need the full
962 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
963 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
964 EkinAveVel because it's needed for the pressure */
965 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
966 top_global, &top, state,
967 &shouldCheckNumberOfBondedInteractions);
968 wallcycle_start(wcycle, ewcUPDATE);
970 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
975 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
976 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
978 /* TODO This is only needed when we're about to write
979 * a checkpoint, because we use it after the restart
980 * (in a kludge?). But what should we be doing if
981 * startingFromCheckpoint or bInitStep are true? */
982 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
984 copy_mat(shake_vir, state->svir_prev);
985 copy_mat(force_vir, state->fvir_prev);
987 if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
989 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
990 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
991 enerd->term[F_EKIN] = trace(ekind->ekin);
996 wallcycle_stop(wcycle, ewcUPDATE);
997 /* We need the kinetic energy at minus the half step for determining
998 * the full step kinetic energy and possibly for T-coupling.*/
999 /* This may not be quite working correctly yet . . . . */
1000 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1001 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1002 constr, &nullSignaller, state->box,
1003 nullptr, &bSumEkinhOld,
1004 CGLO_GSTAT | CGLO_TEMPERATURE);
1005 wallcycle_start(wcycle, ewcUPDATE);
1008 /* if it's the initial step, we performed this first step just to get the constraint virial */
1009 if (ir->eI == eiVV && bInitStep)
1011 copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
1014 wallcycle_stop(wcycle, ewcUPDATE);
1017 /* compute the conserved quantity */
1020 saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1023 last_ekin = enerd->term[F_EKIN];
1025 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1027 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1029 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1030 if (ir->efep != efepNO)
1032 sum_dhdl(enerd, state->lambda, ir->fepvals);
1036 /* ######## END FIRST UPDATE STEP ############## */
1037 /* ######## If doing VV, we now have v(dt) ###### */
1040 /* perform extended ensemble sampling in lambda - we don't
1041 actually move to the new state before outputting
1042 statistics, but if performing simulated tempering, we
1043 do update the velocities and the tau_t. */
1045 lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v.rvec_array(), mdatoms);
1046 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1049 copy_df_history(state_global->dfhist, state->dfhist);
1053 /* Now we have the energies and forces corresponding to the
1054 * coordinates at time t. We must output all of this before
1057 do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1058 ir, state, state_global, observablesHistory,
1060 outf, mdebin, ekind, f,
1061 checkpointHandler->isCheckpointingStep(),
1062 bRerunMD, bLastStep,
1063 mdrunOptions.writeConfout,
1065 /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1066 bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x.rvec_array(), ir, t, wcycle);
1068 /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1069 if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1071 copy_mat(state->svir_prev, shake_vir);
1072 copy_mat(state->fvir_prev, force_vir);
1075 stopHandler->setSignal();
1076 resetHandler->setSignal(walltime_accounting);
1078 if (bGStat || !PAR(cr))
1080 /* In parallel we only have to check for checkpointing in steps
1081 * where we do global communication,
1082 * otherwise the other nodes don't know.
1084 checkpointHandler->setSignal(walltime_accounting);
1087 /* ######### START SECOND UPDATE STEP ################# */
1089 /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1092 if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1094 gmx_bool bIfRandomize;
1095 bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
1096 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1097 if (constr && bIfRandomize)
1099 constrain_velocities(step, nullptr,
1103 bCalcVir, do_log, do_ene);
1106 /* Box is changed in update() when we do pressure coupling,
1107 * but we should still use the old box for energy corrections and when
1108 * writing it to the energy file, so it matches the trajectory files for
1109 * the same timestep above. Make a copy in a separate array.
1111 copy_mat(state->box, lastbox);
1115 wallcycle_start(wcycle, ewcUPDATE);
1116 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1119 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1120 /* We can only do Berendsen coupling after we have summed
1121 * the kinetic energy or virial. Since the happens
1122 * in global_state after update, we should only do it at
1123 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1128 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1129 update_pcouple_before_coordinates(fplog, step, ir, state,
1130 parrinellorahmanMu, M,
1136 /* velocity half-step update */
1137 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1138 ekind, M, &upd, etrtVELOCITY2,
1142 /* Above, initialize just copies ekinh into ekin,
1143 * it doesn't copy position (for VV),
1144 * and entire integrator for MD.
1147 if (ir->eI == eiVVAK)
1149 /* We probably only need md->homenr, not state->natoms */
1150 if (state->natoms > cbuf_nalloc)
1152 cbuf_nalloc = state->natoms;
1153 srenew(cbuf, cbuf_nalloc);
1155 copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1158 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1159 ekind, M, &upd, etrtPOSITION, cr, constr);
1160 wallcycle_stop(wcycle, ewcUPDATE);
1162 constrain_coordinates(step, &dvdl_constr, state,
1165 bCalcVir, do_log, do_ene);
1166 update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
1167 cr, nrnb, wcycle, &upd, constr, do_log, do_ene);
1168 finish_update(ir, mdatoms,
1170 nrnb, wcycle, &upd, constr);
1172 if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1174 updatePrevStepPullCom(ir->pull_work, state);
1177 if (ir->eI == eiVVAK)
1179 /* erase F_EKIN and F_TEMP here? */
1180 /* just compute the kinetic energy at the half step to perform a trotter step */
1181 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1182 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1183 constr, &nullSignaller, lastbox,
1184 nullptr, &bSumEkinhOld,
1185 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1187 wallcycle_start(wcycle, ewcUPDATE);
1188 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1189 /* now we know the scaling, we can compute the positions again again */
1190 copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1192 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1193 ekind, M, &upd, etrtPOSITION, cr, constr);
1194 wallcycle_stop(wcycle, ewcUPDATE);
1196 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1197 /* are the small terms in the shake_vir here due
1198 * to numerical errors, or are they important
1199 * physically? I'm thinking they are just errors, but not completely sure.
1200 * For now, will call without actually constraining, constr=NULL*/
1201 finish_update(ir, mdatoms,
1203 nrnb, wcycle, &upd, nullptr);
1207 /* this factor or 2 correction is necessary
1208 because half of the constraint force is removed
1209 in the vv step, so we have to double it. See
1210 the Redmine issue #1255. It is not yet clear
1211 if the factor of 2 is exact, or just a very
1212 good approximation, and this will be
1213 investigated. The next step is to see if this
1214 can be done adding a dhdl contribution from the
1215 rattle step, but this is somewhat more
1216 complicated with the current code. Will be
1217 investigated, hopefully for 4.6.3. However,
1218 this current solution is much better than
1219 having it completely wrong.
1221 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1225 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1228 if (vsite != nullptr)
1230 wallcycle_start(wcycle, ewcVSITECONSTR);
1231 if (graph != nullptr)
1233 shift_self(graph, state->box, state->x.rvec_array());
1235 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
1236 top.idef.iparams, top.idef.il,
1237 fr->ePBC, fr->bMolPBC, cr, state->box);
1239 if (graph != nullptr)
1241 unshift_self(graph, state->box, state->x.rvec_array());
1243 wallcycle_stop(wcycle, ewcVSITECONSTR);
1246 /* ############## IF NOT VV, Calculate globals HERE ############ */
1247 /* With Leap-Frog we can skip compute_globals at
1248 * non-communication steps, but we need to calculate
1249 * the kinetic energy one step before communication.
1252 // Organize to do inter-simulation signalling on steps if
1253 // and when algorithms require it.
1254 bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1256 if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1258 // Since we're already communicating at this step, we
1259 // can propagate intra-simulation signals. Note that
1260 // check_nstglobalcomm has the responsibility for
1261 // choosing the value of nstglobalcomm that is one way
1262 // bGStat becomes true, so we can't get into a
1263 // situation where e.g. checkpointing can't be
1265 bool doIntraSimSignal = true;
1266 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1268 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1269 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1272 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1273 (bGStat ? CGLO_GSTAT : 0)
1274 | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1275 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1276 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1277 | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
1279 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1281 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1282 top_global, &top, state,
1283 &shouldCheckNumberOfBondedInteractions);
1287 /* ############# END CALC EKIN AND PRESSURE ################# */
1289 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1290 the virial that should probably be addressed eventually. state->veta has better properies,
1291 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1292 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1294 if (ir->efep != efepNO && !EI_VV(ir->eI))
1296 /* Sum up the foreign energy and dhdl terms for md and sd.
1297 Currently done every step so that dhdl is correct in the .edr */
1298 sum_dhdl(enerd, state->lambda, ir->fepvals);
1301 update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1302 pres, force_vir, shake_vir,
1306 /* ################# END UPDATE STEP 2 ################# */
1307 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1309 /* The coordinates (x) were unshifted in update */
1312 /* We will not sum ekinh_old,
1313 * so signal that we still have to do it.
1315 bSumEkinhOld = TRUE;
1320 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1322 /* use the directly determined last velocity, not actually the averaged half steps */
1323 if (bTrotter && ir->eI == eiVV)
1325 enerd->term[F_EKIN] = last_ekin;
1327 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1329 if (integratorHasConservedEnergyQuantity(ir))
1333 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1337 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1340 /* ######### END PREPARING EDR OUTPUT ########### */
1346 if (fplog && do_log && bDoExpanded)
1348 /* only needed if doing expanded ensemble */
1349 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1350 state_global->dfhist, state->fep_state, ir->nstlog, step);
1354 upd_mdebin(mdebin, bDoDHDL, bCalcEnerStep,
1355 t, mdatoms->tmass, enerd, state,
1356 ir->fepvals, ir->expandedvals, lastbox,
1357 shake_vir, force_vir, total_vir, pres,
1358 ekind, mu_tot, constr);
1362 upd_mdebin_step(mdebin);
1365 gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
1366 gmx_bool do_or = do_per_step(step, ir->nstorireout);
1368 print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
1370 eprNORMAL, mdebin, 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) ? fr->nbv : 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 print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
1505 eprAVER, mdebin, fcd, groups, &(ir->opts), awh.get());
1508 done_mdebin(mdebin);
1513 pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
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);