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/compat/make_unique.h"
57 #include "gromacs/domdec/collect.h"
58 #include "gromacs/domdec/domdec.h"
59 #include "gromacs/domdec/domdec_network.h"
60 #include "gromacs/domdec/domdec_struct.h"
61 #include "gromacs/domdec/partition.h"
62 #include "gromacs/essentialdynamics/edsam.h"
63 #include "gromacs/ewald/pme.h"
64 #include "gromacs/ewald/pme-load-balancing.h"
65 #include "gromacs/fileio/trxio.h"
66 #include "gromacs/gmxlib/network.h"
67 #include "gromacs/gmxlib/nrnb.h"
68 #include "gromacs/gpu_utils/gpu_utils.h"
69 #include "gromacs/imd/imd.h"
70 #include "gromacs/listed_forces/manage-threading.h"
71 #include "gromacs/math/functions.h"
72 #include "gromacs/math/utilities.h"
73 #include "gromacs/math/vec.h"
74 #include "gromacs/math/vectypes.h"
75 #include "gromacs/mdlib/checkpointhandler.h"
76 #include "gromacs/mdlib/compute_io.h"
77 #include "gromacs/mdlib/constr.h"
78 #include "gromacs/mdlib/ebin.h"
79 #include "gromacs/mdlib/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/mdebin.h"
86 #include "gromacs/mdlib/mdoutf.h"
87 #include "gromacs/mdlib/mdrun.h"
88 #include "gromacs/mdlib/mdsetup.h"
89 #include "gromacs/mdlib/membed.h"
90 #include "gromacs/mdlib/nb_verlet.h"
91 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
92 #include "gromacs/mdlib/ns.h"
93 #include "gromacs/mdlib/resethandler.h"
94 #include "gromacs/mdlib/shellfc.h"
95 #include "gromacs/mdlib/sighandler.h"
96 #include "gromacs/mdlib/sim_util.h"
97 #include "gromacs/mdlib/simulationsignal.h"
98 #include "gromacs/mdlib/stophandler.h"
99 #include "gromacs/mdlib/tgroup.h"
100 #include "gromacs/mdlib/trajectory_writing.h"
101 #include "gromacs/mdlib/update.h"
102 #include "gromacs/mdlib/vcm.h"
103 #include "gromacs/mdlib/vsite.h"
104 #include "gromacs/mdtypes/awh-history.h"
105 #include "gromacs/mdtypes/awh-params.h"
106 #include "gromacs/mdtypes/commrec.h"
107 #include "gromacs/mdtypes/df_history.h"
108 #include "gromacs/mdtypes/energyhistory.h"
109 #include "gromacs/mdtypes/fcdata.h"
110 #include "gromacs/mdtypes/forcerec.h"
111 #include "gromacs/mdtypes/group.h"
112 #include "gromacs/mdtypes/inputrec.h"
113 #include "gromacs/mdtypes/interaction_const.h"
114 #include "gromacs/mdtypes/md_enums.h"
115 #include "gromacs/mdtypes/mdatom.h"
116 #include "gromacs/mdtypes/observableshistory.h"
117 #include "gromacs/mdtypes/pullhistory.h"
118 #include "gromacs/mdtypes/state.h"
119 #include "gromacs/pbcutil/mshift.h"
120 #include "gromacs/pbcutil/pbc.h"
121 #include "gromacs/pulling/output.h"
122 #include "gromacs/pulling/pull.h"
123 #include "gromacs/swap/swapcoords.h"
124 #include "gromacs/timing/wallcycle.h"
125 #include "gromacs/timing/walltime_accounting.h"
126 #include "gromacs/topology/atoms.h"
127 #include "gromacs/topology/idef.h"
128 #include "gromacs/topology/mtop_util.h"
129 #include "gromacs/topology/topology.h"
130 #include "gromacs/trajectory/trajectoryframe.h"
131 #include "gromacs/utility/basedefinitions.h"
132 #include "gromacs/utility/cstringutil.h"
133 #include "gromacs/utility/fatalerror.h"
134 #include "gromacs/utility/logger.h"
135 #include "gromacs/utility/real.h"
136 #include "gromacs/utility/smalloc.h"
138 #include "integrator.h"
139 #include "replicaexchange.h"
142 #include "corewrap.h"
145 using gmx::SimulationSignaller;
147 void gmx::Integrator::do_md()
149 // TODO Historically, the EM and MD "integrators" used different
150 // names for the t_inputrec *parameter, but these must have the
151 // same name, now that it's a member of a struct. We use this ir
152 // alias to avoid a large ripple of nearly useless changes.
153 // t_inputrec is being replaced by IMdpOptionsProvider, so this
154 // will go away eventually.
155 t_inputrec *ir = inputrec;
156 int64_t step, step_rel;
157 double t = ir->init_t, t0 = ir->init_t, lam0[efptNR];
158 gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
159 gmx_bool bNS, bNStList, bStopCM,
160 bFirstStep, bInitStep, bLastStep = FALSE;
161 gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
162 gmx_bool do_ene, do_log, do_verbose;
163 gmx_bool bMasterState;
164 int force_flags, cglo_flags;
165 tensor force_vir = {{0}}, shake_vir = {{0}}, total_vir = {{0}},
166 tmp_vir = {{0}}, pres = {{0}};
169 matrix parrinellorahmanMu, M;
170 gmx_repl_ex_t repl_ex = nullptr;
172 gmx_enerdata_t *enerd;
173 PaddedVector<gmx::RVec> f {};
174 gmx_global_stat_t gstat;
175 t_graph *graph = nullptr;
176 gmx_groups_t *groups;
177 gmx_shellfc_t *shellfc;
178 gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
179 gmx_bool bTemp, bPres, bTrotter;
181 rvec *cbuf = nullptr;
188 real saved_conserved_quantity = 0;
191 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
193 /* PME load balancing data for GPU kernels */
194 gmx_bool bPMETune = FALSE;
195 gmx_bool bPMETunePrinting = FALSE;
198 gmx_bool bIMDstep = FALSE;
200 /* Domain decomposition could incorrectly miss a bonded
201 interaction, but checking for that requires a global
202 communication stage, which does not otherwise happen in DD
203 code. So we do that alongside the first global energy reduction
204 after a new DD is made. These variables handle whether the
205 check happens, and the result it returns. */
206 bool shouldCheckNumberOfBondedInteractions = false;
207 int totalNumberOfBondedInteractions = -1;
209 SimulationSignals signals;
210 // Most global communnication stages don't propagate mdrun
211 // signals, and will use this object to achieve that.
212 SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
214 if (!mdrunOptions.writeConfout)
216 // This is on by default, and the main known use case for
217 // turning it off is for convenience in benchmarking, which is
218 // something that should not show up in the general user
220 GMX_LOG(mdlog.info).asParagraph().
221 appendText("The -noconfout functionality is deprecated, and may be removed in a future version.");
224 /* md-vv uses averaged full step velocities for T-control
225 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
226 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
227 bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
229 const bool bRerunMD = false;
230 int nstglobalcomm = mdrunOptions.globalCommunicationInterval;
232 nstglobalcomm = check_nstglobalcomm(mdlog, nstglobalcomm, ir, cr);
233 bGStatEveryStep = (nstglobalcomm == 1);
235 groups = &top_global->groups;
237 std::unique_ptr<EssentialDynamics> ed = nullptr;
238 if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
240 /* Initialize essential dynamics sampling */
241 ed = init_edsam(mdlog,
242 opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
245 state_global, observablesHistory,
246 oenv, mdrunOptions.continuationOptions.appendFiles);
249 initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda, lam0);
250 Update upd(ir, deform);
251 bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
252 if (!mdrunOptions.continuationOptions.appendFiles)
254 pleaseCiteCouplingAlgorithms(fplog, *ir);
257 gmx_mdoutf *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr,
258 outputProvider, ir, top_global, oenv, wcycle);
259 t_mdebin *mdebin = init_mdebin(mdrunOptions.continuationOptions.appendFiles ? nullptr : mdoutf_get_fp_ene(outf),
260 top_global, ir, mdoutf_get_fp_dhdl(outf));
262 /* Energy terms and groups */
264 init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
267 /* Kinetic energy data */
268 std::unique_ptr<gmx_ekindata_t> eKinData = compat::make_unique<gmx_ekindata_t>();
269 gmx_ekindata_t *ekind = eKinData.get();
270 init_ekindata(fplog, top_global, &(ir->opts), ekind);
271 /* Copy the cos acceleration to the groups struct */
272 ekind->cosacc.cos_accel = ir->cos_accel;
274 gstat = global_stat_init(ir);
276 /* Check for polarizable models and flexible constraints */
277 shellfc = init_shell_flexcon(fplog,
278 top_global, constr ? constr->numFlexibleConstraints() : 0,
279 ir->nstcalcenergy, DOMAINDECOMP(cr));
282 double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1);
283 if ((io > 2000) && MASTER(cr))
286 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
291 /* Set up interactive MD (IMD) */
292 init_IMD(ir, cr, ms, top_global, fplog, ir->nstcalcenergy,
293 MASTER(cr) ? state_global->x.rvec_array() : nullptr,
294 nfile, fnm, oenv, mdrunOptions);
296 // Local state only becomes valid now.
297 std::unique_ptr<t_state> stateInstance;
300 if (DOMAINDECOMP(cr))
302 dd_init_local_top(*top_global, &top);
304 stateInstance = compat::make_unique<t_state>();
305 state = stateInstance.get();
306 dd_init_local_state(cr->dd, state_global, state);
308 /* Distribute the charge groups over the nodes from the master node */
309 dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
310 state_global, *top_global, ir,
311 state, &f, mdAtoms, &top, fr,
313 nrnb, nullptr, FALSE);
314 shouldCheckNumberOfBondedInteractions = true;
315 upd.setNumAtoms(state->natoms);
319 state_change_natoms(state_global, state_global->natoms);
320 f.resizeWithPadding(state_global->natoms);
321 /* Copy the pointer to the global state */
322 state = state_global;
324 /* Generate and initialize new topology */
325 mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr,
326 &graph, mdAtoms, constr, vsite, shellfc);
328 upd.setNumAtoms(state->natoms);
331 auto mdatoms = mdAtoms->mdatoms();
333 // NOTE: The global state is no longer used at this point.
334 // But state_global is still used as temporary storage space for writing
335 // the global state to file and potentially for replica exchange.
336 // (Global topology should persist.)
338 update_mdatoms(mdatoms, state->lambda[efptMASS]);
340 const ContinuationOptions &continuationOptions = mdrunOptions.continuationOptions;
341 bool startingFromCheckpoint = continuationOptions.startedFromCheckpoint;
345 /* Check nstexpanded here, because the grompp check was broken */
346 if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
348 gmx_fatal(FARGS, "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
350 init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
355 if (startingFromCheckpoint)
357 /* Update mdebin with energy history if appending to output files */
358 if (continuationOptions.appendFiles)
360 /* If no history is available (because a checkpoint is from before
361 * it was written) make a new one later, otherwise restore it.
363 if (observablesHistory->energyHistory)
365 restore_energyhistory_from_state(mdebin, observablesHistory->energyHistory.get());
368 else if (observablesHistory->energyHistory)
370 /* We might have read an energy history from checkpoint.
371 * As we are not appending, we want to restart the statistics.
372 * Free the allocated memory and reset the counts.
374 observablesHistory->energyHistory = {};
375 /* We might have read a pull history from checkpoint.
376 * We will still want to keep the statistics, so that the files
377 * can be joined and still be meaningful.
378 * This means that observablesHistory->pullHistory
379 * should not be reset.
383 if (!observablesHistory->energyHistory)
385 observablesHistory->energyHistory = compat::make_unique<energyhistory_t>();
387 if (!observablesHistory->pullHistory)
389 observablesHistory->pullHistory = compat::make_unique<PullHistory>();
391 /* Set the initial energy history in state by updating once */
392 update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
395 preparePrevStepPullCom(ir, mdatoms, state, state_global, cr, startingFromCheckpoint);
397 // TODO: Remove this by converting AWH into a ForceProvider
398 auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingFromCheckpoint,
400 opt2fn("-awh", nfile, fnm), ir->pull_work);
402 const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
403 if (useReplicaExchange && MASTER(cr))
405 repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
408 /* PME tuning is only supported in the Verlet scheme, with PME for
409 * Coulomb. It is not supported with only LJ PME. */
410 bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) &&
411 !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
413 pme_load_balancing_t *pme_loadbal = nullptr;
416 pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box,
417 *fr->ic, *fr->nbv->listParams, fr->pmedata, use_GPU(fr->nbv),
421 if (!ir->bContinuation)
423 if (state->flags & (1 << estV))
425 auto v = makeArrayRef(state->v);
426 /* Set the velocities of vsites, shells and frozen atoms to zero */
427 for (i = 0; i < mdatoms->homenr; i++)
429 if (mdatoms->ptype[i] == eptVSite ||
430 mdatoms->ptype[i] == eptShell)
434 else if (mdatoms->cFREEZE)
436 for (m = 0; m < DIM; m++)
438 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
449 /* Constrain the initial coordinates and velocities */
450 do_constrain_first(fplog, constr, ir, mdatoms, state);
454 /* Construct the virtual sites for the initial configuration */
455 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
456 top.idef.iparams, top.idef.il,
457 fr->ePBC, fr->bMolPBC, cr, state->box);
461 if (ir->efep != efepNO)
463 /* Set free energy calculation frequency as the greatest common
464 * denominator of nstdhdl and repl_ex_nst. */
465 nstfep = ir->fepvals->nstdhdl;
468 nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
470 if (useReplicaExchange)
472 nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
476 /* Be REALLY careful about what flags you set here. You CANNOT assume
477 * this is the first step, since we might be restarting from a checkpoint,
478 * and in that case we should not do any modifications to the state.
480 bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
482 if (continuationOptions.haveReadEkin)
484 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
487 cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT
488 | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
489 | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
490 | (continuationOptions.haveReadEkin ? CGLO_READEKIN : 0));
492 bSumEkinhOld = FALSE;
494 t_vcm vcm(top_global->groups, *ir);
495 reportComRemovalInfo(fplog, vcm);
497 /* To minimize communication, compute_globals computes the COM velocity
498 * and the kinetic energy for the velocities without COM motion removed.
499 * Thus to get the kinetic energy without the COM contribution, we need
500 * to call compute_globals twice.
502 for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
504 int cglo_flags_iteration = cglo_flags;
505 if (bStopCM && cgloIteration == 0)
507 cglo_flags_iteration |= CGLO_STOPCM;
508 cglo_flags_iteration &= ~CGLO_TEMPERATURE;
510 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
511 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
512 constr, &nullSignaller, state->box,
513 &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
514 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
516 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
517 top_global, &top, state,
518 &shouldCheckNumberOfBondedInteractions);
519 if (ir->eI == eiVVAK)
521 /* a second call to get the half step temperature initialized as well */
522 /* we do the same call as above, but turn the pressure off -- internally to
523 compute_globals, this is recognized as a velocity verlet half-step
524 kinetic energy calculation. This minimized excess variables, but
525 perhaps loses some logic?*/
527 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
528 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
529 constr, &nullSignaller, state->box,
530 nullptr, &bSumEkinhOld,
531 cglo_flags & ~CGLO_PRESSURE);
534 /* Calculate the initial half step temperature, and save the ekinh_old */
535 if (!continuationOptions.startedFromCheckpoint)
537 for (i = 0; (i < ir->opts.ngtc); i++)
539 copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
543 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
544 temperature control */
545 auto trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
549 if (!ir->bContinuation)
551 if (constr && ir->eConstrAlg == econtLINCS)
554 "RMS relative constraint deviation after constraining: %.2e\n",
557 if (EI_STATE_VELOCITY(ir->eI))
559 real temp = enerd->term[F_TEMP];
562 /* Result of Ekin averaged over velocities of -half
563 * and +half step, while we only have -half step here.
567 fprintf(fplog, "Initial temperature: %g K\n", temp);
572 fprintf(stderr, "starting mdrun '%s'\n",
573 *(top_global->name));
576 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
580 sprintf(tbuf, "%s", "infinite");
582 if (ir->init_step > 0)
584 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
585 gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
586 gmx_step_str(ir->init_step, sbuf2),
587 ir->init_step*ir->delta_t);
591 fprintf(stderr, "%s steps, %s ps.\n",
592 gmx_step_str(ir->nsteps, sbuf), tbuf);
594 fprintf(fplog, "\n");
597 walltime_accounting_start_time(walltime_accounting);
598 wallcycle_start(wcycle, ewcRUN);
599 print_start(fplog, cr, walltime_accounting, "mdrun");
602 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
603 int chkpt_ret = fcCheckPointParallel( cr->nodeid,
607 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
611 /***********************************************************
615 ************************************************************/
618 /* Skip the first Nose-Hoover integration when we get the state from tpx */
619 bInitStep = !startingFromCheckpoint || EI_VV(ir->eI);
620 bSumEkinhOld = FALSE;
622 bNeedRepartition = FALSE;
624 bool simulationsShareState = false;
625 int nstSignalComm = nstglobalcomm;
627 // TODO This implementation of ensemble orientation restraints is nasty because
628 // a user can't just do multi-sim with single-sim orientation restraints.
629 bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
630 bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
632 // Replica exchange, ensemble restraints and AWH need all
633 // simulations to remain synchronized, so they need
634 // checkpoints and stop conditions to act on the same step, so
635 // the propagation of such signals must take place between
636 // simulations, not just within simulations.
637 // TODO: Make algorithm initializers set these flags.
638 simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
640 if (simulationsShareState)
642 // Inter-simulation signal communication does not need to happen
643 // often, so we use a minimum of 200 steps to reduce overhead.
644 const int c_minimumInterSimulationSignallingInterval = 200;
645 nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
649 auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
650 compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
651 MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
652 mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
654 auto checkpointHandler = compat::make_unique<CheckpointHandler>(
655 compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]),
656 simulationsShareState, ir->nstlist == 0, MASTER(cr),
657 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
659 const bool resetCountersIsLocal = true;
660 auto resetHandler = compat::make_unique<ResetHandler>(
661 compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]), !resetCountersIsLocal,
662 ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
663 mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
665 DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
666 DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion = (DOMAINDECOMP(cr) ? DdCloseBalanceRegionAfterForceComputation::yes : DdCloseBalanceRegionAfterForceComputation::no);
668 step = ir->init_step;
671 // TODO extract this to new multi-simulation module
672 if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
674 if (!multisim_int_all_are_equal(ms, ir->nsteps))
676 GMX_LOG(mdlog.warning).appendText(
677 "Note: The number of steps is not consistent across multi simulations,\n"
678 "but we are proceeding anyway!");
680 if (!multisim_int_all_are_equal(ms, ir->init_step))
682 GMX_LOG(mdlog.warning).appendText(
683 "Note: The initial step is not consistent across multi simulations,\n"
684 "but we are proceeding anyway!");
688 /* and stop now if we should */
689 bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
693 /* Determine if this is a neighbor search step */
694 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
696 if (bPMETune && bNStList)
698 /* PME grid + cut-off optimization with GPUs or PME nodes */
699 pme_loadbal_do(pme_loadbal, cr,
700 (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
708 wallcycle_start(wcycle, ewcSTEP);
710 bLastStep = (step_rel == ir->nsteps);
711 t = t0 + step*ir->delta_t;
713 // TODO Refactor this, so that nstfep does not need a default value of zero
714 if (ir->efep != efepNO || ir->bSimTemp)
716 /* find and set the current lambdas */
717 setCurrentLambdasLocal(step, ir->fepvals, lam0, state);
719 bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
720 bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
721 bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
722 && (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
725 bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
726 do_per_step(step, replExParams.exchangeInterval));
728 if (doSimulatedAnnealing)
730 update_annealing_target_temp(ir, t, &upd);
733 /* Stop Center of Mass motion */
734 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
736 /* Determine whether or not to do Neighbour Searching */
737 bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
739 bLastStep = bLastStep || stopHandler->stoppingAfterCurrentStep(bNS);
741 /* do_log triggers energy and virial calculation. Because this leads
742 * to different code paths, forces can be different. Thus for exact
743 * continuation we should avoid extra log output.
744 * Note that the || bLastStep can result in non-exact continuation
745 * beyond the last step. But we don't consider that to be an issue.
747 do_log = do_per_step(step, ir->nstlog) || (bFirstStep && !startingFromCheckpoint) || bLastStep;
748 do_verbose = mdrunOptions.verbose &&
749 (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
751 if (bNS && !(bFirstStep && ir->bContinuation))
753 bMasterState = FALSE;
754 /* Correct the new box if it is too skewed */
755 if (inputrecDynamicBox(ir))
757 if (correct_box(fplog, step, state->box, graph))
762 if (DOMAINDECOMP(cr) && bMasterState)
764 dd_collect_state(cr->dd, state, state_global);
767 if (DOMAINDECOMP(cr))
769 /* Repartition the domain decomposition */
770 dd_partition_system(fplog, mdlog, step, cr,
771 bMasterState, nstglobalcomm,
772 state_global, *top_global, ir,
773 state, &f, mdAtoms, &top, fr,
776 do_verbose && !bPMETunePrinting);
777 shouldCheckNumberOfBondedInteractions = true;
778 upd.setNumAtoms(state->natoms);
782 if (MASTER(cr) && do_log)
784 print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
787 if (ir->efep != efepNO)
789 update_mdatoms(mdatoms, state->lambda[efptMASS]);
795 /* We need the kinetic energy at minus the half step for determining
796 * the full step kinetic energy and possibly for T-coupling.*/
797 /* This may not be quite working correctly yet . . . . */
798 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
799 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
800 constr, &nullSignaller, state->box,
801 &totalNumberOfBondedInteractions, &bSumEkinhOld,
802 CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
803 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
804 top_global, &top, state,
805 &shouldCheckNumberOfBondedInteractions);
807 clear_mat(force_vir);
809 checkpointHandler->decideIfCheckpointingThisStep(bNS, bFirstStep, bLastStep);
811 /* Determine the energy and pressure:
812 * at nstcalcenergy steps and at energy output steps (set below).
814 if (EI_VV(ir->eI) && (!bInitStep))
816 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
817 bCalcVir = bCalcEnerStep ||
818 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
822 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
823 bCalcVir = bCalcEnerStep ||
824 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
826 bCalcEner = bCalcEnerStep;
828 do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
830 if (do_ene || do_log || bDoReplEx)
836 /* Do we need global communication ? */
837 bGStat = (bCalcVir || bCalcEner || bStopCM ||
838 do_per_step(step, nstglobalcomm) ||
839 (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
841 force_flags = (GMX_FORCE_STATECHANGED |
842 ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0) |
843 GMX_FORCE_ALLFORCES |
844 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
845 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
846 (bDoFEP ? GMX_FORCE_DHDL : 0)
851 /* Now is the time to relax the shells */
852 relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
853 enforcedRotation, step,
854 ir, bNS, force_flags, &top,
856 state, f.arrayRefWithPadding(), force_vir, mdatoms,
857 nrnb, wcycle, graph, groups,
858 shellfc, fr, ppForceWorkload, t, mu_tot,
860 ddOpenBalanceRegion, ddCloseBalanceRegion);
864 /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
865 (or the AWH update will be performed twice for one step when continuing). It would be best to
866 call this update function from do_md_trajectory_writing but that would occur after do_force.
867 One would have to divide the update_awh function into one function applying the AWH force
868 and one doing the AWH bias update. The update AWH bias function could then be called after
869 do_md_trajectory_writing (then containing update_awh_history).
870 The checkpointing will in the future probably moved to the start of the md loop which will
871 rid of this issue. */
872 if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
874 awh->updateHistory(state_global->awhHistory.get());
877 /* The coordinates (x) are shifted (to get whole molecules)
879 * This is parallellized as well, and does communication too.
880 * Check comments in sim_util.c
882 do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation,
883 step, nrnb, wcycle, &top, groups,
884 state->box, state->x.arrayRefWithPadding(), &state->hist,
885 f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
886 state->lambda, graph,
887 fr, ppForceWorkload, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
888 (bNS ? GMX_FORCE_NS : 0) | force_flags,
889 ddOpenBalanceRegion, ddCloseBalanceRegion);
892 if (EI_VV(ir->eI) && !startingFromCheckpoint)
893 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
895 rvec *vbuf = nullptr;
897 wallcycle_start(wcycle, ewcUPDATE);
898 if (ir->eI == eiVV && bInitStep)
900 /* if using velocity verlet with full time step Ekin,
901 * take the first half step only to compute the
902 * virial for the first step. From there,
903 * revert back to the initial coordinates
904 * so that the input is actually the initial step.
906 snew(vbuf, state->natoms);
907 copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
911 /* this is for NHC in the Ekin(t+dt/2) version of vv */
912 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
915 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
916 ekind, M, &upd, etrtVELOCITY1,
919 wallcycle_stop(wcycle, ewcUPDATE);
920 constrain_velocities(step, nullptr,
924 bCalcVir, do_log, do_ene);
925 wallcycle_start(wcycle, ewcUPDATE);
926 /* if VV, compute the pressure and constraints */
927 /* For VV2, we strictly only need this if using pressure
928 * control, but we really would like to have accurate pressures
930 * Think about ways around this in the future?
931 * For now, keep this choice in comments.
933 /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
934 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
936 bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
937 if (bCalcEner && ir->eI == eiVVAK)
941 /* for vv, the first half of the integration actually corresponds to the previous step.
942 So we need information from the last step in the first half of the integration */
943 if (bGStat || do_per_step(step-1, nstglobalcomm))
945 wallcycle_stop(wcycle, ewcUPDATE);
946 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
947 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
948 constr, &nullSignaller, state->box,
949 &totalNumberOfBondedInteractions, &bSumEkinhOld,
950 (bGStat ? CGLO_GSTAT : 0)
951 | (bCalcEner ? CGLO_ENERGY : 0)
952 | (bTemp ? CGLO_TEMPERATURE : 0)
953 | (bPres ? CGLO_PRESSURE : 0)
954 | (bPres ? CGLO_CONSTRAINT : 0)
955 | (bStopCM ? CGLO_STOPCM : 0)
956 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
959 /* explanation of above:
960 a) We compute Ekin at the full time step
961 if 1) we are using the AveVel Ekin, and it's not the
962 initial step, or 2) if we are using AveEkin, but need the full
963 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
964 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
965 EkinAveVel because it's needed for the pressure */
966 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
967 top_global, &top, state,
968 &shouldCheckNumberOfBondedInteractions);
969 wallcycle_start(wcycle, ewcUPDATE);
971 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
976 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
977 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
979 /* TODO This is only needed when we're about to write
980 * a checkpoint, because we use it after the restart
981 * (in a kludge?). But what should we be doing if
982 * startingFromCheckpoint or bInitStep are true? */
983 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
985 copy_mat(shake_vir, state->svir_prev);
986 copy_mat(force_vir, state->fvir_prev);
988 if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
990 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
991 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
992 enerd->term[F_EKIN] = trace(ekind->ekin);
997 wallcycle_stop(wcycle, ewcUPDATE);
998 /* We need the kinetic energy at minus the half step for determining
999 * the full step kinetic energy and possibly for T-coupling.*/
1000 /* This may not be quite working correctly yet . . . . */
1001 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1002 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1003 constr, &nullSignaller, state->box,
1004 nullptr, &bSumEkinhOld,
1005 CGLO_GSTAT | CGLO_TEMPERATURE);
1006 wallcycle_start(wcycle, ewcUPDATE);
1009 /* if it's the initial step, we performed this first step just to get the constraint virial */
1010 if (ir->eI == eiVV && bInitStep)
1012 copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
1015 wallcycle_stop(wcycle, ewcUPDATE);
1018 /* compute the conserved quantity */
1021 saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1024 last_ekin = enerd->term[F_EKIN];
1026 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1028 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1030 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1031 if (ir->efep != efepNO)
1033 sum_dhdl(enerd, state->lambda, ir->fepvals);
1037 /* ######## END FIRST UPDATE STEP ############## */
1038 /* ######## If doing VV, we now have v(dt) ###### */
1041 /* perform extended ensemble sampling in lambda - we don't
1042 actually move to the new state before outputting
1043 statistics, but if performing simulated tempering, we
1044 do update the velocities and the tau_t. */
1046 lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v.rvec_array(), mdatoms);
1047 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1050 copy_df_history(state_global->dfhist, state->dfhist);
1054 /* Now we have the energies and forces corresponding to the
1055 * coordinates at time t. We must output all of this before
1058 do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1059 ir, state, state_global, observablesHistory,
1061 outf, mdebin, ekind, f,
1062 checkpointHandler->isCheckpointingStep(),
1063 bRerunMD, bLastStep,
1064 mdrunOptions.writeConfout,
1066 /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1067 bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x.rvec_array(), ir, t, wcycle);
1069 /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1070 if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1072 copy_mat(state->svir_prev, shake_vir);
1073 copy_mat(state->fvir_prev, force_vir);
1076 stopHandler->setSignal();
1077 resetHandler->setSignal(walltime_accounting);
1079 if (bGStat || !PAR(cr))
1081 /* In parallel we only have to check for checkpointing in steps
1082 * where we do global communication,
1083 * otherwise the other nodes don't know.
1085 checkpointHandler->setSignal(walltime_accounting);
1088 /* ######### START SECOND UPDATE STEP ################# */
1090 /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1093 if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1095 gmx_bool bIfRandomize;
1096 bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
1097 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1098 if (constr && bIfRandomize)
1100 constrain_velocities(step, nullptr,
1104 bCalcVir, do_log, do_ene);
1107 /* Box is changed in update() when we do pressure coupling,
1108 * but we should still use the old box for energy corrections and when
1109 * writing it to the energy file, so it matches the trajectory files for
1110 * the same timestep above. Make a copy in a separate array.
1112 copy_mat(state->box, lastbox);
1116 wallcycle_start(wcycle, ewcUPDATE);
1117 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1120 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1121 /* We can only do Berendsen coupling after we have summed
1122 * the kinetic energy or virial. Since the happens
1123 * in global_state after update, we should only do it at
1124 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1129 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1130 update_pcouple_before_coordinates(fplog, step, ir, state,
1131 parrinellorahmanMu, M,
1137 /* velocity half-step update */
1138 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1139 ekind, M, &upd, etrtVELOCITY2,
1143 /* Above, initialize just copies ekinh into ekin,
1144 * it doesn't copy position (for VV),
1145 * and entire integrator for MD.
1148 if (ir->eI == eiVVAK)
1150 /* We probably only need md->homenr, not state->natoms */
1151 if (state->natoms > cbuf_nalloc)
1153 cbuf_nalloc = state->natoms;
1154 srenew(cbuf, cbuf_nalloc);
1156 copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1159 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1160 ekind, M, &upd, etrtPOSITION, cr, constr);
1161 wallcycle_stop(wcycle, ewcUPDATE);
1163 constrain_coordinates(step, &dvdl_constr, state,
1166 bCalcVir, do_log, do_ene);
1167 update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
1168 cr, nrnb, wcycle, &upd, constr, do_log, do_ene);
1169 finish_update(ir, mdatoms,
1171 nrnb, wcycle, &upd, constr);
1173 if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1175 updatePrevStepPullCom(ir->pull_work, state);
1178 if (ir->eI == eiVVAK)
1180 /* erase F_EKIN and F_TEMP here? */
1181 /* just compute the kinetic energy at the half step to perform a trotter step */
1182 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1183 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1184 constr, &nullSignaller, lastbox,
1185 nullptr, &bSumEkinhOld,
1186 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1188 wallcycle_start(wcycle, ewcUPDATE);
1189 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1190 /* now we know the scaling, we can compute the positions again again */
1191 copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1193 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1194 ekind, M, &upd, etrtPOSITION, cr, constr);
1195 wallcycle_stop(wcycle, ewcUPDATE);
1197 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1198 /* are the small terms in the shake_vir here due
1199 * to numerical errors, or are they important
1200 * physically? I'm thinking they are just errors, but not completely sure.
1201 * For now, will call without actually constraining, constr=NULL*/
1202 finish_update(ir, mdatoms,
1204 nrnb, wcycle, &upd, nullptr);
1208 /* this factor or 2 correction is necessary
1209 because half of the constraint force is removed
1210 in the vv step, so we have to double it. See
1211 the Redmine issue #1255. It is not yet clear
1212 if the factor of 2 is exact, or just a very
1213 good approximation, and this will be
1214 investigated. The next step is to see if this
1215 can be done adding a dhdl contribution from the
1216 rattle step, but this is somewhat more
1217 complicated with the current code. Will be
1218 investigated, hopefully for 4.6.3. However,
1219 this current solution is much better than
1220 having it completely wrong.
1222 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1226 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1229 if (vsite != nullptr)
1231 wallcycle_start(wcycle, ewcVSITECONSTR);
1232 if (graph != nullptr)
1234 shift_self(graph, state->box, state->x.rvec_array());
1236 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
1237 top.idef.iparams, top.idef.il,
1238 fr->ePBC, fr->bMolPBC, cr, state->box);
1240 if (graph != nullptr)
1242 unshift_self(graph, state->box, state->x.rvec_array());
1244 wallcycle_stop(wcycle, ewcVSITECONSTR);
1247 /* ############## IF NOT VV, Calculate globals HERE ############ */
1248 /* With Leap-Frog we can skip compute_globals at
1249 * non-communication steps, but we need to calculate
1250 * the kinetic energy one step before communication.
1253 // Organize to do inter-simulation signalling on steps if
1254 // and when algorithms require it.
1255 bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1257 if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1259 // Since we're already communicating at this step, we
1260 // can propagate intra-simulation signals. Note that
1261 // check_nstglobalcomm has the responsibility for
1262 // choosing the value of nstglobalcomm that is one way
1263 // bGStat becomes true, so we can't get into a
1264 // situation where e.g. checkpointing can't be
1266 bool doIntraSimSignal = true;
1267 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1269 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1270 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1273 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1274 (bGStat ? CGLO_GSTAT : 0)
1275 | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1276 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1277 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1278 | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
1280 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1282 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1283 top_global, &top, state,
1284 &shouldCheckNumberOfBondedInteractions);
1288 /* ############# END CALC EKIN AND PRESSURE ################# */
1290 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1291 the virial that should probably be addressed eventually. state->veta has better properies,
1292 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1293 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1295 if (ir->efep != efepNO && !EI_VV(ir->eI))
1297 /* Sum up the foreign energy and dhdl terms for md and sd.
1298 Currently done every step so that dhdl is correct in the .edr */
1299 sum_dhdl(enerd, state->lambda, ir->fepvals);
1302 update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1303 pres, force_vir, shake_vir,
1307 /* ################# END UPDATE STEP 2 ################# */
1308 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1310 /* The coordinates (x) were unshifted in update */
1313 /* We will not sum ekinh_old,
1314 * so signal that we still have to do it.
1316 bSumEkinhOld = TRUE;
1321 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1323 /* use the directly determined last velocity, not actually the averaged half steps */
1324 if (bTrotter && ir->eI == eiVV)
1326 enerd->term[F_EKIN] = last_ekin;
1328 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1330 if (integratorHasConservedEnergyQuantity(ir))
1334 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1338 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1341 /* ######### END PREPARING EDR OUTPUT ########### */
1347 if (fplog && do_log && bDoExpanded)
1349 /* only needed if doing expanded ensemble */
1350 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1351 state_global->dfhist, state->fep_state, ir->nstlog, step);
1355 upd_mdebin(mdebin, bDoDHDL, bCalcEnerStep,
1356 t, mdatoms->tmass, enerd, state,
1357 ir->fepvals, ir->expandedvals, lastbox,
1358 shake_vir, force_vir, total_vir, pres,
1359 ekind, mu_tot, constr);
1363 upd_mdebin_step(mdebin);
1366 gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
1367 gmx_bool do_or = do_per_step(step, ir->nstorireout);
1369 print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
1371 eprNORMAL, mdebin, fcd, groups, &(ir->opts), awh.get());
1375 pull_print_output(ir->pull_work, step, t);
1378 if (do_per_step(step, ir->nstlog))
1380 if (fflush(fplog) != 0)
1382 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1388 /* Have to do this part _after_ outputting the logfile and the edr file */
1389 /* Gets written into the state at the beginning of next loop*/
1390 state->fep_state = lamnew;
1392 /* Print the remaining wall clock time for the run */
1393 if (isMasterSimMasterRank(ms, cr) &&
1394 (do_verbose || gmx_got_usr_signal()) &&
1399 fprintf(stderr, "\n");
1401 print_time(stderr, walltime_accounting, step, ir, cr);
1404 /* Ion/water position swapping.
1405 * Not done in last step since trajectory writing happens before this call
1406 * in the MD loop and exchanges would be lost anyway. */
1407 bNeedRepartition = FALSE;
1408 if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1409 do_per_step(step, ir->swap->nstswap))
1411 bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
1412 as_rvec_array(state->x.data()),
1414 MASTER(cr) && mdrunOptions.verbose,
1417 if (bNeedRepartition && DOMAINDECOMP(cr))
1419 dd_collect_state(cr->dd, state, state_global);
1423 /* Replica exchange */
1427 bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
1428 state_global, enerd,
1432 if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1434 dd_partition_system(fplog, mdlog, step, cr, TRUE, 1,
1435 state_global, *top_global, ir,
1436 state, &f, mdAtoms, &top, fr,
1438 nrnb, wcycle, FALSE);
1439 shouldCheckNumberOfBondedInteractions = true;
1440 upd.setNumAtoms(state->natoms);
1445 startingFromCheckpoint = false;
1447 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1448 /* With all integrators, except VV, we need to retain the pressure
1449 * at the current step for coupling at the next step.
1451 if ((state->flags & (1<<estPRES_PREV)) &&
1453 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1455 /* Store the pressure in t_state for pressure coupling
1456 * at the next MD step.
1458 copy_mat(pres, state->pres_prev);
1461 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1463 if ( (membed != nullptr) && (!bLastStep) )
1465 rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1468 cycles = wallcycle_stop(wcycle, ewcSTEP);
1469 if (DOMAINDECOMP(cr) && wcycle)
1471 dd_cycles_add(cr->dd, cycles, ddCyclStep);
1474 /* increase the MD step number */
1478 resetHandler->resetCounters(
1479 step, step_rel, mdlog, fplog, cr, (use_GPU(fr->nbv) ? fr->nbv : nullptr),
1480 nrnb, fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
1482 /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1483 IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle);
1486 /* End of main MD loop */
1488 /* Closing TNG files can include compressing data. Therefore it is good to do that
1489 * before stopping the time measurements. */
1490 mdoutf_tng_close(outf);
1492 /* Stop measuring walltime */
1493 walltime_accounting_end_time(walltime_accounting);
1495 if (!thisRankHasDuty(cr, DUTY_PME))
1497 /* Tell the PME only node to finish */
1498 gmx_pme_send_finish(cr);
1503 if (ir->nstcalcenergy > 0)
1505 print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
1506 eprAVER, mdebin, fcd, groups, &(ir->opts), awh.get());
1509 done_mdebin(mdebin);
1514 pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
1517 done_shellfc(fplog, shellfc, step_rel);
1519 if (useReplicaExchange && MASTER(cr))
1521 print_replica_exchange_statistics(fplog, repl_ex);
1524 // Clean up swapcoords
1525 if (ir->eSwapCoords != eswapNO)
1527 finish_swapcoords(ir->swap);
1530 /* IMD cleanup, if bIMD is TRUE. */
1531 IMD_finalize(ir->bIMD, ir->imd);
1533 walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1535 destroy_enerdata(enerd);
1538 global_stat_destroy(gstat);