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, 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
56 #include "gromacs/awh/awh.h"
57 #include "gromacs/commandline/filenm.h"
58 #include "gromacs/compat/make_unique.h"
59 #include "gromacs/domdec/collect.h"
60 #include "gromacs/domdec/domdec.h"
61 #include "gromacs/domdec/domdec_network.h"
62 #include "gromacs/domdec/domdec_struct.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/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/repl_ex.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/tgroup.h"
99 #include "gromacs/mdlib/trajectory_writing.h"
100 #include "gromacs/mdlib/update.h"
101 #include "gromacs/mdlib/vcm.h"
102 #include "gromacs/mdlib/vsite.h"
103 #include "gromacs/mdtypes/awh-history.h"
104 #include "gromacs/mdtypes/awh-params.h"
105 #include "gromacs/mdtypes/commrec.h"
106 #include "gromacs/mdtypes/df_history.h"
107 #include "gromacs/mdtypes/energyhistory.h"
108 #include "gromacs/mdtypes/fcdata.h"
109 #include "gromacs/mdtypes/forcerec.h"
110 #include "gromacs/mdtypes/group.h"
111 #include "gromacs/mdtypes/inputrec.h"
112 #include "gromacs/mdtypes/interaction_const.h"
113 #include "gromacs/mdtypes/md_enums.h"
114 #include "gromacs/mdtypes/mdatom.h"
115 #include "gromacs/mdtypes/observableshistory.h"
116 #include "gromacs/mdtypes/state.h"
117 #include "gromacs/pbcutil/mshift.h"
118 #include "gromacs/pbcutil/pbc.h"
119 #include "gromacs/pulling/pull.h"
120 #include "gromacs/swap/swapcoords.h"
121 #include "gromacs/timing/wallcycle.h"
122 #include "gromacs/timing/walltime_accounting.h"
123 #include "gromacs/topology/atoms.h"
124 #include "gromacs/topology/idef.h"
125 #include "gromacs/topology/mtop_util.h"
126 #include "gromacs/topology/topology.h"
127 #include "gromacs/trajectory/trajectoryframe.h"
128 #include "gromacs/utility/basedefinitions.h"
129 #include "gromacs/utility/cstringutil.h"
130 #include "gromacs/utility/fatalerror.h"
131 #include "gromacs/utility/logger.h"
132 #include "gromacs/utility/real.h"
133 #include "gromacs/utility/smalloc.h"
135 #include "integrator.h"
138 #include "corewrap.h"
141 using gmx::SimulationSignaller
;
143 //! Resets all the counters.
144 static void reset_all_counters(FILE *fplog
, const gmx::MDLogger
&mdlog
, t_commrec
*cr
,
146 gmx_int64_t
*step_rel
, t_inputrec
*ir
,
147 gmx_wallcycle_t wcycle
, t_nrnb
*nrnb
,
148 gmx_walltime_accounting_t walltime_accounting
,
149 struct nonbonded_verlet_t
*nbv
,
150 struct gmx_pme_t
*pme
)
152 char sbuf
[STEPSTRSIZE
];
154 /* Reset all the counters related to performance over the run */
155 GMX_LOG(mdlog
.warning
).asParagraph().appendTextFormatted(
156 "step %s: resetting all time and cycle counters",
157 gmx_step_str(step
, sbuf
));
161 nbnxn_gpu_reset_timings(nbv
);
164 if (pme_gpu_task_enabled(pme
))
166 pme_gpu_reset_timings(pme
);
169 if (use_GPU(nbv
) || pme_gpu_task_enabled(pme
))
174 wallcycle_stop(wcycle
, ewcRUN
);
175 wallcycle_reset_all(wcycle
);
176 if (DOMAINDECOMP(cr
))
178 reset_dd_statistics_counters(cr
->dd
);
181 ir
->init_step
+= *step_rel
;
182 ir
->nsteps
-= *step_rel
;
184 wallcycle_start(wcycle
, ewcRUN
);
185 walltime_accounting_start(walltime_accounting
);
186 print_date_and_time(fplog
, cr
->nodeid
, "Restarted time", gmx_gettime());
189 /*! \brief Copy the state from \p rerunFrame to \p globalState and, if requested, construct vsites
191 * \param[in] rerunFrame The trajectory frame to compute energy/forces for
192 * \param[in,out] globalState The global state container
193 * \param[in] constructVsites When true, vsite coordinates are constructed
194 * \param[in] vsite Vsite setup, can be nullptr when \p constructVsites = false
195 * \param[in] idef Topology parameters, used for constructing vsites
196 * \param[in] timeStep Time step, used for constructing vsites
197 * \param[in] forceRec Force record, used for constructing vsites
198 * \param[in,out] graph The molecular graph, used for constructing vsites when != nullptr
199 * \param[in,out] warnWhenNoV When true, issue a warning when no velocities are present in \p rerunFrame; is set to false when a warning was issued
201 static void prepareRerunState(const t_trxframe
&rerunFrame
,
202 t_state
*globalState
,
203 bool constructVsites
,
204 const gmx_vsite_t
*vsite
,
207 const t_forcerec
&forceRec
,
209 gmx_bool
*warnWhenNoV
)
211 for (int i
= 0; i
< globalState
->natoms
; i
++)
213 copy_rvec(rerunFrame
.x
[i
], globalState
->x
[i
]);
217 for (int i
= 0; i
< globalState
->natoms
; i
++)
219 copy_rvec(rerunFrame
.v
[i
], globalState
->v
[i
]);
224 for (int i
= 0; i
< globalState
->natoms
; i
++)
226 clear_rvec(globalState
->v
[i
]);
230 fprintf(stderr
, "\nWARNING: Some frames do not contain velocities.\n"
231 " Ekin, temperature and pressure are incorrect,\n"
232 " the virial will be incorrect when constraints are present.\n"
234 *warnWhenNoV
= FALSE
;
237 copy_mat(rerunFrame
.box
, globalState
->box
);
241 GMX_ASSERT(vsite
, "Need valid vsite for constructing vsites");
245 /* Following is necessary because the graph may get out of sync
246 * with the coordinates if we only have every N'th coordinate set
248 mk_mshift(nullptr, graph
, forceRec
.ePBC
, globalState
->box
, as_rvec_array(globalState
->x
.data()));
249 shift_self(graph
, globalState
->box
, as_rvec_array(globalState
->x
.data()));
251 construct_vsites(vsite
, as_rvec_array(globalState
->x
.data()), timeStep
, as_rvec_array(globalState
->v
.data()),
252 idef
.iparams
, idef
.il
,
253 forceRec
.ePBC
, forceRec
.bMolPBC
, nullptr, globalState
->box
);
256 unshift_self(graph
, globalState
->box
, as_rvec_array(globalState
->x
.data()));
261 void gmx::Integrator::do_md()
263 // TODO Historically, the EM and MD "integrators" used different
264 // names for the t_inputrec *parameter, but these must have the
265 // same name, now that it's a member of a struct. We use this ir
266 // alias to avoid a large ripple of nearly useless changes.
267 // t_inputrec is being replaced by IMdpOptionsProvider, so this
268 // will go away eventually.
269 t_inputrec
*ir
= inputrec
;
270 gmx_mdoutf
*outf
= nullptr;
271 gmx_int64_t step
, step_rel
;
273 double t
, t0
, lam0
[efptNR
];
274 gmx_bool bGStatEveryStep
, bGStat
, bCalcVir
, bCalcEnerStep
, bCalcEner
;
275 gmx_bool bNS
, bNStList
, bSimAnn
, bStopCM
,
276 bFirstStep
, bInitStep
, bLastStep
= FALSE
;
277 gmx_bool bDoDHDL
= FALSE
, bDoFEP
= FALSE
, bDoExpanded
= FALSE
;
278 gmx_bool do_ene
, do_log
, do_verbose
, bRerunWarnNoV
= TRUE
,
279 bForceUpdate
= FALSE
, bCPT
;
280 gmx_bool bMasterState
;
281 int force_flags
, cglo_flags
;
282 tensor force_vir
, shake_vir
, total_vir
, tmp_vir
, pres
;
287 matrix parrinellorahmanMu
, M
;
289 gmx_repl_ex_t repl_ex
= nullptr;
292 t_mdebin
*mdebin
= nullptr;
293 gmx_enerdata_t
*enerd
;
294 PaddedRVecVector f
{};
295 gmx_global_stat_t gstat
;
296 gmx_update_t
*upd
= nullptr;
297 t_graph
*graph
= nullptr;
298 gmx_groups_t
*groups
;
299 gmx_ekindata_t
*ekind
;
300 gmx_shellfc_t
*shellfc
;
301 gmx_bool bSumEkinhOld
, bDoReplEx
, bExchanged
, bNeedRepartition
;
302 gmx_bool bResetCountersHalfMaxH
= FALSE
;
303 gmx_bool bTemp
, bPres
, bTrotter
;
305 rvec
*cbuf
= nullptr;
312 real saved_conserved_quantity
= 0;
316 char sbuf
[STEPSTRSIZE
], sbuf2
[STEPSTRSIZE
];
317 int handled_stop_condition
= gmx_stop_cond_none
; /* compare to get_stop_condition*/
320 /* PME load balancing data for GPU kernels */
321 pme_load_balancing_t
*pme_loadbal
= nullptr;
322 gmx_bool bPMETune
= FALSE
;
323 gmx_bool bPMETunePrinting
= FALSE
;
326 gmx_bool bIMDstep
= FALSE
;
329 /* Temporary addition for FAHCORE checkpointing */
332 /* Domain decomposition could incorrectly miss a bonded
333 interaction, but checking for that requires a global
334 communication stage, which does not otherwise happen in DD
335 code. So we do that alongside the first global energy reduction
336 after a new DD is made. These variables handle whether the
337 check happens, and the result it returns. */
338 bool shouldCheckNumberOfBondedInteractions
= false;
339 int totalNumberOfBondedInteractions
= -1;
341 SimulationSignals signals
;
342 // Most global communnication stages don't propagate mdrun
343 // signals, and will use this object to achieve that.
344 SimulationSignaller
nullSignaller(nullptr, nullptr, nullptr, false, false);
346 if (!mdrunOptions
.writeConfout
)
348 // This is on by default, and the main known use case for
349 // turning it off is for convenience in benchmarking, which is
350 // something that should not show up in the general user
352 GMX_LOG(mdlog
.info
).asParagraph().
353 appendText("The -noconfout functionality is deprecated, and may be removed in a future version.");
355 if (mdrunOptions
.timingOptions
.resetHalfway
)
357 GMX_LOG(mdlog
.info
).asParagraph().
358 appendText("The -resethway functionality is deprecated, and may be removed in a future version.");
361 /* Signal to reset the counters half the simulation steps. */
362 wcycle_set_reset_counters(wcycle
, ir
->nsteps
/2);
364 /* Signal to reset the counters halfway the simulation time. */
365 bResetCountersHalfMaxH
= (mdrunOptions
.maximumHoursToRun
> 0);
368 /* md-vv uses averaged full step velocities for T-control
369 md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
370 md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
371 bTrotter
= (EI_VV(ir
->eI
) && (inputrecNptTrotter(ir
) || inputrecNphTrotter(ir
) || inputrecNvtTrotter(ir
)));
373 const gmx_bool bRerunMD
= mdrunOptions
.rerun
;
374 int nstglobalcomm
= mdrunOptions
.globalCommunicationInterval
;
378 /* Since we don't know if the frames read are related in any way,
379 * rebuild the neighborlist at every step.
382 ir
->nstcalcenergy
= 1;
386 nstglobalcomm
= check_nstglobalcomm(mdlog
, nstglobalcomm
, ir
);
387 bGStatEveryStep
= (nstglobalcomm
== 1);
391 ir
->nstxout_compressed
= 0;
393 groups
= &top_global
->groups
;
395 gmx_edsam
*ed
= nullptr;
396 if (opt2bSet("-ei", nfile
, fnm
) || observablesHistory
->edsamHistory
!= nullptr)
398 /* Initialize essential dynamics sampling */
399 ed
= init_edsam(opt2fn_null("-ei", nfile
, fnm
), opt2fn("-eo", nfile
, fnm
),
402 state_global
, observablesHistory
,
403 oenv
, mdrunOptions
.continuationOptions
.appendFiles
);
406 if (ir
->eSwapCoords
!= eswapNO
)
408 /* Initialize ion swapping code */
409 init_swapcoords(fplog
, ir
, opt2fn_master("-swap", nfile
, fnm
, cr
),
411 state_global
, observablesHistory
,
412 cr
, oenv
, mdrunOptions
);
416 init_md(fplog
, cr
, outputProvider
, ir
, oenv
, mdrunOptions
,
417 &t
, &t0
, state_global
, lam0
,
418 nrnb
, top_global
, &upd
, deform
,
419 nfile
, fnm
, &outf
, &mdebin
,
420 force_vir
, shake_vir
, mu_tot
, &bSimAnn
, &vcm
, wcycle
);
422 clear_mat(total_vir
);
424 /* Energy terms and groups */
426 init_enerdata(top_global
->groups
.grps
[egcENER
].nr
, ir
->fepvals
->n_lambda
,
429 /* Kinetic energy data */
431 init_ekindata(fplog
, top_global
, &(ir
->opts
), ekind
);
432 /* Copy the cos acceleration to the groups struct */
433 ekind
->cosacc
.cos_accel
= ir
->cos_accel
;
435 gstat
= global_stat_init(ir
);
437 /* Check for polarizable models and flexible constraints */
438 shellfc
= init_shell_flexcon(fplog
,
439 top_global
, constr
? constr
->numFlexibleConstraints() : 0,
440 ir
->nstcalcenergy
, DOMAINDECOMP(cr
));
443 double io
= compute_io(ir
, top_global
->natoms
, groups
, mdebin
->ebin
->nener
, 1);
444 if ((io
> 2000) && MASTER(cr
))
447 "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
452 /* Set up interactive MD (IMD) */
453 init_IMD(ir
, cr
, ms
, top_global
, fplog
, ir
->nstcalcenergy
,
454 MASTER(cr
) ? as_rvec_array(state_global
->x
.data()) : nullptr,
455 nfile
, fnm
, oenv
, mdrunOptions
);
457 // Local state only becomes valid now.
458 std::unique_ptr
<t_state
> stateInstance
;
461 if (DOMAINDECOMP(cr
))
463 top
= dd_init_local_top(top_global
);
465 stateInstance
= gmx::compat::make_unique
<t_state
>();
466 state
= stateInstance
.get();
467 dd_init_local_state(cr
->dd
, state_global
, state
);
469 /* Distribute the charge groups over the nodes from the master node */
470 dd_partition_system(fplog
, ir
->init_step
, cr
, TRUE
, 1,
471 state_global
, top_global
, ir
,
472 state
, &f
, mdAtoms
, top
, fr
,
474 nrnb
, nullptr, FALSE
);
475 shouldCheckNumberOfBondedInteractions
= true;
476 update_realloc(upd
, state
->natoms
);
480 state_change_natoms(state_global
, state_global
->natoms
);
481 /* We need to allocate one element extra, since we might use
482 * (unaligned) 4-wide SIMD loads to access rvec entries.
484 f
.resize(gmx::paddedRVecVectorSize(state_global
->natoms
));
485 /* Copy the pointer to the global state */
486 state
= state_global
;
489 mdAlgorithmsSetupAtomData(cr
, ir
, top_global
, top
, fr
,
490 &graph
, mdAtoms
, constr
, vsite
, shellfc
);
492 update_realloc(upd
, state
->natoms
);
495 auto mdatoms
= mdAtoms
->mdatoms();
497 // NOTE: The global state is no longer used at this point.
498 // But state_global is still used as temporary storage space for writing
499 // the global state to file and potentially for replica exchange.
500 // (Global topology should persist.)
502 update_mdatoms(mdatoms
, state
->lambda
[efptMASS
]);
504 const ContinuationOptions
&continuationOptions
= mdrunOptions
.continuationOptions
;
505 bool startingFromCheckpoint
= continuationOptions
.startedFromCheckpoint
;
509 init_expanded_ensemble(startingFromCheckpoint
, ir
, state
->dfhist
);
514 if (startingFromCheckpoint
)
516 /* Update mdebin with energy history if appending to output files */
517 if (continuationOptions
.appendFiles
)
519 restore_energyhistory_from_state(mdebin
, observablesHistory
->energyHistory
.get());
521 else if (observablesHistory
->energyHistory
.get() != nullptr)
523 /* We might have read an energy history from checkpoint.
524 * As we are not appending, we want to restart the statistics.
525 * Free the allocated memory and reset the counts.
527 observablesHistory
->energyHistory
= {};
530 if (observablesHistory
->energyHistory
.get() == nullptr)
532 observablesHistory
->energyHistory
= std::unique_ptr
<energyhistory_t
>(new energyhistory_t
{});
534 /* Set the initial energy history in state by updating once */
535 update_energyhistory(observablesHistory
->energyHistory
.get(), mdebin
);
538 // TODO: Remove this by converting AWH into a ForceProvider
539 auto awh
= prepareAwhModule(fplog
, *ir
, state_global
, cr
, ms
, startingFromCheckpoint
,
541 opt2fn("-awh", nfile
, fnm
), ir
->pull_work
);
543 const bool useReplicaExchange
= (replExParams
.exchangeInterval
> 0);
544 if (useReplicaExchange
&& MASTER(cr
))
546 repl_ex
= init_replica_exchange(fplog
, ms
, top_global
->natoms
, ir
,
549 /* PME tuning is only supported in the Verlet scheme, with PME for
550 * Coulomb. It is not supported with only LJ PME, or for
552 bPMETune
= (mdrunOptions
.tunePme
&& EEL_PME(fr
->ic
->eeltype
) && !bRerunMD
&&
553 !mdrunOptions
.reproducible
&& ir
->cutoff_scheme
!= ecutsGROUP
);
556 pme_loadbal_init(&pme_loadbal
, cr
, mdlog
, ir
, state
->box
,
557 fr
->ic
, fr
->nbv
->listParams
.get(), fr
->pmedata
, use_GPU(fr
->nbv
),
561 if (!ir
->bContinuation
&& !bRerunMD
)
563 if (state
->flags
& (1 << estV
))
565 /* Set the velocities of vsites, shells and frozen atoms to zero */
566 for (i
= 0; i
< mdatoms
->homenr
; i
++)
568 if (mdatoms
->ptype
[i
] == eptVSite
||
569 mdatoms
->ptype
[i
] == eptShell
)
571 clear_rvec(state
->v
[i
]);
573 else if (mdatoms
->cFREEZE
)
575 for (m
= 0; m
< DIM
; m
++)
577 if (ir
->opts
.nFreeze
[mdatoms
->cFREEZE
[i
]][m
])
588 /* Constrain the initial coordinates and velocities */
589 do_constrain_first(fplog
, constr
, ir
, mdatoms
, state
);
593 /* Construct the virtual sites for the initial configuration */
594 construct_vsites(vsite
, as_rvec_array(state
->x
.data()), ir
->delta_t
, nullptr,
595 top
->idef
.iparams
, top
->idef
.il
,
596 fr
->ePBC
, fr
->bMolPBC
, cr
, state
->box
);
600 if (ir
->efep
!= efepNO
)
602 /* Set free energy calculation frequency as the greatest common
603 * denominator of nstdhdl and repl_ex_nst. */
604 nstfep
= ir
->fepvals
->nstdhdl
;
607 nstfep
= gmx_greatest_common_divisor(ir
->expandedvals
->nstexpanded
, nstfep
);
609 if (useReplicaExchange
)
611 nstfep
= gmx_greatest_common_divisor(replExParams
.exchangeInterval
, nstfep
);
615 /* Be REALLY careful about what flags you set here. You CANNOT assume
616 * this is the first step, since we might be restarting from a checkpoint,
617 * and in that case we should not do any modifications to the state.
619 bStopCM
= (ir
->comm_mode
!= ecmNO
&& !ir
->bContinuation
);
621 if (continuationOptions
.haveReadEkin
)
623 restore_ekinstate_from_state(cr
, ekind
, &state_global
->ekinstate
);
626 cglo_flags
= (CGLO_INITIALIZATION
| CGLO_TEMPERATURE
| CGLO_GSTAT
627 | (EI_VV(ir
->eI
) ? CGLO_PRESSURE
: 0)
628 | (EI_VV(ir
->eI
) ? CGLO_CONSTRAINT
: 0)
629 | (continuationOptions
.haveReadEkin
? CGLO_READEKIN
: 0));
631 bSumEkinhOld
= FALSE
;
632 /* To minimize communication, compute_globals computes the COM velocity
633 * and the kinetic energy for the velocities without COM motion removed.
634 * Thus to get the kinetic energy without the COM contribution, we need
635 * to call compute_globals twice.
637 for (int cgloIteration
= 0; cgloIteration
< (bStopCM
? 2 : 1); cgloIteration
++)
639 int cglo_flags_iteration
= cglo_flags
;
640 if (bStopCM
&& cgloIteration
== 0)
642 cglo_flags_iteration
|= CGLO_STOPCM
;
643 cglo_flags_iteration
&= ~CGLO_TEMPERATURE
;
645 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, mdatoms
, nrnb
, vcm
,
646 nullptr, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
647 constr
, &nullSignaller
, state
->box
,
648 &totalNumberOfBondedInteractions
, &bSumEkinhOld
, cglo_flags_iteration
649 | (shouldCheckNumberOfBondedInteractions
? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
: 0));
651 checkNumberOfBondedInteractions(fplog
, cr
, totalNumberOfBondedInteractions
,
652 top_global
, top
, state
,
653 &shouldCheckNumberOfBondedInteractions
);
654 if (ir
->eI
== eiVVAK
)
656 /* a second call to get the half step temperature initialized as well */
657 /* we do the same call as above, but turn the pressure off -- internally to
658 compute_globals, this is recognized as a velocity verlet half-step
659 kinetic energy calculation. This minimized excess variables, but
660 perhaps loses some logic?*/
662 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, mdatoms
, nrnb
, vcm
,
663 nullptr, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
664 constr
, &nullSignaller
, state
->box
,
665 nullptr, &bSumEkinhOld
,
666 cglo_flags
& ~CGLO_PRESSURE
);
669 /* Calculate the initial half step temperature, and save the ekinh_old */
670 if (!continuationOptions
.startedFromCheckpoint
)
672 for (i
= 0; (i
< ir
->opts
.ngtc
); i
++)
674 copy_mat(ekind
->tcstat
[i
].ekinh
, ekind
->tcstat
[i
].ekinh_old
);
678 /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
679 temperature control */
680 trotter_seq
= init_npt_vars(ir
, state
, &MassQ
, bTrotter
);
684 if (!ir
->bContinuation
)
686 if (constr
&& ir
->eConstrAlg
== econtLINCS
)
689 "RMS relative constraint deviation after constraining: %.2e\n",
692 if (EI_STATE_VELOCITY(ir
->eI
))
694 real temp
= enerd
->term
[F_TEMP
];
697 /* Result of Ekin averaged over velocities of -half
698 * and +half step, while we only have -half step here.
702 fprintf(fplog
, "Initial temperature: %g K\n", temp
);
708 fprintf(stderr
, "starting md rerun '%s', reading coordinates from"
709 " input trajectory '%s'\n\n",
710 *(top_global
->name
), opt2fn("-rerun", nfile
, fnm
));
711 if (mdrunOptions
.verbose
)
713 fprintf(stderr
, "Calculated time to finish depends on nsteps from "
714 "run input file,\nwhich may not correspond to the time "
715 "needed to process input trajectory.\n\n");
721 fprintf(stderr
, "starting mdrun '%s'\n",
722 *(top_global
->name
));
725 sprintf(tbuf
, "%8.1f", (ir
->init_step
+ir
->nsteps
)*ir
->delta_t
);
729 sprintf(tbuf
, "%s", "infinite");
731 if (ir
->init_step
> 0)
733 fprintf(stderr
, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
734 gmx_step_str(ir
->init_step
+ir
->nsteps
, sbuf
), tbuf
,
735 gmx_step_str(ir
->init_step
, sbuf2
),
736 ir
->init_step
*ir
->delta_t
);
740 fprintf(stderr
, "%s steps, %s ps.\n",
741 gmx_step_str(ir
->nsteps
, sbuf
), tbuf
);
744 fprintf(fplog
, "\n");
747 walltime_accounting_start(walltime_accounting
);
748 wallcycle_start(wcycle
, ewcRUN
);
749 print_start(fplog
, cr
, walltime_accounting
, "mdrun");
751 /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
753 chkpt_ret
= fcCheckPointParallel( cr
->nodeid
,
757 gmx_fatal( 3, __FILE__
, __LINE__
, "Checkpoint error on step %d\n", 0 );
761 /***********************************************************
765 ************************************************************/
767 /* if rerunMD then read coordinates and velocities from input trajectory */
770 if (getenv("GMX_FORCE_UPDATE"))
778 bLastStep
= !read_first_frame(oenv
, &status
,
779 opt2fn("-rerun", nfile
, fnm
),
780 &rerun_fr
, TRX_NEED_X
| TRX_READ_V
);
781 if (rerun_fr
.natoms
!= top_global
->natoms
)
784 "Number of atoms in trajectory (%d) does not match the "
785 "run input file (%d)\n",
786 rerun_fr
.natoms
, top_global
->natoms
);
788 if (ir
->ePBC
!= epbcNONE
)
792 gmx_fatal(FARGS
, "Rerun trajectory frame step %d time %f does not contain a box, while pbc is used", rerun_fr
.step
, rerun_fr
.time
);
794 if (max_cutoff2(ir
->ePBC
, rerun_fr
.box
) < gmx::square(fr
->rlist
))
796 gmx_fatal(FARGS
, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr
.step
, rerun_fr
.time
);
803 rerun_parallel_comm(cr
, &rerun_fr
, &bLastStep
);
806 if (ir
->ePBC
!= epbcNONE
)
808 /* Set the shift vectors.
809 * Necessary here when have a static box different from the tpr box.
811 calc_shifts(rerun_fr
.box
, fr
->shift_vec
);
815 /* Loop over MD steps or if rerunMD to end of input trajectory,
816 * or, if max_hours>0, until max_hours is reached.
818 real max_hours
= mdrunOptions
.maximumHoursToRun
;
820 /* Skip the first Nose-Hoover integration when we get the state from tpx */
821 bInitStep
= !startingFromCheckpoint
|| EI_VV(ir
->eI
);
822 bSumEkinhOld
= FALSE
;
824 bNeedRepartition
= FALSE
;
826 bool simulationsShareState
= false;
827 int nstSignalComm
= nstglobalcomm
;
829 // TODO This implementation of ensemble orientation restraints is nasty because
830 // a user can't just do multi-sim with single-sim orientation restraints.
831 bool usingEnsembleRestraints
= (fcd
->disres
.nsystems
> 1) || (ms
&& fcd
->orires
.nr
);
832 bool awhUsesMultiSim
= (ir
->bDoAwh
&& ir
->awhParams
->shareBiasMultisim
&& ms
);
834 // Replica exchange, ensemble restraints and AWH need all
835 // simulations to remain synchronized, so they need
836 // checkpoints and stop conditions to act on the same step, so
837 // the propagation of such signals must take place between
838 // simulations, not just within simulations.
839 // TODO: Make algorithm initializers set these flags.
840 simulationsShareState
= useReplicaExchange
|| usingEnsembleRestraints
|| awhUsesMultiSim
;
841 bool resetCountersIsLocal
= true;
842 signals
[eglsCHKPT
] = SimulationSignal(!simulationsShareState
);
843 signals
[eglsSTOPCOND
] = SimulationSignal(!simulationsShareState
);
844 signals
[eglsRESETCOUNTERS
] = SimulationSignal(resetCountersIsLocal
);
846 if (simulationsShareState
)
848 // Inter-simulation signal communication does not need to happen
849 // often, so we use a minimum of 200 steps to reduce overhead.
850 const int c_minimumInterSimulationSignallingInterval
= 200;
851 nstSignalComm
= ((c_minimumInterSimulationSignallingInterval
+ nstglobalcomm
- 1)/nstglobalcomm
)*nstglobalcomm
;
855 DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion
= (DOMAINDECOMP(cr
) ? DdOpenBalanceRegionBeforeForceComputation::yes
: DdOpenBalanceRegionBeforeForceComputation::no
);
856 DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion
= (DOMAINDECOMP(cr
) ? DdCloseBalanceRegionAfterForceComputation::yes
: DdCloseBalanceRegionAfterForceComputation::no
);
858 step
= ir
->init_step
;
861 // TODO extract this to new multi-simulation module
862 if (MASTER(cr
) && isMultiSim(ms
) && !useReplicaExchange
)
864 if (!multisim_int_all_are_equal(ms
, ir
->nsteps
))
866 GMX_LOG(mdlog
.warning
).appendText(
867 "Note: The number of steps is not consistent across multi simulations,\n"
868 "but we are proceeding anyway!");
870 if (!multisim_int_all_are_equal(ms
, ir
->init_step
))
872 GMX_LOG(mdlog
.warning
).appendText(
873 "Note: The initial step is not consistent across multi simulations,\n"
874 "but we are proceeding anyway!");
878 /* and stop now if we should */
879 bLastStep
= (bLastStep
|| (ir
->nsteps
>= 0 && step_rel
> ir
->nsteps
));
883 /* Determine if this is a neighbor search step */
884 bNStList
= (ir
->nstlist
> 0 && step
% ir
->nstlist
== 0);
886 if (bPMETune
&& bNStList
)
888 /* PME grid + cut-off optimization with GPUs or PME nodes */
889 pme_loadbal_do(pme_loadbal
, cr
,
890 (mdrunOptions
.verbose
&& MASTER(cr
)) ? stderr
: nullptr,
898 wallcycle_start(wcycle
, ewcSTEP
);
904 step
= rerun_fr
.step
;
905 step_rel
= step
- ir
->init_step
;
918 bLastStep
= (step_rel
== ir
->nsteps
);
919 t
= t0
+ step
*ir
->delta_t
;
922 // TODO Refactor this, so that nstfep does not need a default value of zero
923 if (ir
->efep
!= efepNO
|| ir
->bSimTemp
)
925 /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value,
926 requiring different logic. */
931 setCurrentLambdasRerun(step
, ir
->fepvals
, &rerun_fr
, lam0
, state_global
);
936 setCurrentLambdasLocal(step
, ir
->fepvals
, lam0
, state
);
938 bDoDHDL
= do_per_step(step
, ir
->fepvals
->nstdhdl
);
939 bDoFEP
= ((ir
->efep
!= efepNO
) && do_per_step(step
, nstfep
));
940 bDoExpanded
= (do_per_step(step
, ir
->expandedvals
->nstexpanded
)
941 && (ir
->bExpanded
) && (step
> 0) && (!startingFromCheckpoint
));
944 bDoReplEx
= (useReplicaExchange
&& (step
> 0) && !bLastStep
&&
945 do_per_step(step
, replExParams
.exchangeInterval
));
949 update_annealing_target_temp(ir
, t
, upd
);
952 if (bRerunMD
&& MASTER(cr
))
954 const bool constructVsites
= (vsite
&& mdrunOptions
.rerunConstructVsites
);
955 if (constructVsites
&& DOMAINDECOMP(cr
))
957 gmx_fatal(FARGS
, "Vsite recalculation with -rerun is not implemented with domain decomposition, use a single rank");
959 prepareRerunState(rerun_fr
, state_global
, constructVsites
, vsite
, top
->idef
, ir
->delta_t
, *fr
, graph
, &bRerunWarnNoV
);
962 /* Stop Center of Mass motion */
963 bStopCM
= (ir
->comm_mode
!= ecmNO
&& do_per_step(step
, ir
->nstcomm
));
967 /* for rerun MD always do Neighbour Searching */
968 bNS
= (bFirstStep
|| ir
->nstlist
!= 0);
972 /* Determine whether or not to do Neighbour Searching */
973 bNS
= (bFirstStep
|| bNStList
|| bExchanged
|| bNeedRepartition
);
976 /* < 0 means stop at next step, > 0 means stop at next NS step */
977 if ( (signals
[eglsSTOPCOND
].set
< 0) ||
978 ( (signals
[eglsSTOPCOND
].set
> 0 ) && ( bNS
|| ir
->nstlist
== 0)))
983 /* do_log triggers energy and virial calculation. Because this leads
984 * to different code paths, forces can be different. Thus for exact
985 * continuation we should avoid extra log output.
986 * Note that the || bLastStep can result in non-exact continuation
987 * beyond the last step. But we don't consider that to be an issue.
989 do_log
= do_per_step(step
, ir
->nstlog
) || (bFirstStep
&& !startingFromCheckpoint
) || bLastStep
|| bRerunMD
;
990 do_verbose
= mdrunOptions
.verbose
&&
991 (step
% mdrunOptions
.verboseStepPrintInterval
== 0 || bFirstStep
|| bLastStep
|| bRerunMD
);
993 if (bNS
&& !(bFirstStep
&& ir
->bContinuation
&& !bRerunMD
))
1001 bMasterState
= FALSE
;
1002 /* Correct the new box if it is too skewed */
1003 if (inputrecDynamicBox(ir
))
1005 if (correct_box(fplog
, step
, state
->box
, graph
))
1007 bMasterState
= TRUE
;
1010 if (DOMAINDECOMP(cr
) && bMasterState
)
1012 dd_collect_state(cr
->dd
, state
, state_global
);
1016 if (DOMAINDECOMP(cr
))
1018 /* Repartition the domain decomposition */
1019 dd_partition_system(fplog
, step
, cr
,
1020 bMasterState
, nstglobalcomm
,
1021 state_global
, top_global
, ir
,
1022 state
, &f
, mdAtoms
, top
, fr
,
1025 do_verbose
&& !bPMETunePrinting
);
1026 shouldCheckNumberOfBondedInteractions
= true;
1027 update_realloc(upd
, state
->natoms
);
1031 if (MASTER(cr
) && do_log
)
1033 print_ebin_header(fplog
, step
, t
); /* can we improve the information printed here? */
1036 if (ir
->efep
!= efepNO
)
1038 update_mdatoms(mdatoms
, state
->lambda
[efptMASS
]);
1041 if ((bRerunMD
&& rerun_fr
.bV
) || bExchanged
)
1044 /* We need the kinetic energy at minus the half step for determining
1045 * the full step kinetic energy and possibly for T-coupling.*/
1046 /* This may not be quite working correctly yet . . . . */
1047 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, mdatoms
, nrnb
, vcm
,
1048 wcycle
, enerd
, nullptr, nullptr, nullptr, nullptr, mu_tot
,
1049 constr
, &nullSignaller
, state
->box
,
1050 &totalNumberOfBondedInteractions
, &bSumEkinhOld
,
1051 CGLO_GSTAT
| CGLO_TEMPERATURE
| CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
);
1052 checkNumberOfBondedInteractions(fplog
, cr
, totalNumberOfBondedInteractions
,
1053 top_global
, top
, state
,
1054 &shouldCheckNumberOfBondedInteractions
);
1056 clear_mat(force_vir
);
1058 /* We write a checkpoint at this MD step when:
1059 * either at an NS step when we signalled through gs,
1060 * or at the last step (but not when we do not want confout),
1061 * but never at the first step or with rerun.
1063 bCPT
= (((signals
[eglsCHKPT
].set
&& (bNS
|| ir
->nstlist
== 0)) ||
1064 (bLastStep
&& mdrunOptions
.writeConfout
)) &&
1065 step
> ir
->init_step
&& !bRerunMD
);
1068 signals
[eglsCHKPT
].set
= 0;
1071 /* Determine the energy and pressure:
1072 * at nstcalcenergy steps and at energy output steps (set below).
1074 if (EI_VV(ir
->eI
) && (!bInitStep
))
1076 /* for vv, the first half of the integration actually corresponds
1077 to the previous step. bCalcEner is only required to be evaluated on the 'next' step,
1078 but the virial needs to be calculated on both the current step and the 'next' step. Future
1079 reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */
1081 /* TODO: This is probably not what we want, we will write to energy file one step after nstcalcenergy steps. */
1082 bCalcEnerStep
= do_per_step(step
- 1, ir
->nstcalcenergy
);
1083 bCalcVir
= bCalcEnerStep
||
1084 (ir
->epc
!= epcNO
&& (do_per_step(step
, ir
->nstpcouple
) || do_per_step(step
-1, ir
->nstpcouple
)));
1088 bCalcEnerStep
= do_per_step(step
, ir
->nstcalcenergy
);
1089 bCalcVir
= bCalcEnerStep
||
1090 (ir
->epc
!= epcNO
&& do_per_step(step
, ir
->nstpcouple
));
1092 bCalcEner
= bCalcEnerStep
;
1094 do_ene
= (do_per_step(step
, ir
->nstenergy
) || bLastStep
|| bRerunMD
);
1096 if (do_ene
|| do_log
|| bDoReplEx
)
1102 /* Do we need global communication ? */
1103 bGStat
= (bCalcVir
|| bCalcEner
|| bStopCM
||
1104 do_per_step(step
, nstglobalcomm
) ||
1105 (EI_VV(ir
->eI
) && inputrecNvtTrotter(ir
) && do_per_step(step
-1, nstglobalcomm
)));
1107 force_flags
= (GMX_FORCE_STATECHANGED
|
1108 ((inputrecDynamicBox(ir
) || bRerunMD
) ? GMX_FORCE_DYNAMICBOX
: 0) |
1109 GMX_FORCE_ALLFORCES
|
1110 (bCalcVir
? GMX_FORCE_VIRIAL
: 0) |
1111 (bCalcEner
? GMX_FORCE_ENERGY
: 0) |
1112 (bDoFEP
? GMX_FORCE_DHDL
: 0)
1117 /* Now is the time to relax the shells */
1118 relax_shell_flexcon(fplog
, cr
, ms
, mdrunOptions
.verbose
, step
,
1119 ir
, bNS
, force_flags
, top
,
1121 state
, f
, force_vir
, mdatoms
,
1122 nrnb
, wcycle
, graph
, groups
,
1123 shellfc
, fr
, t
, mu_tot
,
1125 ddOpenBalanceRegion
, ddCloseBalanceRegion
);
1129 /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
1130 (or the AWH update will be performed twice for one step when continuing). It would be best to
1131 call this update function from do_md_trajectory_writing but that would occur after do_force.
1132 One would have to divide the update_awh function into one function applying the AWH force
1133 and one doing the AWH bias update. The update AWH bias function could then be called after
1134 do_md_trajectory_writing (then containing update_awh_history).
1135 The checkpointing will in the future probably moved to the start of the md loop which will
1136 rid of this issue. */
1137 if (awh
&& bCPT
&& MASTER(cr
))
1139 awh
->updateHistory(state_global
->awhHistory
.get());
1142 /* The coordinates (x) are shifted (to get whole molecules)
1144 * This is parallellized as well, and does communication too.
1145 * Check comments in sim_util.c
1147 do_force(fplog
, cr
, ms
, ir
, awh
.get(),
1148 step
, nrnb
, wcycle
, top
, groups
,
1149 state
->box
, state
->x
, &state
->hist
,
1150 f
, force_vir
, mdatoms
, enerd
, fcd
,
1151 state
->lambda
, graph
,
1152 fr
, vsite
, mu_tot
, t
, ed
,
1153 (bNS
? GMX_FORCE_NS
: 0) | force_flags
,
1154 ddOpenBalanceRegion
, ddCloseBalanceRegion
);
1157 if (EI_VV(ir
->eI
) && !startingFromCheckpoint
&& !bRerunMD
)
1158 /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
1160 rvec
*vbuf
= nullptr;
1162 wallcycle_start(wcycle
, ewcUPDATE
);
1163 if (ir
->eI
== eiVV
&& bInitStep
)
1165 /* if using velocity verlet with full time step Ekin,
1166 * take the first half step only to compute the
1167 * virial for the first step. From there,
1168 * revert back to the initial coordinates
1169 * so that the input is actually the initial step.
1171 snew(vbuf
, state
->natoms
);
1172 copy_rvecn(as_rvec_array(state
->v
.data()), vbuf
, 0, state
->natoms
); /* should make this better for parallelizing? */
1176 /* this is for NHC in the Ekin(t+dt/2) version of vv */
1177 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ1
);
1180 update_coords(step
, ir
, mdatoms
, state
, f
, fcd
,
1181 ekind
, M
, upd
, etrtVELOCITY1
,
1184 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
) /* Why is rerun_fr.bV here? Unclear. */
1186 wallcycle_stop(wcycle
, ewcUPDATE
);
1187 constrain_velocities(step
, nullptr,
1191 bCalcVir
, do_log
, do_ene
);
1192 wallcycle_start(wcycle
, ewcUPDATE
);
1196 /* Need to unshift here if a do_force has been
1197 called in the previous step */
1198 unshift_self(graph
, state
->box
, as_rvec_array(state
->x
.data()));
1200 /* if VV, compute the pressure and constraints */
1201 /* For VV2, we strictly only need this if using pressure
1202 * control, but we really would like to have accurate pressures
1204 * Think about ways around this in the future?
1205 * For now, keep this choice in comments.
1207 /*bPres = (ir->eI==eiVV || inputrecNptTrotter(ir)); */
1208 /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && inputrecNptTrotter(ir)));*/
1210 bTemp
= ((ir
->eI
== eiVV
&& (!bInitStep
)) || (ir
->eI
== eiVVAK
));
1211 if (bCalcEner
&& ir
->eI
== eiVVAK
)
1213 bSumEkinhOld
= TRUE
;
1215 /* for vv, the first half of the integration actually corresponds to the previous step.
1216 So we need information from the last step in the first half of the integration */
1217 if (bGStat
|| do_per_step(step
-1, nstglobalcomm
))
1219 wallcycle_stop(wcycle
, ewcUPDATE
);
1220 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, mdatoms
, nrnb
, vcm
,
1221 wcycle
, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
1222 constr
, &nullSignaller
, state
->box
,
1223 &totalNumberOfBondedInteractions
, &bSumEkinhOld
,
1224 (bGStat
? CGLO_GSTAT
: 0)
1226 | (bTemp
? CGLO_TEMPERATURE
: 0)
1227 | (bPres
? CGLO_PRESSURE
: 0)
1228 | (bPres
? CGLO_CONSTRAINT
: 0)
1229 | (bStopCM
? CGLO_STOPCM
: 0)
1230 | (shouldCheckNumberOfBondedInteractions
? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
: 0)
1233 /* explanation of above:
1234 a) We compute Ekin at the full time step
1235 if 1) we are using the AveVel Ekin, and it's not the
1236 initial step, or 2) if we are using AveEkin, but need the full
1237 time step kinetic energy for the pressure (always true now, since we want accurate statistics).
1238 b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
1239 EkinAveVel because it's needed for the pressure */
1240 checkNumberOfBondedInteractions(fplog
, cr
, totalNumberOfBondedInteractions
,
1241 top_global
, top
, state
,
1242 &shouldCheckNumberOfBondedInteractions
);
1243 wallcycle_start(wcycle
, ewcUPDATE
);
1245 /* temperature scaling and pressure scaling to produce the extended variables at t+dt */
1250 m_add(force_vir
, shake_vir
, total_vir
); /* we need the un-dispersion corrected total vir here */
1251 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ2
);
1253 /* TODO This is only needed when we're about to write
1254 * a checkpoint, because we use it after the restart
1255 * (in a kludge?). But what should we be doing if
1256 * startingFromCheckpoint or bInitStep are true? */
1257 if (inputrecNptTrotter(ir
) || inputrecNphTrotter(ir
))
1259 copy_mat(shake_vir
, state
->svir_prev
);
1260 copy_mat(force_vir
, state
->fvir_prev
);
1262 if (inputrecNvtTrotter(ir
) && ir
->eI
== eiVV
)
1264 /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
1265 enerd
->term
[F_TEMP
] = sum_ekin(&(ir
->opts
), ekind
, nullptr, (ir
->eI
== eiVV
), FALSE
);
1266 enerd
->term
[F_EKIN
] = trace(ekind
->ekin
);
1269 else if (bExchanged
)
1271 wallcycle_stop(wcycle
, ewcUPDATE
);
1272 /* We need the kinetic energy at minus the half step for determining
1273 * the full step kinetic energy and possibly for T-coupling.*/
1274 /* This may not be quite working correctly yet . . . . */
1275 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, mdatoms
, nrnb
, vcm
,
1276 wcycle
, enerd
, nullptr, nullptr, nullptr, nullptr, mu_tot
,
1277 constr
, &nullSignaller
, state
->box
,
1278 nullptr, &bSumEkinhOld
,
1279 CGLO_GSTAT
| CGLO_TEMPERATURE
);
1280 wallcycle_start(wcycle
, ewcUPDATE
);
1283 /* if it's the initial step, we performed this first step just to get the constraint virial */
1284 if (ir
->eI
== eiVV
&& bInitStep
)
1286 copy_rvecn(vbuf
, as_rvec_array(state
->v
.data()), 0, state
->natoms
);
1289 wallcycle_stop(wcycle
, ewcUPDATE
);
1292 /* compute the conserved quantity */
1295 saved_conserved_quantity
= NPT_energy(ir
, state
, &MassQ
);
1298 last_ekin
= enerd
->term
[F_EKIN
];
1300 if ((ir
->eDispCorr
!= edispcEnerPres
) && (ir
->eDispCorr
!= edispcAllEnerPres
))
1302 saved_conserved_quantity
-= enerd
->term
[F_DISPCORR
];
1304 /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
1305 if (ir
->efep
!= efepNO
&& !bRerunMD
)
1307 sum_dhdl(enerd
, state
->lambda
, ir
->fepvals
);
1311 /* ######## END FIRST UPDATE STEP ############## */
1312 /* ######## If doing VV, we now have v(dt) ###### */
1315 /* perform extended ensemble sampling in lambda - we don't
1316 actually move to the new state before outputting
1317 statistics, but if performing simulated tempering, we
1318 do update the velocities and the tau_t. */
1320 lamnew
= ExpandedEnsembleDynamics(fplog
, ir
, enerd
, state
, &MassQ
, state
->fep_state
, state
->dfhist
, step
, as_rvec_array(state
->v
.data()), mdatoms
);
1321 /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
1324 copy_df_history(state_global
->dfhist
, state
->dfhist
);
1328 /* Now we have the energies and forces corresponding to the
1329 * coordinates at time t. We must output all of this before
1332 do_md_trajectory_writing(fplog
, cr
, nfile
, fnm
, step
, step_rel
, t
,
1333 ir
, state
, state_global
, observablesHistory
,
1335 outf
, mdebin
, ekind
, f
,
1337 bCPT
, bRerunMD
, bLastStep
,
1338 mdrunOptions
.writeConfout
,
1340 /* Check if IMD step and do IMD communication, if bIMD is TRUE. */
1341 bIMDstep
= do_IMD(ir
->bIMD
, step
, cr
, bNS
, state
->box
, as_rvec_array(state
->x
.data()), ir
, t
, wcycle
);
1343 /* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
1344 if (startingFromCheckpoint
&& (inputrecNptTrotter(ir
) || inputrecNphTrotter(ir
)))
1346 copy_mat(state
->svir_prev
, shake_vir
);
1347 copy_mat(state
->fvir_prev
, force_vir
);
1350 elapsed_time
= walltime_accounting_get_current_elapsed_time(walltime_accounting
);
1352 /* Check whether everything is still allright */
1353 if (((int)gmx_get_stop_condition() > handled_stop_condition
)
1359 int nsteps_stop
= -1;
1361 /* this just makes signals[].sig compatible with the hack
1362 of sending signals around by MPI_Reduce together with
1364 if ((gmx_get_stop_condition() == gmx_stop_cond_next_ns
) ||
1365 (mdrunOptions
.reproducible
&&
1366 gmx_get_stop_condition() == gmx_stop_cond_next
))
1368 /* We need at least two global communication steps to pass
1369 * around the signal. We stop at a pair-list creation step
1370 * to allow for exact continuation, when possible.
1372 signals
[eglsSTOPCOND
].sig
= 1;
1373 nsteps_stop
= std::max(ir
->nstlist
, 2*nstSignalComm
);
1375 else if (gmx_get_stop_condition() == gmx_stop_cond_next
)
1377 /* Stop directly after the next global communication step.
1378 * This breaks exact continuation.
1380 signals
[eglsSTOPCOND
].sig
= -1;
1381 nsteps_stop
= nstSignalComm
+ 1;
1386 "\n\nReceived the %s signal, stopping within %d steps\n\n",
1387 gmx_get_signal_name(), nsteps_stop
);
1391 "\n\nReceived the %s signal, stopping within %d steps\n\n",
1392 gmx_get_signal_name(), nsteps_stop
);
1394 handled_stop_condition
= (int)gmx_get_stop_condition();
1396 else if (MASTER(cr
) && (bNS
|| ir
->nstlist
<= 0) &&
1397 (max_hours
> 0 && elapsed_time
> max_hours
*60.0*60.0*0.99) &&
1398 signals
[eglsSTOPCOND
].sig
== 0 && signals
[eglsSTOPCOND
].set
== 0)
1400 /* Signal to terminate the run */
1401 signals
[eglsSTOPCOND
].sig
= 1;
1404 fprintf(fplog
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step
, sbuf
), max_hours
*0.99);
1406 fprintf(stderr
, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step
, sbuf
), max_hours
*0.99);
1409 if (bResetCountersHalfMaxH
&& MASTER(cr
) &&
1410 elapsed_time
> max_hours
*60.0*60.0*0.495)
1412 /* Set flag that will communicate the signal to all ranks in the simulation */
1413 signals
[eglsRESETCOUNTERS
].sig
= 1;
1416 /* In parallel we only have to check for checkpointing in steps
1417 * where we do global communication,
1418 * otherwise the other nodes don't know.
1420 const real cpt_period
= mdrunOptions
.checkpointOptions
.period
;
1421 if (MASTER(cr
) && ((bGStat
|| !PAR(cr
)) &&
1424 elapsed_time
>= nchkpt
*cpt_period
*60.0)) &&
1425 signals
[eglsCHKPT
].set
== 0)
1427 signals
[eglsCHKPT
].sig
= 1;
1430 /* ######### START SECOND UPDATE STEP ################# */
1432 /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
1435 if (ETC_ANDERSEN(ir
->etc
)) /* keep this outside of update_tcouple because of the extra info required to pass */
1437 gmx_bool bIfRandomize
;
1438 bIfRandomize
= update_randomize_velocities(ir
, step
, cr
, mdatoms
, state
, upd
, constr
);
1439 /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
1440 if (constr
&& bIfRandomize
)
1442 constrain_velocities(step
, nullptr,
1446 bCalcVir
, do_log
, do_ene
);
1449 /* Box is changed in update() when we do pressure coupling,
1450 * but we should still use the old box for energy corrections and when
1451 * writing it to the energy file, so it matches the trajectory files for
1452 * the same timestep above. Make a copy in a separate array.
1454 copy_mat(state
->box
, lastbox
);
1458 if (!bRerunMD
|| rerun_fr
.bV
|| bForceUpdate
)
1460 wallcycle_start(wcycle
, ewcUPDATE
);
1461 /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
1464 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ3
);
1465 /* We can only do Berendsen coupling after we have summed
1466 * the kinetic energy or virial. Since the happens
1467 * in global_state after update, we should only do it at
1468 * step % nstlist = 1 with bGStatEveryStep=FALSE.
1473 update_tcouple(step
, ir
, state
, ekind
, &MassQ
, mdatoms
);
1474 update_pcouple_before_coordinates(fplog
, step
, ir
, state
,
1475 parrinellorahmanMu
, M
,
1481 /* velocity half-step update */
1482 update_coords(step
, ir
, mdatoms
, state
, f
, fcd
,
1483 ekind
, M
, upd
, etrtVELOCITY2
,
1487 /* Above, initialize just copies ekinh into ekin,
1488 * it doesn't copy position (for VV),
1489 * and entire integrator for MD.
1492 if (ir
->eI
== eiVVAK
)
1494 /* We probably only need md->homenr, not state->natoms */
1495 if (state
->natoms
> cbuf_nalloc
)
1497 cbuf_nalloc
= state
->natoms
;
1498 srenew(cbuf
, cbuf_nalloc
);
1500 copy_rvecn(as_rvec_array(state
->x
.data()), cbuf
, 0, state
->natoms
);
1503 update_coords(step
, ir
, mdatoms
, state
, f
, fcd
,
1504 ekind
, M
, upd
, etrtPOSITION
, cr
, constr
);
1505 wallcycle_stop(wcycle
, ewcUPDATE
);
1507 constrain_coordinates(step
, &dvdl_constr
, state
,
1509 wcycle
, upd
, constr
,
1510 bCalcVir
, do_log
, do_ene
);
1511 update_sd_second_half(step
, &dvdl_constr
, ir
, mdatoms
, state
,
1512 cr
, nrnb
, wcycle
, upd
, constr
, do_log
, do_ene
);
1513 finish_update(ir
, mdatoms
,
1515 nrnb
, wcycle
, upd
, constr
);
1517 if (ir
->eI
== eiVVAK
)
1519 /* erase F_EKIN and F_TEMP here? */
1520 /* just compute the kinetic energy at the half step to perform a trotter step */
1521 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, mdatoms
, nrnb
, vcm
,
1522 wcycle
, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
1523 constr
, &nullSignaller
, lastbox
,
1524 nullptr, &bSumEkinhOld
,
1525 (bGStat
? CGLO_GSTAT
: 0) | CGLO_TEMPERATURE
1527 wallcycle_start(wcycle
, ewcUPDATE
);
1528 trotter_update(ir
, step
, ekind
, enerd
, state
, total_vir
, mdatoms
, &MassQ
, trotter_seq
, ettTSEQ4
);
1529 /* now we know the scaling, we can compute the positions again again */
1530 copy_rvecn(cbuf
, as_rvec_array(state
->x
.data()), 0, state
->natoms
);
1532 update_coords(step
, ir
, mdatoms
, state
, f
, fcd
,
1533 ekind
, M
, upd
, etrtPOSITION
, cr
, constr
);
1534 wallcycle_stop(wcycle
, ewcUPDATE
);
1536 /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
1537 /* are the small terms in the shake_vir here due
1538 * to numerical errors, or are they important
1539 * physically? I'm thinking they are just errors, but not completely sure.
1540 * For now, will call without actually constraining, constr=NULL*/
1541 finish_update(ir
, mdatoms
,
1543 nrnb
, wcycle
, upd
, nullptr);
1547 /* this factor or 2 correction is necessary
1548 because half of the constraint force is removed
1549 in the vv step, so we have to double it. See
1550 the Redmine issue #1255. It is not yet clear
1551 if the factor of 2 is exact, or just a very
1552 good approximation, and this will be
1553 investigated. The next step is to see if this
1554 can be done adding a dhdl contribution from the
1555 rattle step, but this is somewhat more
1556 complicated with the current code. Will be
1557 investigated, hopefully for 4.6.3. However,
1558 this current solution is much better than
1559 having it completely wrong.
1561 enerd
->term
[F_DVDL_CONSTR
] += 2*dvdl_constr
;
1565 enerd
->term
[F_DVDL_CONSTR
] += dvdl_constr
;
1570 /* Need to unshift here */
1571 unshift_self(graph
, state
->box
, as_rvec_array(state
->x
.data()));
1574 if (vsite
!= nullptr)
1576 wallcycle_start(wcycle
, ewcVSITECONSTR
);
1577 if (graph
!= nullptr)
1579 shift_self(graph
, state
->box
, as_rvec_array(state
->x
.data()));
1581 construct_vsites(vsite
, as_rvec_array(state
->x
.data()), ir
->delta_t
, as_rvec_array(state
->v
.data()),
1582 top
->idef
.iparams
, top
->idef
.il
,
1583 fr
->ePBC
, fr
->bMolPBC
, cr
, state
->box
);
1585 if (graph
!= nullptr)
1587 unshift_self(graph
, state
->box
, as_rvec_array(state
->x
.data()));
1589 wallcycle_stop(wcycle
, ewcVSITECONSTR
);
1592 /* ############## IF NOT VV, Calculate globals HERE ############ */
1593 /* With Leap-Frog we can skip compute_globals at
1594 * non-communication steps, but we need to calculate
1595 * the kinetic energy one step before communication.
1598 // Organize to do inter-simulation signalling on steps if
1599 // and when algorithms require it.
1600 bool doInterSimSignal
= (simulationsShareState
&& do_per_step(step
, nstSignalComm
));
1602 if (bGStat
|| (!EI_VV(ir
->eI
) && do_per_step(step
+1, nstglobalcomm
)) || doInterSimSignal
)
1604 // Since we're already communicating at this step, we
1605 // can propagate intra-simulation signals. Note that
1606 // check_nstglobalcomm has the responsibility for
1607 // choosing the value of nstglobalcomm that is one way
1608 // bGStat becomes true, so we can't get into a
1609 // situation where e.g. checkpointing can't be
1611 bool doIntraSimSignal
= true;
1612 SimulationSignaller
signaller(&signals
, cr
, ms
, doInterSimSignal
, doIntraSimSignal
);
1614 compute_globals(fplog
, gstat
, cr
, ir
, fr
, ekind
, state
, mdatoms
, nrnb
, vcm
,
1615 wcycle
, enerd
, force_vir
, shake_vir
, total_vir
, pres
, mu_tot
,
1618 &totalNumberOfBondedInteractions
, &bSumEkinhOld
,
1619 (bGStat
? CGLO_GSTAT
: 0)
1620 | (!EI_VV(ir
->eI
) || bRerunMD
? CGLO_ENERGY
: 0)
1621 | (!EI_VV(ir
->eI
) && bStopCM
? CGLO_STOPCM
: 0)
1622 | (!EI_VV(ir
->eI
) ? CGLO_TEMPERATURE
: 0)
1623 | (!EI_VV(ir
->eI
) || bRerunMD
? CGLO_PRESSURE
: 0)
1625 | (shouldCheckNumberOfBondedInteractions
? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
: 0)
1627 checkNumberOfBondedInteractions(fplog
, cr
, totalNumberOfBondedInteractions
,
1628 top_global
, top
, state
,
1629 &shouldCheckNumberOfBondedInteractions
);
1633 /* ############# END CALC EKIN AND PRESSURE ################# */
1635 /* Note: this is OK, but there are some numerical precision issues with using the convergence of
1636 the virial that should probably be addressed eventually. state->veta has better properies,
1637 but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
1638 generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
1640 if (ir
->efep
!= efepNO
&& (!EI_VV(ir
->eI
) || bRerunMD
))
1642 /* Sum up the foreign energy and dhdl terms for md and sd.
1643 Currently done every step so that dhdl is correct in the .edr */
1644 sum_dhdl(enerd
, state
->lambda
, ir
->fepvals
);
1647 update_pcouple_after_coordinates(fplog
, step
, ir
, mdatoms
,
1648 pres
, force_vir
, shake_vir
,
1652 /* ################# END UPDATE STEP 2 ################# */
1653 /* #### We now have r(t+dt) and v(t+dt/2) ############# */
1655 /* The coordinates (x) were unshifted in update */
1658 /* We will not sum ekinh_old,
1659 * so signal that we still have to do it.
1661 bSumEkinhOld
= TRUE
;
1666 /* ######### BEGIN PREPARING EDR OUTPUT ########### */
1668 /* use the directly determined last velocity, not actually the averaged half steps */
1669 if (bTrotter
&& ir
->eI
== eiVV
)
1671 enerd
->term
[F_EKIN
] = last_ekin
;
1673 enerd
->term
[F_ETOT
] = enerd
->term
[F_EPOT
] + enerd
->term
[F_EKIN
];
1675 if (integratorHasConservedEnergyQuantity(ir
))
1679 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + saved_conserved_quantity
;
1683 enerd
->term
[F_ECONSERVED
] = enerd
->term
[F_ETOT
] + NPT_energy(ir
, state
, &MassQ
);
1686 /* ######### END PREPARING EDR OUTPUT ########### */
1692 if (fplog
&& do_log
&& bDoExpanded
)
1694 /* only needed if doing expanded ensemble */
1695 PrintFreeEnergyInfoToFile(fplog
, ir
->fepvals
, ir
->expandedvals
, ir
->bSimTemp
? ir
->simtempvals
: nullptr,
1696 state_global
->dfhist
, state
->fep_state
, ir
->nstlog
, step
);
1700 upd_mdebin(mdebin
, bDoDHDL
, bCalcEnerStep
,
1701 t
, mdatoms
->tmass
, enerd
, state
,
1702 ir
->fepvals
, ir
->expandedvals
, lastbox
,
1703 shake_vir
, force_vir
, total_vir
, pres
,
1704 ekind
, mu_tot
, constr
);
1708 upd_mdebin_step(mdebin
);
1711 gmx_bool do_dr
= do_per_step(step
, ir
->nstdisreout
);
1712 gmx_bool do_or
= do_per_step(step
, ir
->nstorireout
);
1714 print_ebin(mdoutf_get_fp_ene(outf
), do_ene
, do_dr
, do_or
, do_log
? fplog
: nullptr,
1716 eprNORMAL
, mdebin
, fcd
, groups
, &(ir
->opts
), awh
.get());
1720 pull_print_output(ir
->pull_work
, step
, t
);
1723 if (do_per_step(step
, ir
->nstlog
))
1725 if (fflush(fplog
) != 0)
1727 gmx_fatal(FARGS
, "Cannot flush logfile - maybe you are out of disk space?");
1733 /* Have to do this part _after_ outputting the logfile and the edr file */
1734 /* Gets written into the state at the beginning of next loop*/
1735 state
->fep_state
= lamnew
;
1737 /* Print the remaining wall clock time for the run */
1738 if (isMasterSimMasterRank(ms
, cr
) &&
1739 (do_verbose
|| gmx_got_usr_signal()) &&
1744 fprintf(stderr
, "\n");
1746 print_time(stderr
, walltime_accounting
, step
, ir
, cr
);
1749 /* Ion/water position swapping.
1750 * Not done in last step since trajectory writing happens before this call
1751 * in the MD loop and exchanges would be lost anyway. */
1752 bNeedRepartition
= FALSE
;
1753 if ((ir
->eSwapCoords
!= eswapNO
) && (step
> 0) && !bLastStep
&&
1754 do_per_step(step
, ir
->swap
->nstswap
))
1756 bNeedRepartition
= do_swapcoords(cr
, step
, t
, ir
, wcycle
,
1757 bRerunMD
? rerun_fr
.x
: as_rvec_array(state
->x
.data()),
1758 bRerunMD
? rerun_fr
.box
: state
->box
,
1759 MASTER(cr
) && mdrunOptions
.verbose
,
1762 if (bNeedRepartition
&& DOMAINDECOMP(cr
))
1764 dd_collect_state(cr
->dd
, state
, state_global
);
1768 /* Replica exchange */
1772 bExchanged
= replica_exchange(fplog
, cr
, ms
, repl_ex
,
1773 state_global
, enerd
,
1777 if ( (bExchanged
|| bNeedRepartition
) && DOMAINDECOMP(cr
) )
1779 dd_partition_system(fplog
, step
, cr
, TRUE
, 1,
1780 state_global
, top_global
, ir
,
1781 state
, &f
, mdAtoms
, top
, fr
,
1783 nrnb
, wcycle
, FALSE
);
1784 shouldCheckNumberOfBondedInteractions
= true;
1785 update_realloc(upd
, state
->natoms
);
1790 startingFromCheckpoint
= false;
1792 /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
1793 /* With all integrators, except VV, we need to retain the pressure
1794 * at the current step for coupling at the next step.
1796 if ((state
->flags
& (1<<estPRES_PREV
)) &&
1798 (ir
->nstpcouple
> 0 && step
% ir
->nstpcouple
== 0)))
1800 /* Store the pressure in t_state for pressure coupling
1801 * at the next MD step.
1803 copy_mat(pres
, state
->pres_prev
);
1806 /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
1808 if ( (membed
!= nullptr) && (!bLastStep
) )
1810 rescale_membed(step_rel
, membed
, as_rvec_array(state_global
->x
.data()));
1817 /* read next frame from input trajectory */
1818 bLastStep
= !read_next_frame(oenv
, status
, &rerun_fr
);
1823 rerun_parallel_comm(cr
, &rerun_fr
, &bLastStep
);
1827 cycles
= wallcycle_stop(wcycle
, ewcSTEP
);
1828 if (DOMAINDECOMP(cr
) && wcycle
)
1830 dd_cycles_add(cr
->dd
, cycles
, ddCyclStep
);
1833 if (!bRerunMD
|| !rerun_fr
.bStep
)
1835 /* increase the MD step number */
1840 /* TODO make a counter-reset module */
1841 /* If it is time to reset counters, set a flag that remains
1842 true until counters actually get reset */
1843 if (step_rel
== wcycle_get_reset_counters(wcycle
) ||
1844 signals
[eglsRESETCOUNTERS
].set
!= 0)
1846 if (pme_loadbal_is_active(pme_loadbal
))
1848 /* Do not permit counter reset while PME load
1849 * balancing is active. The only purpose for resetting
1850 * counters is to measure reliable performance data,
1851 * and that can't be done before balancing
1854 * TODO consider fixing this by delaying the reset
1855 * until after load balancing completes,
1856 * e.g. https://gerrit.gromacs.org/#/c/4964/2 */
1857 gmx_fatal(FARGS
, "PME tuning was still active when attempting to "
1858 "reset mdrun counters at step %" GMX_PRId64
". Try "
1859 "resetting counters later in the run, e.g. with gmx "
1860 "mdrun -resetstep.", step
);
1862 reset_all_counters(fplog
, mdlog
, cr
, step
, &step_rel
, ir
, wcycle
, nrnb
, walltime_accounting
,
1863 use_GPU(fr
->nbv
) ? fr
->nbv
: nullptr, fr
->pmedata
);
1864 wcycle_set_reset_counters(wcycle
, -1);
1865 if (!thisRankHasDuty(cr
, DUTY_PME
))
1867 /* Tell our PME node to reset its counters */
1868 gmx_pme_send_resetcounters(cr
, step
);
1870 /* Correct max_hours for the elapsed time */
1871 max_hours
-= elapsed_time
/(60.0*60.0);
1872 /* If mdrun -maxh -resethway was active, it can only trigger once */
1873 bResetCountersHalfMaxH
= FALSE
; /* TODO move this to where signals[eglsRESETCOUNTERS].sig is set */
1874 /* Reset can only happen once, so clear the triggering flag. */
1875 signals
[eglsRESETCOUNTERS
].set
= 0;
1878 /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
1879 IMD_prep_energies_send_positions(ir
->bIMD
&& MASTER(cr
), bIMDstep
, ir
->imd
, enerd
, step
, bCalcEner
, wcycle
);
1882 /* End of main MD loop */
1884 /* Closing TNG files can include compressing data. Therefore it is good to do that
1885 * before stopping the time measurements. */
1886 mdoutf_tng_close(outf
);
1888 /* Stop measuring walltime */
1889 walltime_accounting_end(walltime_accounting
);
1891 if (bRerunMD
&& MASTER(cr
))
1896 if (!thisRankHasDuty(cr
, DUTY_PME
))
1898 /* Tell the PME only node to finish */
1899 gmx_pme_send_finish(cr
);
1904 if (ir
->nstcalcenergy
> 0 && !bRerunMD
)
1906 print_ebin(mdoutf_get_fp_ene(outf
), FALSE
, FALSE
, FALSE
, fplog
, step
, t
,
1907 eprAVER
, mdebin
, fcd
, groups
, &(ir
->opts
), awh
.get());
1910 done_mdebin(mdebin
);
1915 pme_loadbal_done(pme_loadbal
, fplog
, mdlog
, use_GPU(fr
->nbv
));
1918 done_shellfc(fplog
, shellfc
, step_rel
);
1920 if (useReplicaExchange
&& MASTER(cr
))
1922 print_replica_exchange_statistics(fplog
, repl_ex
);
1925 // Clean up swapcoords
1926 if (ir
->eSwapCoords
!= eswapNO
)
1928 finish_swapcoords(ir
->swap
);
1931 /* Do essential dynamics cleanup if needed. Close .edo file */
1934 /* IMD cleanup, if bIMD is TRUE. */
1935 IMD_finalize(ir
->bIMD
, ir
->imd
);
1937 walltime_accounting_set_nsteps_done(walltime_accounting
, step_rel
);
1938 if (step_rel
>= wcycle_get_reset_counters(wcycle
) &&
1939 signals
[eglsRESETCOUNTERS
].set
== 0 &&
1940 !bResetCountersHalfMaxH
)
1942 walltime_accounting_set_valid_finish(walltime_accounting
);
1945 destroy_enerdata(enerd
);