2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2011,2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 * \brief Implements the integrator for normal molecular dynamics simulations
41 * \author David van der Spoel <david.vanderspoel@icm.uu.se>
42 * \ingroup module_mdrun
54 #include "gromacs/awh/awh.h"
55 #include "gromacs/commandline/filenm.h"
56 #include "gromacs/domdec/collect.h"
57 #include "gromacs/domdec/dlbtiming.h"
58 #include "gromacs/domdec/domdec.h"
59 #include "gromacs/domdec/domdec_network.h"
60 #include "gromacs/domdec/domdec_struct.h"
61 #include "gromacs/domdec/mdsetup.h"
62 #include "gromacs/domdec/partition.h"
63 #include "gromacs/essentialdynamics/edsam.h"
64 #include "gromacs/ewald/pme.h"
65 #include "gromacs/ewald/pme_load_balancing.h"
66 #include "gromacs/fileio/trxio.h"
67 #include "gromacs/gmxlib/network.h"
68 #include "gromacs/gmxlib/nrnb.h"
69 #include "gromacs/gpu_utils/gpu_utils.h"
70 #include "gromacs/imd/imd.h"
71 #include "gromacs/listed_forces/manage_threading.h"
72 #include "gromacs/math/functions.h"
73 #include "gromacs/math/utilities.h"
74 #include "gromacs/math/vec.h"
75 #include "gromacs/math/vectypes.h"
76 #include "gromacs/mdlib/checkpointhandler.h"
77 #include "gromacs/mdlib/compute_io.h"
78 #include "gromacs/mdlib/constr.h"
79 #include "gromacs/mdlib/ebin.h"
80 #include "gromacs/mdlib/enerdata_utils.h"
81 #include "gromacs/mdlib/energyoutput.h"
82 #include "gromacs/mdlib/expanded.h"
83 #include "gromacs/mdlib/force.h"
84 #include "gromacs/mdlib/force_flags.h"
85 #include "gromacs/mdlib/forcerec.h"
86 #include "gromacs/mdlib/md_support.h"
87 #include "gromacs/mdlib/mdatoms.h"
88 #include "gromacs/mdlib/mdoutf.h"
89 #include "gromacs/mdlib/membed.h"
90 #include "gromacs/mdlib/resethandler.h"
91 #include "gromacs/mdlib/sighandler.h"
92 #include "gromacs/mdlib/simulationsignal.h"
93 #include "gromacs/mdlib/stat.h"
94 #include "gromacs/mdlib/stophandler.h"
95 #include "gromacs/mdlib/tgroup.h"
96 #include "gromacs/mdlib/trajectory_writing.h"
97 #include "gromacs/mdlib/update.h"
98 #include "gromacs/mdlib/update_constrain_cuda.h"
99 #include "gromacs/mdlib/vcm.h"
100 #include "gromacs/mdlib/vsite.h"
101 #include "gromacs/mdrunutility/handlerestart.h"
102 #include "gromacs/mdrunutility/multisim.h"
103 #include "gromacs/mdrunutility/printtime.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/mdrunoptions.h"
117 #include "gromacs/mdtypes/observableshistory.h"
118 #include "gromacs/mdtypes/pullhistory.h"
119 #include "gromacs/mdtypes/state.h"
120 #include "gromacs/nbnxm/nbnxm.h"
121 #include "gromacs/pbcutil/mshift.h"
122 #include "gromacs/pbcutil/pbc.h"
123 #include "gromacs/pulling/output.h"
124 #include "gromacs/pulling/pull.h"
125 #include "gromacs/swap/swapcoords.h"
126 #include "gromacs/timing/wallcycle.h"
127 #include "gromacs/timing/walltime_accounting.h"
128 #include "gromacs/topology/atoms.h"
129 #include "gromacs/topology/idef.h"
130 #include "gromacs/topology/mtop_util.h"
131 #include "gromacs/topology/topology.h"
132 #include "gromacs/trajectory/trajectoryframe.h"
133 #include "gromacs/utility/basedefinitions.h"
134 #include "gromacs/utility/cstringutil.h"
135 #include "gromacs/utility/fatalerror.h"
136 #include "gromacs/utility/logger.h"
137 #include "gromacs/utility/real.h"
138 #include "gromacs/utility/smalloc.h"
140 #include "replicaexchange.h"
142 #include "simulator.h"
145 #include "corewrap.h"
148 using gmx::SimulationSignaller;
150 //! Whether the GPU versions of Leap-Frog integrator and LINCS and SHAKE constraints
151 static const bool c_useGpuUpdateConstrain = (getenv("GMX_UPDATE_CONSTRAIN_GPU") != nullptr);
153 void gmx::Simulator::do_md()
155 // TODO Historically, the EM and MD "integrators" used different
156 // names for the t_inputrec *parameter, but these must have the
157 // same name, now that it's a member of a struct. We use this ir
158 // alias to avoid a large ripple of nearly useless changes.
159 // t_inputrec is being replaced by IMdpOptionsProvider, so this
160 // will go away eventually.
161 t_inputrec *ir = inputrec;
162 int64_t step, step_rel;
163 double t, t0 = ir->init_t, lam0[efptNR];
164 gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
165 gmx_bool bNS, bNStList, bStopCM,
166 bFirstStep, bInitStep, bLastStep = FALSE;
167 gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
168 gmx_bool do_ene, do_log, do_verbose;
169 gmx_bool bMasterState;
170 int force_flags, cglo_flags;
171 tensor force_vir = {{0}}, shake_vir = {{0}}, total_vir = {{0}},
172 tmp_vir = {{0}}, pres = {{0}};
175 matrix parrinellorahmanMu, M;
176 gmx_repl_ex_t repl_ex = nullptr;
178 PaddedVector<gmx::RVec> f {};
179 gmx_global_stat_t gstat;
180 t_graph *graph = nullptr;
181 gmx_shellfc_t *shellfc;
182 gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
183 gmx_bool bTemp, bPres, bTrotter;
185 rvec *cbuf = nullptr;
192 real saved_conserved_quantity = 0;
195 char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
197 /* PME load balancing data for GPU kernels */
198 gmx_bool bPMETune = FALSE;
199 gmx_bool bPMETunePrinting = FALSE;
201 bool bInteractiveMDstep = false;
203 /* Domain decomposition could incorrectly miss a bonded
204 interaction, but checking for that requires a global
205 communication stage, which does not otherwise happen in DD
206 code. So we do that alongside the first global energy reduction
207 after a new DD is made. These variables handle whether the
208 check happens, and the result it returns. */
209 bool shouldCheckNumberOfBondedInteractions = false;
210 int totalNumberOfBondedInteractions = -1;
212 SimulationSignals signals;
213 // Most global communnication stages don't propagate mdrun
214 // signals, and will use this object to achieve that.
215 SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
217 if (!mdrunOptions.writeConfout)
219 // This is on by default, and the main known use case for
220 // turning it off is for convenience in benchmarking, which is
221 // something that should not show up in the general user
223 GMX_LOG(mdlog.info).asParagraph().
224 appendText("The -noconfout functionality is deprecated, and may be removed in a future version.");
227 /* md-vv uses averaged full step velocities for T-control
228 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
229 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
230 bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
232 const bool bRerunMD = false;
234 int nstglobalcomm = computeGlobalCommunicationPeriod(mdlog, ir, cr);
235 bGStatEveryStep = (nstglobalcomm == 1);
237 SimulationGroups *groups = &top_global->groups;
239 std::unique_ptr<EssentialDynamics> ed = nullptr;
240 if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
242 /* Initialize essential dynamics sampling */
243 ed = init_edsam(mdlog,
244 opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
247 state_global, observablesHistory,
252 initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda, lam0);
253 Update upd(ir, deform);
254 bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
255 if (startingBehavior != StartingBehavior::RestartWithAppending)
257 pleaseCiteCouplingAlgorithms(fplog, *ir);
260 gmx_mdoutf *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, ir, top_global, oenv, wcycle,
261 StartingBehavior::NewSimulation);
262 gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work, mdoutf_get_fp_dhdl(outf), false);
264 /* Kinetic energy data */
265 std::unique_ptr<gmx_ekindata_t> eKinData = std::make_unique<gmx_ekindata_t>();
266 gmx_ekindata_t *ekind = eKinData.get();
267 init_ekindata(fplog, top_global, &(ir->opts), ekind);
268 /* Copy the cos acceleration to the groups struct */
269 ekind->cosacc.cos_accel = ir->cos_accel;
271 gstat = global_stat_init(ir);
273 /* Check for polarizable models and flexible constraints */
274 shellfc = init_shell_flexcon(fplog,
275 top_global, constr ? constr->numFlexibleConstraints() : 0,
276 ir->nstcalcenergy, DOMAINDECOMP(cr));
279 double io = compute_io(ir, top_global->natoms, *groups, energyOutput.numEnergyTerms(), 1);
280 if ((io > 2000) && MASTER(cr))
283 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
288 // Local state only becomes valid now.
289 std::unique_ptr<t_state> stateInstance;
293 auto mdatoms = mdAtoms->mdatoms();
295 std::unique_ptr<UpdateConstrainCuda> integrator;
297 if (DOMAINDECOMP(cr))
299 dd_init_local_top(*top_global, &top);
301 stateInstance = std::make_unique<t_state>();
302 state = stateInstance.get();
303 dd_init_local_state(cr->dd, state_global, state);
305 /* Distribute the charge groups over the nodes from the master node */
306 dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
307 state_global, *top_global, ir, imdSession,
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);
328 if (c_useGpuUpdateConstrain)
330 GMX_RELEASE_ASSERT(ir->eI == eiMD, "Only md integrator is supported on the GPU.");
331 GMX_RELEASE_ASSERT(ir->etc == etcNO, "Temperature coupling is not supported on the GPU.");
332 GMX_RELEASE_ASSERT(ir->epc == epcNO, "Pressure coupling is not supported on the GPU.");
333 GMX_RELEASE_ASSERT(!mdatoms->haveVsites, "Virtual sites are not supported on the GPU");
334 GMX_LOG(mdlog.info).asParagraph().
335 appendText("Using CUDA GPU-based update and constraints module.");
336 integrator = std::make_unique<UpdateConstrainCuda>(state->natoms, *ir, *top_global);
337 integrator->set(top.idef, *mdatoms);
339 set_pbc(&pbc, epbcXYZ, state->box);
340 integrator->setPbc(&pbc);
345 if (fr->nbv->useGpu())
347 changePinningPolicy(&state->x, gmx::PinningPolicy::PinnedIfSupported);
350 // NOTE: The global state is no longer used at this point.
351 // But state_global is still used as temporary storage space for writing
352 // the global state to file and potentially for replica exchange.
353 // (Global topology should persist.)
355 update_mdatoms(mdatoms, state->lambda[efptMASS]);
359 /* Check nstexpanded here, because the grompp check was broken */
360 if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
362 gmx_fatal(FARGS, "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
364 init_expanded_ensemble(startingBehavior != StartingBehavior::NewSimulation,
370 if (startingBehavior != StartingBehavior::NewSimulation)
372 /* Restore from energy history if appending to output files */
373 if (startingBehavior == StartingBehavior::RestartWithAppending)
375 /* If no history is available (because a checkpoint is from before
376 * it was written) make a new one later, otherwise restore it.
378 if (observablesHistory->energyHistory)
380 energyOutput.restoreFromEnergyHistory(*observablesHistory->energyHistory);
383 else if (observablesHistory->energyHistory)
385 /* We might have read an energy history from checkpoint.
386 * As we are not appending, we want to restart the statistics.
387 * Free the allocated memory and reset the counts.
389 observablesHistory->energyHistory = {};
390 /* We might have read a pull history from checkpoint.
391 * We will still want to keep the statistics, so that the files
392 * can be joined and still be meaningful.
393 * This means that observablesHistory->pullHistory
394 * should not be reset.
398 if (!observablesHistory->energyHistory)
400 observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
402 if (!observablesHistory->pullHistory)
404 observablesHistory->pullHistory = std::make_unique<PullHistory>();
406 /* Set the initial energy history */
407 energyOutput.fillEnergyHistory(observablesHistory->energyHistory.get());
410 preparePrevStepPullCom(ir, pull_work, mdatoms, state, state_global, cr,
411 startingBehavior != StartingBehavior::NewSimulation);
413 // TODO: Remove this by converting AWH into a ForceProvider
414 auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingBehavior != StartingBehavior::NewSimulation,
416 opt2fn("-awh", nfile, fnm), pull_work);
418 const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
419 if (useReplicaExchange && MASTER(cr))
421 repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
424 /* PME tuning is only supported in the Verlet scheme, with PME for
425 * Coulomb. It is not supported with only LJ PME. */
426 bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) &&
427 !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
429 pme_load_balancing_t *pme_loadbal = nullptr;
432 pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box,
433 *fr->ic, *fr->nbv, fr->pmedata, fr->nbv->useGpu(),
437 if (!ir->bContinuation)
439 if (state->flags & (1 << estV))
441 auto v = makeArrayRef(state->v);
442 /* Set the velocities of vsites, shells and frozen atoms to zero */
443 for (i = 0; i < mdatoms->homenr; i++)
445 if (mdatoms->ptype[i] == eptVSite ||
446 mdatoms->ptype[i] == eptShell)
450 else if (mdatoms->cFREEZE)
452 for (m = 0; m < DIM; m++)
454 if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
465 /* Constrain the initial coordinates and velocities */
466 do_constrain_first(fplog, constr, ir, mdatoms, state);
470 /* Construct the virtual sites for the initial configuration */
471 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
472 top.idef.iparams, top.idef.il,
473 fr->ePBC, fr->bMolPBC, cr, state->box);
477 if (ir->efep != efepNO)
479 /* Set free energy calculation frequency as the greatest common
480 * denominator of nstdhdl and repl_ex_nst. */
481 nstfep = ir->fepvals->nstdhdl;
484 nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
486 if (useReplicaExchange)
488 nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
492 /* Be REALLY careful about what flags you set here. You CANNOT assume
493 * this is the first step, since we might be restarting from a checkpoint,
494 * and in that case we should not do any modifications to the state.
496 bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
498 // When restarting from a checkpoint, it can be appropriate to
499 // initialize ekind from quantities in the checkpoint. Otherwise,
500 // compute_globals must initialize ekind before the simulation
501 // starts/restarts. However, only the master rank knows what was
502 // found in the checkpoint file, so we have to communicate in
503 // order to coordinate the restart.
505 // TODO Consider removing this communication if/when checkpoint
506 // reading directly follows .tpr reading, because all ranks can
507 // agree on hasReadEkinState at that time.
508 bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
511 gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr);
513 if (hasReadEkinState)
515 restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
518 cglo_flags = (CGLO_INITIALIZATION | CGLO_TEMPERATURE | CGLO_GSTAT
519 | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
520 | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
521 | (hasReadEkinState ? CGLO_READEKIN : 0));
523 bSumEkinhOld = FALSE;
525 t_vcm vcm(top_global->groups, *ir);
526 reportComRemovalInfo(fplog, vcm);
528 /* To minimize communication, compute_globals computes the COM velocity
529 * and the kinetic energy for the velocities without COM motion removed.
530 * Thus to get the kinetic energy without the COM contribution, we need
531 * to call compute_globals twice.
533 for (int cgloIteration = 0; cgloIteration < (bStopCM ? 2 : 1); cgloIteration++)
535 int cglo_flags_iteration = cglo_flags;
536 if (bStopCM && cgloIteration == 0)
538 cglo_flags_iteration |= CGLO_STOPCM;
539 cglo_flags_iteration &= ~CGLO_TEMPERATURE;
541 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
542 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
543 constr, &nullSignaller, state->box,
544 &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
545 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
547 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
548 top_global, &top, state,
549 &shouldCheckNumberOfBondedInteractions);
550 if (ir->eI == eiVVAK)
552 /* a second call to get the half step temperature initialized as well */
553 /* we do the same call as above, but turn the pressure off -- internally to
554 compute_globals, this is recognized as a velocity verlet half-step
555 kinetic energy calculation. This minimized excess variables, but
556 perhaps loses some logic?*/
558 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
559 nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
560 constr, &nullSignaller, state->box,
561 nullptr, &bSumEkinhOld,
562 cglo_flags & ~CGLO_PRESSURE);
565 /* Calculate the initial half step temperature, and save the ekinh_old */
566 if (startingBehavior == StartingBehavior::NewSimulation)
568 for (i = 0; (i < ir->opts.ngtc); i++)
570 copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old);
574 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
575 temperature control */
576 auto trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
580 if (!ir->bContinuation)
582 if (constr && ir->eConstrAlg == econtLINCS)
585 "RMS relative constraint deviation after constraining: %.2e\n",
588 if (EI_STATE_VELOCITY(ir->eI))
590 real temp = enerd->term[F_TEMP];
593 /* Result of Ekin averaged over velocities of -half
594 * and +half step, while we only have -half step here.
598 fprintf(fplog, "Initial temperature: %g K\n", temp);
603 fprintf(stderr, "starting mdrun '%s'\n",
604 *(top_global->name));
607 sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
611 sprintf(tbuf, "%s", "infinite");
613 if (ir->init_step > 0)
615 fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
616 gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
617 gmx_step_str(ir->init_step, sbuf2),
618 ir->init_step*ir->delta_t);
622 fprintf(stderr, "%s steps, %s ps.\n",
623 gmx_step_str(ir->nsteps, sbuf), tbuf);
625 fprintf(fplog, "\n");
628 walltime_accounting_start_time(walltime_accounting);
629 wallcycle_start(wcycle, ewcRUN);
630 print_start(fplog, cr, walltime_accounting, "mdrun");
633 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
634 int chkpt_ret = fcCheckPointParallel( cr->nodeid,
638 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
642 /***********************************************************
646 ************************************************************/
649 /* Skip the first Nose-Hoover integration when we get the state from tpx */
650 bInitStep = startingBehavior == StartingBehavior::NewSimulation || EI_VV(ir->eI);
651 bSumEkinhOld = FALSE;
653 bNeedRepartition = FALSE;
655 bool simulationsShareState = false;
656 int nstSignalComm = nstglobalcomm;
658 // TODO This implementation of ensemble orientation restraints is nasty because
659 // a user can't just do multi-sim with single-sim orientation restraints.
660 bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
661 bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
663 // Replica exchange, ensemble restraints and AWH need all
664 // simulations to remain synchronized, so they need
665 // checkpoints and stop conditions to act on the same step, so
666 // the propagation of such signals must take place between
667 // simulations, not just within simulations.
668 // TODO: Make algorithm initializers set these flags.
669 simulationsShareState = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
671 if (simulationsShareState)
673 // Inter-simulation signal communication does not need to happen
674 // often, so we use a minimum of 200 steps to reduce overhead.
675 const int c_minimumInterSimulationSignallingInterval = 200;
676 nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
680 auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
681 compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
682 MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
683 mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
685 auto checkpointHandler = std::make_unique<CheckpointHandler>(
686 compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]),
687 simulationsShareState, ir->nstlist == 0, MASTER(cr),
688 mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
690 const bool resetCountersIsLocal = true;
691 auto resetHandler = std::make_unique<ResetHandler>(
692 compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]), !resetCountersIsLocal,
693 ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
694 mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
696 const DDBalanceRegionHandler ddBalanceRegionHandler(cr);
698 step = ir->init_step;
701 // TODO extract this to new multi-simulation module
702 if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
704 if (!multisim_int_all_are_equal(ms, ir->nsteps))
706 GMX_LOG(mdlog.warning).appendText(
707 "Note: The number of steps is not consistent across multi simulations,\n"
708 "but we are proceeding anyway!");
710 if (!multisim_int_all_are_equal(ms, ir->init_step))
712 GMX_LOG(mdlog.warning).appendText(
713 "Note: The initial step is not consistent across multi simulations,\n"
714 "but we are proceeding anyway!");
718 /* and stop now if we should */
719 bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
723 /* Determine if this is a neighbor search step */
724 bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
726 if (bPMETune && bNStList)
728 /* PME grid + cut-off optimization with GPUs or PME nodes */
729 pme_loadbal_do(pme_loadbal, cr,
730 (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
738 wallcycle_start(wcycle, ewcSTEP);
740 bLastStep = (step_rel == ir->nsteps);
741 t = t0 + step*ir->delta_t;
743 // TODO Refactor this, so that nstfep does not need a default value of zero
744 if (ir->efep != efepNO || ir->bSimTemp)
746 /* find and set the current lambdas */
747 setCurrentLambdasLocal(step, ir->fepvals, lam0, state);
749 bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
750 bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
751 bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
752 && (ir->bExpanded) && (step > 0) &&
753 (startingBehavior == StartingBehavior::NewSimulation));
756 bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
757 do_per_step(step, replExParams.exchangeInterval));
759 if (doSimulatedAnnealing)
761 update_annealing_target_temp(ir, t, &upd);
764 /* Stop Center of Mass motion */
765 bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
767 /* Determine whether or not to do Neighbour Searching */
768 bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
770 bLastStep = bLastStep || stopHandler->stoppingAfterCurrentStep(bNS);
772 /* do_log triggers energy and virial calculation. Because this leads
773 * to different code paths, forces can be different. Thus for exact
774 * continuation we should avoid extra log output.
775 * Note that the || bLastStep can result in non-exact continuation
776 * beyond the last step. But we don't consider that to be an issue.
779 (do_per_step(step, ir->nstlog) ||
780 (bFirstStep && startingBehavior == StartingBehavior::NewSimulation) ||
782 do_verbose = mdrunOptions.verbose &&
783 (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
785 if (bNS && !(bFirstStep && ir->bContinuation))
787 bMasterState = FALSE;
788 /* Correct the new box if it is too skewed */
789 if (inputrecDynamicBox(ir))
791 if (correct_box(fplog, step, state->box, graph))
796 if (DOMAINDECOMP(cr) && bMasterState)
798 dd_collect_state(cr->dd, state, state_global);
801 if (DOMAINDECOMP(cr))
803 /* Repartition the domain decomposition */
804 dd_partition_system(fplog, mdlog, step, cr,
805 bMasterState, nstglobalcomm,
806 state_global, *top_global, ir, imdSession,
808 state, &f, mdAtoms, &top, fr,
811 do_verbose && !bPMETunePrinting);
812 shouldCheckNumberOfBondedInteractions = true;
813 upd.setNumAtoms(state->natoms);
817 if (MASTER(cr) && do_log)
819 energyOutput.printHeader(fplog, step, t); /* can we improve the information printed here? */
822 if (ir->efep != efepNO)
824 update_mdatoms(mdatoms, state->lambda[efptMASS]);
830 /* We need the kinetic energy at minus the half step for determining
831 * the full step kinetic energy and possibly for T-coupling.*/
832 /* This may not be quite working correctly yet . . . . */
833 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
834 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
835 constr, &nullSignaller, state->box,
836 &totalNumberOfBondedInteractions, &bSumEkinhOld,
837 CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
838 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
839 top_global, &top, state,
840 &shouldCheckNumberOfBondedInteractions);
842 clear_mat(force_vir);
844 checkpointHandler->decideIfCheckpointingThisStep(bNS, bFirstStep, bLastStep);
846 /* Determine the energy and pressure:
847 * at nstcalcenergy steps and at energy output steps (set below).
849 if (EI_VV(ir->eI) && (!bInitStep))
851 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
852 bCalcVir = bCalcEnerStep ||
853 (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
857 bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
858 bCalcVir = bCalcEnerStep ||
859 (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
861 bCalcEner = bCalcEnerStep;
863 do_ene = (do_per_step(step, ir->nstenergy) || bLastStep);
865 if (do_ene || do_log || bDoReplEx)
871 /* Do we need global communication ? */
872 bGStat = (bCalcVir || bCalcEner || bStopCM ||
873 do_per_step(step, nstglobalcomm) ||
874 (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
876 force_flags = (GMX_FORCE_STATECHANGED |
877 ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0) |
878 GMX_FORCE_ALLFORCES |
879 (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
880 (bCalcEner ? GMX_FORCE_ENERGY : 0) |
881 (bDoFEP ? GMX_FORCE_DHDL : 0)
886 /* Now is the time to relax the shells */
887 relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
888 enforcedRotation, step,
889 ir, imdSession, pull_work, bNS, force_flags, &top,
891 state, f.arrayRefWithPadding(), force_vir, mdatoms,
893 shellfc, fr, ppForceWorkload, t, mu_tot,
895 ddBalanceRegionHandler);
899 /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
900 (or the AWH update will be performed twice for one step when continuing). It would be best to
901 call this update function from do_md_trajectory_writing but that would occur after do_force.
902 One would have to divide the update_awh function into one function applying the AWH force
903 and one doing the AWH bias update. The update AWH bias function could then be called after
904 do_md_trajectory_writing (then containing update_awh_history).
905 The checkpointing will in the future probably moved to the start of the md loop which will
906 rid of this issue. */
907 if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
909 awh->updateHistory(state_global->awhHistory.get());
912 /* The coordinates (x) are shifted (to get whole molecules)
914 * This is parallellized as well, and does communication too.
915 * Check comments in sim_util.c
917 do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession,
919 step, nrnb, wcycle, &top,
920 state->box, state->x.arrayRefWithPadding(), &state->hist,
921 f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
922 state->lambda, graph,
923 fr, ppForceWorkload, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
924 (bNS ? GMX_FORCE_NS : 0) | force_flags,
925 ddBalanceRegionHandler);
928 if (EI_VV(ir->eI) && startingBehavior == StartingBehavior::NewSimulation)
929 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
931 rvec *vbuf = nullptr;
933 wallcycle_start(wcycle, ewcUPDATE);
934 if (ir->eI == eiVV && bInitStep)
936 /* if using velocity verlet with full time step Ekin,
937 * take the first half step only to compute the
938 * virial for the first step. From there,
939 * revert back to the initial coordinates
940 * so that the input is actually the initial step.
942 snew(vbuf, state->natoms);
943 copy_rvecn(state->v.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
947 /* this is for NHC in the Ekin(t+dt/2) version of vv */
948 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
951 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
952 ekind, M, &upd, etrtVELOCITY1,
955 wallcycle_stop(wcycle, ewcUPDATE);
956 constrain_velocities(step, nullptr,
960 bCalcVir, do_log, do_ene);
961 wallcycle_start(wcycle, ewcUPDATE);
962 /* if VV, compute the pressure and constraints */
963 /* For VV2, we strictly only need this if using pressure
964 * control, but we really would like to have accurate pressures
966 * Think about ways around this in the future?
967 * For now, keep this choice in comments.
969 /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
970 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
972 bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK));
973 if (bCalcEner && ir->eI == eiVVAK)
977 /* for vv, the first half of the integration actually corresponds to the previous step.
978 So we need information from the last step in the first half of the integration */
979 if (bGStat || do_per_step(step-1, nstglobalcomm))
981 wallcycle_stop(wcycle, ewcUPDATE);
982 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
983 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
984 constr, &nullSignaller, state->box,
985 &totalNumberOfBondedInteractions, &bSumEkinhOld,
986 (bGStat ? CGLO_GSTAT : 0)
987 | (bCalcEner ? CGLO_ENERGY : 0)
988 | (bTemp ? CGLO_TEMPERATURE : 0)
989 | (bPres ? CGLO_PRESSURE : 0)
990 | (bPres ? CGLO_CONSTRAINT : 0)
991 | (bStopCM ? CGLO_STOPCM : 0)
992 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
995 /* explanation of above:
996 a) We compute Ekin at the full time step
997 if 1) we are using the AveVel Ekin, and it's not the
998 initial step, or 2) if we are using AveEkin, but need the full
999 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1000 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1001 EkinAveVel because it's needed for the pressure */
1002 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1003 top_global, &top, state,
1004 &shouldCheckNumberOfBondedInteractions);
1005 wallcycle_start(wcycle, ewcUPDATE);
1007 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1012 m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
1013 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
1015 /* TODO This is only needed when we're about to write
1016 * a checkpoint, because we use it after the restart
1017 * (in a kludge?). But what should we be doing if
1018 * the startingBehavior is NewSimulation or bInitStep are true? */
1019 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1021 copy_mat(shake_vir, state->svir_prev);
1022 copy_mat(force_vir, state->fvir_prev);
1024 if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
1026 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1027 enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
1028 enerd->term[F_EKIN] = trace(ekind->ekin);
1031 else if (bExchanged)
1033 wallcycle_stop(wcycle, ewcUPDATE);
1034 /* We need the kinetic energy at minus the half step for determining
1035 * the full step kinetic energy and possibly for T-coupling.*/
1036 /* This may not be quite working correctly yet . . . . */
1037 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1038 wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
1039 constr, &nullSignaller, state->box,
1040 nullptr, &bSumEkinhOld,
1041 CGLO_GSTAT | CGLO_TEMPERATURE);
1042 wallcycle_start(wcycle, ewcUPDATE);
1045 /* if it's the initial step, we performed this first step just to get the constraint virial */
1046 if (ir->eI == eiVV && bInitStep)
1048 copy_rvecn(vbuf, state->v.rvec_array(), 0, state->natoms);
1051 wallcycle_stop(wcycle, ewcUPDATE);
1054 /* compute the conserved quantity */
1057 saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
1060 last_ekin = enerd->term[F_EKIN];
1062 if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres))
1064 saved_conserved_quantity -= enerd->term[F_DISPCORR];
1066 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1067 if (ir->efep != efepNO)
1069 sum_dhdl(enerd, state->lambda, ir->fepvals);
1073 /* ######## END FIRST UPDATE STEP ############## */
1074 /* ######## If doing VV, we now have v(dt) ###### */
1077 /* perform extended ensemble sampling in lambda - we don't
1078 actually move to the new state before outputting
1079 statistics, but if performing simulated tempering, we
1080 do update the velocities and the tau_t. */
1082 lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, state->v.rvec_array(), mdatoms);
1083 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1086 copy_df_history(state_global->dfhist, state->dfhist);
1090 /* Now we have the energies and forces corresponding to the
1091 * coordinates at time t. We must output all of this before
1094 do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
1095 ir, state, state_global, observablesHistory,
1097 outf, energyOutput, ekind, f,
1098 checkpointHandler->isCheckpointingStep(),
1099 bRerunMD, bLastStep,
1100 mdrunOptions.writeConfout,
1102 /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1103 bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x.rvec_array(), t);
1105 /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1106 if (startingBehavior != StartingBehavior::NewSimulation &&
1107 (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
1109 copy_mat(state->svir_prev, shake_vir);
1110 copy_mat(state->fvir_prev, force_vir);
1113 stopHandler->setSignal();
1114 resetHandler->setSignal(walltime_accounting);
1116 if (bGStat || !PAR(cr))
1118 /* In parallel we only have to check for checkpointing in steps
1119 * where we do global communication,
1120 * otherwise the other nodes don't know.
1122 checkpointHandler->setSignal(walltime_accounting);
1125 /* ######### START SECOND UPDATE STEP ################# */
1127 /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1130 if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
1132 gmx_bool bIfRandomize;
1133 bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
1134 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1135 if (constr && bIfRandomize)
1137 constrain_velocities(step, nullptr,
1141 bCalcVir, do_log, do_ene);
1144 /* Box is changed in update() when we do pressure coupling,
1145 * but we should still use the old box for energy corrections and when
1146 * writing it to the energy file, so it matches the trajectory files for
1147 * the same timestep above. Make a copy in a separate array.
1149 copy_mat(state->box, lastbox);
1153 wallcycle_start(wcycle, ewcUPDATE);
1154 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1157 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
1158 /* We can only do Berendsen coupling after we have summed
1159 * the kinetic energy or virial. Since the happens
1160 * in global_state after update, we should only do it at
1161 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1166 update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
1167 update_pcouple_before_coordinates(fplog, step, ir, state,
1168 parrinellorahmanMu, M,
1174 /* velocity half-step update */
1175 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1176 ekind, M, &upd, etrtVELOCITY2,
1180 /* Above, initialize just copies ekinh into ekin,
1181 * it doesn't copy position (for VV),
1182 * and entire integrator for MD.
1185 if (ir->eI == eiVVAK)
1187 /* We probably only need md->homenr, not state->natoms */
1188 if (state->natoms > cbuf_nalloc)
1190 cbuf_nalloc = state->natoms;
1191 srenew(cbuf, cbuf_nalloc);
1193 copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
1196 if (c_useGpuUpdateConstrain)
1198 integrator->copyCoordinatesToGpu(state->x.rvec_array());
1199 integrator->copyVelocitiesToGpu(state->v.rvec_array());
1200 integrator->copyForcesToGpu(as_rvec_array(f.data()));
1202 // This applies Leap-Frog, LINCS and SETTLE in a succession
1203 integrator->integrate(ir->delta_t, true, bCalcVir, shake_vir);
1205 integrator->copyCoordinatesFromGpu(state->x.rvec_array());
1206 integrator->copyVelocitiesFromGpu(state->v.rvec_array());
1210 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1211 ekind, M, &upd, etrtPOSITION, cr, constr);
1213 wallcycle_stop(wcycle, ewcUPDATE);
1215 constrain_coordinates(step, &dvdl_constr, state,
1218 bCalcVir, do_log, do_ene);
1220 update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
1221 cr, nrnb, wcycle, &upd, constr, do_log, do_ene);
1222 finish_update(ir, mdatoms,
1224 nrnb, wcycle, &upd, constr);
1227 if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
1229 updatePrevStepPullCom(pull_work, state);
1232 if (ir->eI == eiVVAK)
1234 /* erase F_EKIN and F_TEMP here? */
1235 /* just compute the kinetic energy at the half step to perform a trotter step */
1236 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1237 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1238 constr, &nullSignaller, lastbox,
1239 nullptr, &bSumEkinhOld,
1240 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
1242 wallcycle_start(wcycle, ewcUPDATE);
1243 trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
1244 /* now we know the scaling, we can compute the positions again again */
1245 copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
1247 update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
1248 ekind, M, &upd, etrtPOSITION, cr, constr);
1249 wallcycle_stop(wcycle, ewcUPDATE);
1251 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1252 /* are the small terms in the shake_vir here due
1253 * to numerical errors, or are they important
1254 * physically? I'm thinking they are just errors, but not completely sure.
1255 * For now, will call without actually constraining, constr=NULL*/
1256 finish_update(ir, mdatoms,
1258 nrnb, wcycle, &upd, nullptr);
1262 /* this factor or 2 correction is necessary
1263 because half of the constraint force is removed
1264 in the vv step, so we have to double it. See
1265 the Redmine issue #1255. It is not yet clear
1266 if the factor of 2 is exact, or just a very
1267 good approximation, and this will be
1268 investigated. The next step is to see if this
1269 can be done adding a dhdl contribution from the
1270 rattle step, but this is somewhat more
1271 complicated with the current code. Will be
1272 investigated, hopefully for 4.6.3. However,
1273 this current solution is much better than
1274 having it completely wrong.
1276 enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
1280 enerd->term[F_DVDL_CONSTR] += dvdl_constr;
1283 if (vsite != nullptr)
1285 wallcycle_start(wcycle, ewcVSITECONSTR);
1286 if (graph != nullptr)
1288 shift_self(graph, state->box, state->x.rvec_array());
1290 construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
1291 top.idef.iparams, top.idef.il,
1292 fr->ePBC, fr->bMolPBC, cr, state->box);
1294 if (graph != nullptr)
1296 unshift_self(graph, state->box, state->x.rvec_array());
1298 wallcycle_stop(wcycle, ewcVSITECONSTR);
1301 /* ############## IF NOT VV, Calculate globals HERE ############ */
1302 /* With Leap-Frog we can skip compute_globals at
1303 * non-communication steps, but we need to calculate
1304 * the kinetic energy one step before communication.
1307 // Organize to do inter-simulation signalling on steps if
1308 // and when algorithms require it.
1309 bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
1311 if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
1313 // Since we're already communicating at this step, we
1314 // can propagate intra-simulation signals. Note that
1315 // check_nstglobalcomm has the responsibility for
1316 // choosing the value of nstglobalcomm that is one way
1317 // bGStat becomes true, so we can't get into a
1318 // situation where e.g. checkpointing can't be
1320 bool doIntraSimSignal = true;
1321 SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
1323 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, &vcm,
1324 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
1327 &totalNumberOfBondedInteractions, &bSumEkinhOld,
1328 (bGStat ? CGLO_GSTAT : 0)
1329 | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
1330 | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
1331 | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
1332 | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
1334 | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
1336 checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
1337 top_global, &top, state,
1338 &shouldCheckNumberOfBondedInteractions);
1342 /* ############# END CALC EKIN AND PRESSURE ################# */
1344 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1345 the virial that should probably be addressed eventually. state->veta has better properies,
1346 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1347 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1349 if (ir->efep != efepNO && !EI_VV(ir->eI))
1351 /* Sum up the foreign energy and dhdl terms for md and sd.
1352 Currently done every step so that dhdl is correct in the .edr */
1353 sum_dhdl(enerd, state->lambda, ir->fepvals);
1356 update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
1357 pres, force_vir, shake_vir,
1361 /* ################# END UPDATE STEP 2 ################# */
1362 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1364 /* The coordinates (x) were unshifted in update */
1367 /* We will not sum ekinh_old,
1368 * so signal that we still have to do it.
1370 bSumEkinhOld = TRUE;
1375 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1377 /* use the directly determined last velocity, not actually the averaged half steps */
1378 if (bTrotter && ir->eI == eiVV)
1380 enerd->term[F_EKIN] = last_ekin;
1382 enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
1384 if (integratorHasConservedEnergyQuantity(ir))
1388 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
1392 enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
1395 /* ######### END PREPARING EDR OUTPUT ########### */
1401 if (fplog && do_log && bDoExpanded)
1403 /* only needed if doing expanded ensemble */
1404 PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
1405 state_global->dfhist, state->fep_state, ir->nstlog, step);
1409 energyOutput.addDataAtEnergyStep(bDoDHDL, bCalcEnerStep,
1410 t, mdatoms->tmass, enerd, state,
1411 ir->fepvals, ir->expandedvals, lastbox,
1412 shake_vir, force_vir, total_vir, pres,
1413 ekind, mu_tot, constr);
1417 energyOutput.recordNonEnergyStep();
1420 gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
1421 gmx_bool do_or = do_per_step(step, ir->nstorireout);
1423 energyOutput.printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
1424 energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
1425 do_log ? fplog : nullptr,
1431 pull_print_output(pull_work, step, t);
1434 if (do_per_step(step, ir->nstlog))
1436 if (fflush(fplog) != 0)
1438 gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
1444 /* Have to do this part _after_ outputting the logfile and the edr file */
1445 /* Gets written into the state at the beginning of next loop*/
1446 state->fep_state = lamnew;
1448 /* Print the remaining wall clock time for the run */
1449 if (isMasterSimMasterRank(ms, MASTER(cr)) &&
1450 (do_verbose || gmx_got_usr_signal()) &&
1455 fprintf(stderr, "\n");
1457 print_time(stderr, walltime_accounting, step, ir, cr);
1460 /* Ion/water position swapping.
1461 * Not done in last step since trajectory writing happens before this call
1462 * in the MD loop and exchanges would be lost anyway. */
1463 bNeedRepartition = FALSE;
1464 if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
1465 do_per_step(step, ir->swap->nstswap))
1467 bNeedRepartition = do_swapcoords(cr, step, t, ir, swap, wcycle,
1468 as_rvec_array(state->x.data()),
1470 MASTER(cr) && mdrunOptions.verbose,
1473 if (bNeedRepartition && DOMAINDECOMP(cr))
1475 dd_collect_state(cr->dd, state, state_global);
1479 /* Replica exchange */
1483 bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
1484 state_global, enerd,
1488 if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
1490 dd_partition_system(fplog, mdlog, step, cr, TRUE, 1,
1491 state_global, *top_global, ir, imdSession,
1493 state, &f, mdAtoms, &top, fr,
1495 nrnb, wcycle, FALSE);
1496 shouldCheckNumberOfBondedInteractions = true;
1497 upd.setNumAtoms(state->natoms);
1503 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1504 /* With all integrators, except VV, we need to retain the pressure
1505 * at the current step for coupling at the next step.
1507 if ((state->flags & (1<<estPRES_PREV)) &&
1509 (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
1511 /* Store the pressure in t_state for pressure coupling
1512 * at the next MD step.
1514 copy_mat(pres, state->pres_prev);
1517 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1519 if ( (membed != nullptr) && (!bLastStep) )
1521 rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
1524 cycles = wallcycle_stop(wcycle, ewcSTEP);
1525 if (DOMAINDECOMP(cr) && wcycle)
1527 dd_cycles_add(cr->dd, cycles, ddCyclStep);
1530 /* increase the MD step number */
1534 resetHandler->resetCounters(
1535 step, step_rel, mdlog, fplog, cr, fr->nbv.get(),
1536 nrnb, fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
1538 /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1539 imdSession->updateEnergyRecordAndSendPositionsAndEnergies(bInteractiveMDstep, step, bCalcEner);
1542 /* End of main MD loop */
1544 /* Closing TNG files can include compressing data. Therefore it is good to do that
1545 * before stopping the time measurements. */
1546 mdoutf_tng_close(outf);
1548 /* Stop measuring walltime */
1549 walltime_accounting_end_time(walltime_accounting);
1551 if (!thisRankHasDuty(cr, DUTY_PME))
1553 /* Tell the PME only node to finish */
1554 gmx_pme_send_finish(cr);
1559 if (ir->nstcalcenergy > 0)
1561 energyOutput.printAnnealingTemperatures(fplog, groups, &(ir->opts));
1562 energyOutput.printAverages(fplog, groups);
1569 pme_loadbal_done(pme_loadbal, fplog, mdlog, fr->nbv->useGpu());
1572 done_shellfc(fplog, shellfc, step_rel);
1574 if (useReplicaExchange && MASTER(cr))
1576 print_replica_exchange_statistics(fplog, repl_ex);
1579 walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
1581 global_stat_destroy(gstat);