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/energyoutput.h"
79 #include "gromacs/mdlib/expanded.h"
80 #include "gromacs/mdlib/force.h"
81 #include "gromacs/mdlib/force_flags.h"
82 #include "gromacs/mdlib/forcerec.h"
83 #include "gromacs/mdlib/md_support.h"
84 #include "gromacs/mdlib/mdatoms.h"
85 #include "gromacs/mdlib/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, outputProvider, ir, top_global, oenv, wcycle);
257 gmx::EnergyOutput energyOutput;
258 energyOutput.prepare(mdoutf_get_fp_ene(outf), top_global, ir, mdoutf_get_fp_dhdl(outf));
260 /* Energy terms and groups */
262 init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
265 /* Kinetic energy data */
266 std::unique_ptr<gmx_ekindata_t> eKinData = std::make_unique<gmx_ekindata_t>();
267 gmx_ekindata_t *ekind = eKinData.get();
268 init_ekindata(fplog, top_global, &(ir->opts), ekind);
269 /* Copy the cos acceleration to the groups struct */
270 ekind->cosacc.cos_accel = ir->cos_accel;
272 gstat = global_stat_init(ir);
274 /* Check for polarizable models and flexible constraints */
275 shellfc = init_shell_flexcon(fplog,
276 top_global, constr ? constr->numFlexibleConstraints() : 0,
277 ir->nstcalcenergy, DOMAINDECOMP(cr));
280 double io = compute_io(ir, top_global->natoms, groups, energyOutput.numEnergyTerms(), 1);
281 if ((io > 2000) && MASTER(cr))
284 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
289 /* Set up interactive MD (IMD) */
290 init_IMD(ir, cr, ms, top_global, fplog, ir->nstcalcenergy,
291 MASTER(cr) ? state_global->x.rvec_array() : nullptr,
292 nfile, fnm, oenv, mdrunOptions);
294 // Local state only becomes valid now.
295 std::unique_ptr<t_state> stateInstance;
298 if (DOMAINDECOMP(cr))
300 dd_init_local_top(*top_global, &top);
302 stateInstance = std::make_unique<t_state>();
303 state = stateInstance.get();
304 dd_init_local_state(cr->dd, state_global, state);
306 /* Distribute the charge groups over the nodes from the master node */
307 dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
308 state_global, *top_global, ir,
309 state, &f, mdAtoms, &top, fr,
311 nrnb, nullptr, FALSE);
312 shouldCheckNumberOfBondedInteractions = true;
313 upd.setNumAtoms(state->natoms);
317 state_change_natoms(state_global, state_global->natoms);
318 f.resizeWithPadding(state_global->natoms);
319 /* Copy the pointer to the global state */
320 state = state_global;
322 /* Generate and initialize new topology */
323 mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr,
324 &graph, mdAtoms, constr, vsite, shellfc);
326 upd.setNumAtoms(state->natoms);
329 auto mdatoms = mdAtoms->mdatoms();
331 // NOTE: The global state is no longer used at this point.
332 // But state_global is still used as temporary storage space for writing
333 // the global state to file and potentially for replica exchange.
334 // (Global topology should persist.)
336 update_mdatoms(mdatoms, state->lambda[efptMASS]);
338 const ContinuationOptions &continuationOptions = mdrunOptions.continuationOptions;
339 bool startingFromCheckpoint = continuationOptions.startedFromCheckpoint;
343 /* Check nstexpanded here, because the grompp check was broken */
344 if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
346 gmx_fatal(FARGS, "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
348 init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
353 if (startingFromCheckpoint)
355 /* Restore from energy history if appending to output files */
356 if (continuationOptions.appendFiles)
358 /* If no history is available (because a checkpoint is from before
359 * it was written) make a new one later, otherwise restore it.
361 if (observablesHistory->energyHistory)
363 energyOutput.restoreFromEnergyHistory(*observablesHistory->energyHistory);
366 else if (observablesHistory->energyHistory)
368 /* We might have read an energy history from checkpoint.
369 * As we are not appending, we want to restart the statistics.
370 * Free the allocated memory and reset the counts.
372 observablesHistory->energyHistory = {};
373 /* We might have read a pull history from checkpoint.
374 * We will still want to keep the statistics, so that the files
375 * can be joined and still be meaningful.
376 * This means that observablesHistory->pullHistory
377 * should not be reset.
381 if (!observablesHistory->energyHistory)
383 observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
385 if (!observablesHistory->pullHistory)
387 observablesHistory->pullHistory = std::make_unique<PullHistory>();
389 /* Set the initial energy history */
390 energyOutput.fillEnergyHistory(observablesHistory->energyHistory.get());
393 preparePrevStepPullCom(ir, mdatoms, state, state_global, cr, startingFromCheckpoint);
395 // TODO: Remove this by converting AWH into a ForceProvider
396 auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingFromCheckpoint,
398 opt2fn("-awh", nfile, fnm), ir->pull_work);
400 const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
401 if (useReplicaExchange && MASTER(cr))
403 repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
406 /* PME tuning is only supported in the Verlet scheme, with PME for
407 * Coulomb. It is not supported with only LJ PME. */
408 bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) &&
409 !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
411 pme_load_balancing_t *pme_loadbal = nullptr;
414 pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box,
415 *fr->ic, *fr->nbv->listParams, fr->pmedata, use_GPU(fr->nbv),
419 if (!ir->bContinuation)
421 if (state->flags & (1 << estV))
423 auto v = makeArrayRef(state->v);
424 /* Set the velocities of vsites, shells and frozen atoms to zero */
425 for (i = 0; i < mdatoms->homenr; i++)
427 if (mdatoms->ptype[i] == eptVSite ||
428 mdatoms->ptype[i] == eptShell)
432 else if (mdatoms->cFREEZE)
434 for (m = 0; m < DIM; m++)
436 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
447 /* Constrain the initial coordinates and velocities */
448 do_constrain_first(fplog, constr, ir, mdatoms, state);
452 /* Construct the virtual sites for the initial configuration */
453 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
454 top.idef.iparams, top.idef.il,
455 fr->ePBC, fr->bMolPBC, cr, state->box);
459 if (ir->efep != efepNO)
461 /* Set free energy calculation frequency as the greatest common
462 * denominator of nstdhdl and repl_ex_nst. */
463 nstfep = ir->fepvals->nstdhdl;
466 nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
468 if (useReplicaExchange)
470 nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
474 /* Be REALLY careful about what flags you set here. You CANNOT assume
475 * this is the first step, since we might be restarting from a checkpoint,
476 * and in that case we should not do any modifications to the state.
478 bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
480 if (continuationOptions.haveReadEkin)
482 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
485 cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT
486 | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
487 | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
488 | (continuationOptions.haveReadEkin ? CGLO_READEKIN : 0));
490 bSumEkinhOld = FALSE;
492 t_vcm vcm(top_global->groups, *ir);
493 reportComRemovalInfo(fplog, vcm);
495 /* To minimize communication, compute_globals computes the COM velocity
496 * and the kinetic energy for the velocities without COM motion removed.
497 * Thus to get the kinetic energy without the COM contribution, we need
498 * to call compute_globals twice.
500 for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
502 int cglo_flags_iteration = cglo_flags;
503 if (bStopCM && cgloIteration == 0)
505 cglo_flags_iteration |= CGLO_STOPCM;
506 cglo_flags_iteration &= ~CGLO_TEMPERATURE;
508 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
509 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
510 constr, &nullSignaller, state->box,
511 &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
512 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
514 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
515 top_global, &top, state,
516 &shouldCheckNumberOfBondedInteractions);
517 if (ir->eI == eiVVAK)
519 /* a second call to get the half step temperature initialized as well */
520 /* we do the same call as above, but turn the pressure off -- internally to
521 compute_globals, this is recognized as a velocity verlet half-step
522 kinetic energy calculation. This minimized excess variables, but
523 perhaps loses some logic?*/
525 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
526 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
527 constr, &nullSignaller, state->box,
528 nullptr, &bSumEkinhOld,
529 cglo_flags & ~CGLO_PRESSURE);
532 /* Calculate the initial half step temperature, and save the ekinh_old */
533 if (!continuationOptions.startedFromCheckpoint)
535 for (i = 0; (i < ir->opts.ngtc); i++)
537 copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
541 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
542 temperature control */
543 auto trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
547 if (!ir->bContinuation)
549 if (constr && ir->eConstrAlg == econtLINCS)
552 "RMS relative constraint deviation after constraining: %.2e\n",
555 if (EI_STATE_VELOCITY(ir->eI))
557 real temp = enerd->term[F_TEMP];
560 /* Result of Ekin averaged over velocities of -half
561 * and +half step, while we only have -half step here.
565 fprintf(fplog, "Initial temperature: %g K\n", temp);
570 fprintf(stderr, "starting mdrun '%s'\n",
571 *(top_global->name));
574 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
578 sprintf(tbuf, "%s", "infinite");
580 if (ir->init_step > 0)
582 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
583 gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
584 gmx_step_str(ir->init_step, sbuf2),
585 ir->init_step*ir->delta_t);
589 fprintf(stderr, "%s steps, %s ps.\n",
590 gmx_step_str(ir->nsteps, sbuf), tbuf);
592 fprintf(fplog, "\n");
595 walltime_accounting_start_time(walltime_accounting);
596 wallcycle_start(wcycle, ewcRUN);
597 print_start(fplog, cr, walltime_accounting, "mdrun");
600 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
601 int chkpt_ret = fcCheckPointParallel( cr->nodeid,
605 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
609 /***********************************************************
613 ************************************************************/
616 /* Skip the first Nose-Hoover integration when we get the state from tpx */
617 bInitStep = !startingFromCheckpoint || EI_VV(ir->eI);
618 bSumEkinhOld = FALSE;
620 bNeedRepartition = FALSE;
622 bool simulationsShareState = false;
623 int nstSignalComm = nstglobalcomm;
625 // TODO This implementation of ensemble orientation restraints is nasty because
626 // a user can't just do multi-sim with single-sim orientation restraints.
627 bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
628 bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
630 // Replica exchange, ensemble restraints and AWH need all
631 // simulations to remain synchronized, so they need
632 // checkpoints and stop conditions to act on the same step, so
633 // the propagation of such signals must take place between
634 // simulations, not just within simulations.
635 // TODO: Make algorithm initializers set these flags.
636 simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
638 if (simulationsShareState)
640 // Inter-simulation signal communication does not need to happen
641 // often, so we use a minimum of 200 steps to reduce overhead.
642 const int c_minimumInterSimulationSignallingInterval = 200;
643 nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
647 auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
648 compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
649 MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
650 mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
652 auto checkpointHandler = std::make_unique<CheckpointHandler>(
653 compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]),
654 simulationsShareState, ir->nstlist == 0, MASTER(cr),
655 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
657 const bool resetCountersIsLocal = true;
658 auto resetHandler = std::make_unique<ResetHandler>(
659 compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]), !resetCountersIsLocal,
660 ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
661 mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
663 DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
664 DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion = (DOMAINDECOMP(cr) ? DdCloseBalanceRegionAfterForceComputation::yes : DdCloseBalanceRegionAfterForceComputation::no);
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 ddOpenBalanceRegion, ddCloseBalanceRegion);
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 ddOpenBalanceRegion, ddCloseBalanceRegion);
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) ? 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 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));
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);