#include <stdlib.h>
#include <algorithm>
+#include <memory>
#include "thread_mpi/threads.h"
#include "gromacs/domdec/domdec.h"
#include "gromacs/domdec/domdec_network.h"
#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/essentialdynamics/edsam.h"
#include "gromacs/ewald/pme.h"
#include "gromacs/ewald/pme-load-balancing.h"
#include "gromacs/fileio/trxio.h"
-#include "gromacs/gmxlib/md_logging.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/gmxlib/nrnb.h"
#include "gromacs/gpu_utils/gpu_utils.h"
#include "gromacs/mdlib/mdebin.h"
#include "gromacs/mdlib/mdoutf.h"
#include "gromacs/mdlib/mdrun.h"
+#include "gromacs/mdlib/mdsetup.h"
#include "gromacs/mdlib/nb_verlet.h"
#include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
#include "gromacs/mdlib/ns.h"
#include "gromacs/mdtypes/interaction_const.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/mdtypes/observableshistory.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/logger.h"
#include "gromacs/utility/real.h"
#include "gromacs/utility/smalloc.h"
}
}
-static void reset_all_counters(FILE *fplog, t_commrec *cr,
+static void reset_all_counters(FILE *fplog, const gmx::MDLogger &mdlog, t_commrec *cr,
gmx_int64_t step,
gmx_int64_t *step_rel, t_inputrec *ir,
gmx_wallcycle_t wcycle, t_nrnb *nrnb,
char sbuf[STEPSTRSIZE];
/* Reset all the counters related to performance over the run */
- md_print_warn(cr, fplog, "step %s: resetting all time and cycle counters\n",
- gmx_step_str(step, sbuf));
+ GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
+ "step %s: resetting all time and cycle counters",
+ gmx_step_str(step, sbuf));
if (use_GPU(nbv))
{
}
/*! \libinternal
- \copydoc integrator_t (FILE *fplog, t_commrec *cr,
+ \copydoc integrator_t (FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
int nfile, const t_filenm fnm[],
const gmx_output_env_t *oenv, gmx_bool bVerbose,
int nstglobalcomm,
gmx_vsite_t *vsite, gmx_constr_t constr,
int stepout,
+ gmx::IMDOutputProvider *outputProvider,
t_inputrec *inputrec,
gmx_mtop_t *top_global, t_fcdata *fcd,
t_state *state_global,
t_mdatoms *mdatoms,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
real cpt_period, real max_hours,
unsigned long Flags,
gmx_walltime_accounting_t walltime_accounting)
*/
-double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
+double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
+ int nfile, const t_filenm fnm[],
const gmx_output_env_t *oenv, gmx_bool bVerbose,
int nstglobalcomm,
gmx_vsite_t *vsite, gmx_constr_t constr,
- int stepout, t_inputrec *ir,
+ int stepout, gmx::IMDOutputProvider *outputProvider,
+ t_inputrec *ir,
gmx_mtop_t *top_global,
t_fcdata *fcd,
t_state *state_global,
+ ObservablesHistory *observablesHistory,
t_mdatoms *mdatoms,
t_nrnb *nrnb, gmx_wallcycle_t wcycle,
- gmx_edsam_t ed, t_forcerec *fr,
- int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+ t_forcerec *fr,
+ const ReplicaExchangeParameters &replExParams,
gmx_membed_t *membed,
real cpt_period, real max_hours,
int imdport,
unsigned long Flags,
gmx_walltime_accounting_t walltime_accounting)
{
- gmx_mdoutf_t outf = NULL;
+ gmx_mdoutf_t outf = nullptr;
gmx_int64_t step, step_rel;
double elapsed_time;
double t, t0, lam0[efptNR];
t_trxstatus *status;
rvec mu_tot;
t_vcm *vcm;
- matrix pcoupl_mu, M;
+ matrix parrinellorahmanMu, M;
t_trxframe rerun_fr;
- gmx_repl_ex_t repl_ex = NULL;
+ gmx_repl_ex_t repl_ex = nullptr;
int nchkpt = 1;
gmx_localtop_t *top;
- t_mdebin *mdebin = NULL;
- t_state *state = NULL;
+ t_mdebin *mdebin = nullptr;
gmx_enerdata_t *enerd;
- rvec *f = NULL;
+ PaddedRVecVector f {};
gmx_global_stat_t gstat;
- gmx_update_t *upd = NULL;
- t_graph *graph = NULL;
+ gmx_update_t *upd = nullptr;
+ t_graph *graph = nullptr;
gmx_groups_t *groups;
gmx_ekindata_t *ekind;
gmx_shellfc_t *shellfc;
gmx_bool bResetCountersHalfMaxH = FALSE;
gmx_bool bTemp, bPres, bTrotter;
real dvdl_constr;
- rvec *cbuf = NULL;
+ rvec *cbuf = nullptr;
int cbuf_nalloc = 0;
matrix lastbox;
int lamnew = 0;
/* PME load balancing data for GPU kernels */
- pme_load_balancing_t *pme_loadbal = NULL;
+ pme_load_balancing_t *pme_loadbal = nullptr;
gmx_bool bPMETune = FALSE;
gmx_bool bPMETunePrinting = FALSE;
nstglobalcomm = 1;
}
- nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir);
+ nstglobalcomm = check_nstglobalcomm(mdlog, nstglobalcomm, ir);
bGStatEveryStep = (nstglobalcomm == 1);
if (bRerunMD)
}
groups = &top_global->groups;
+ gmx_edsam *ed = nullptr;
+ if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
+ {
+ /* Initialize essential dynamics sampling */
+ ed = init_edsam(opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
+ top_global, ir, cr, constr,
+ as_rvec_array(state_global->x.data()),
+ state_global->box, observablesHistory,
+ oenv, Flags & MD_APPENDFILES);
+ }
+
if (ir->eSwapCoords != eswapNO)
{
/* Initialize ion swapping code */
init_swapcoords(fplog, bVerbose, ir, opt2fn_master("-swap", nfile, fnm, cr),
- top_global, state_global->x, state_global->box, &state_global->swapstate, cr, oenv, Flags);
+ top_global, as_rvec_array(state_global->x.data()), state_global->box, observablesHistory, cr, oenv, Flags);
}
/* Initial values */
- init_md(fplog, cr, ir, oenv, &t, &t0, state_global->lambda,
+ init_md(fplog, cr, outputProvider, ir, oenv, &t, &t0, state_global->lambda,
&(state_global->fep_state), lam0,
nrnb, top_global, &upd,
nfile, fnm, &outf, &mdebin,
snew(enerd, 1);
init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
enerd);
- if (DOMAINDECOMP(cr))
- {
- f = NULL;
- }
- else
- {
- snew(f, top_global->natoms);
- }
/* Kinetic energy data */
snew(ekind, 1);
}
}
+ std::unique_ptr<t_state> stateInstance;
+ t_state * state;
+
if (DOMAINDECOMP(cr))
{
top = dd_init_local_top(top_global);
- snew(state, 1);
+ stateInstance = std::unique_ptr<t_state>(new t_state);
+ state = stateInstance.get();
dd_init_local_state(cr->dd, state_global, state);
}
else
{
- top = gmx_mtop_generate_local_top(top_global, ir->efep != efepNO);
-
- state = serial_init_local_state(state_global);
-
- atoms2md(top_global, ir, 0, NULL, top_global->natoms, mdatoms);
-
- if (vsite)
- {
- set_vsite_top(vsite, top, mdatoms, cr);
- }
-
- if (ir->ePBC != epbcNONE && !fr->bMolPBC)
- {
- graph = mk_graph(fplog, &(top->idef), 0, top_global->natoms, FALSE, FALSE);
- }
-
- if (shellfc)
- {
- make_local_shells(cr, mdatoms, shellfc);
- }
+ state_change_natoms(state_global, state_global->natoms);
+ /* We need to allocate one element extra, since we might use
+ * (unaligned) 4-wide SIMD loads to access rvec entries.
+ */
+ f.resize(state_global->natoms + 1);
+ /* Copy the pointer to the global state */
+ state = state_global;
- setup_bonded_threading(fr, &top->idef);
+ snew(top, 1);
+ mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
+ &graph, mdatoms, vsite, shellfc);
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
/* Set up interactive MD (IMD) */
- init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, state_global->x,
+ init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, as_rvec_array(state_global->x.data()),
nfile, fnm, oenv, imdport, Flags);
if (DOMAINDECOMP(cr))
state_global, top_global, ir,
state, &f, mdatoms, top, fr,
vsite, constr,
- nrnb, NULL, FALSE);
+ nrnb, nullptr, FALSE);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
update_mdatoms(mdatoms, state->lambda[efptMASS]);
if (ir->bExpanded)
{
- init_expanded_ensemble(startingFromCheckpoint, ir, &state->dfhist);
+ init_expanded_ensemble(startingFromCheckpoint, ir, state->dfhist);
}
if (MASTER(cr))
/* Update mdebin with energy history if appending to output files */
if (Flags & MD_APPENDFILES)
{
- restore_energyhistory_from_state(mdebin, state_global->enerhist);
+ restore_energyhistory_from_state(mdebin, observablesHistory->energyHistory.get());
}
- else
+ else if (observablesHistory->energyHistory.get() != nullptr)
{
- /* We might have read an energy history from checkpoint,
- * free the allocated memory and reset the counts.
+ /* We might have read an energy history from checkpoint.
+ * As we are not appending, we want to restart the statistics.
+ * Free the allocated memory and reset the counts.
*/
- done_energyhistory(state_global->enerhist);
- init_energyhistory(state_global->enerhist);
+ observablesHistory->energyHistory = {};
}
}
+ if (observablesHistory->energyHistory.get() == nullptr)
+ {
+ observablesHistory->energyHistory = std::unique_ptr<energyhistory_t>(new energyhistory_t {});
+ }
/* Set the initial energy history in state by updating once */
- update_energyhistory(state_global->enerhist, mdebin);
+ update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
}
/* Initialize constraints */
set_constraints(constr, top, ir, mdatoms, cr);
}
- if (repl_ex_nst > 0 && MASTER(cr))
+ const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
+ if (useReplicaExchange && MASTER(cr))
{
repl_ex = init_replica_exchange(fplog, cr->ms, state_global, ir,
- repl_ex_nst, repl_ex_nex, repl_ex_seed);
+ replExParams);
}
/* PME tuning is only supported in the Verlet scheme, with PME for
!(Flags & MD_REPRODUCIBLE) && ir->cutoff_scheme != ecutsGROUP);
if (bPMETune)
{
- pme_loadbal_init(&pme_loadbal, cr, fplog, ir, state->box,
- fr->ic, fr->pmedata, use_GPU(fr->nbv),
+ pme_loadbal_init(&pme_loadbal, cr, mdlog, ir, state->box,
+ fr->ic, fr->nbv->listParams.get(), fr->pmedata, use_GPU(fr->nbv),
&bPMETunePrinting);
}
if (!ir->bContinuation && !bRerunMD)
{
- if (mdatoms->cFREEZE && (state->flags & (1<<estV)))
+ if (state->flags & (1 << estV))
{
- /* Set the velocities of frozen particles to zero */
+ /* Set the velocities of vsites, shells and frozen atoms to zero */
for (i = 0; i < mdatoms->homenr; i++)
{
- for (m = 0; m < DIM; m++)
+ if (mdatoms->ptype[i] == eptVSite ||
+ mdatoms->ptype[i] == eptShell)
+ {
+ clear_rvec(state->v[i]);
+ }
+ else if (mdatoms->cFREEZE)
{
- if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+ for (m = 0; m < DIM; m++)
{
- state->v[i][m] = 0;
+ if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+ {
+ state->v[i][m] = 0;
+ }
}
}
}
if (vsite)
{
/* Construct the virtual sites for the initial configuration */
- construct_vsites(vsite, state->x, ir->delta_t, NULL,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, nullptr,
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
}
{
nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
}
- if (repl_ex_nst > 0)
+ if (useReplicaExchange)
{
- nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep);
+ nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
}
}
bSumEkinhOld = FALSE;
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
- NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, state->box,
&totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags
| (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
perhaps loses some logic?*/
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
- NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, state->box,
- NULL, &bSumEkinhOld,
+ nullptr, &bSumEkinhOld,
cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
}
// checkpoints and stop conditions to act on the same step, so
// the propagation of such signals must take place between
// simulations, not just within simulations.
- bool checkpointIsLocal = (repl_ex_nst <= 0) && !bUsingEnsembleRestraints;
- bool stopConditionIsLocal = (repl_ex_nst <= 0) && !bUsingEnsembleRestraints;
+ bool checkpointIsLocal = !useReplicaExchange && !bUsingEnsembleRestraints;
+ bool stopConditionIsLocal = !useReplicaExchange && !bUsingEnsembleRestraints;
bool resetCountersIsLocal = true;
signals[eglsCHKPT] = SimulationSignal(checkpointIsLocal);
signals[eglsSTOPCOND] = SimulationSignal(stopConditionIsLocal);
signals[eglsRESETCOUNTERS] = SimulationSignal(resetCountersIsLocal);
}
+ DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
+ DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion = (DOMAINDECOMP(cr) ? DdCloseBalanceRegionAfterForceComputation::yes : DdCloseBalanceRegionAfterForceComputation::no);
+
step = ir->init_step;
step_rel = 0;
// TODO extract this to new multi-simulation module
- if (MASTER(cr) && MULTISIM(cr) && (repl_ex_nst <= 0 ))
+ if (MASTER(cr) && MULTISIM(cr) && !useReplicaExchange)
{
if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
{
- md_print_info(cr, fplog,
- "Note: The number of steps is not consistent across multi simulations,\n"
- "but we are proceeding anyway!\n");
+ GMX_LOG(mdlog.warning).appendText(
+ "Note: The number of steps is not consistent across multi simulations,\n"
+ "but we are proceeding anyway!");
}
if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
{
- md_print_info(cr, fplog,
- "Note: The initial step is not consistent across multi simulations,\n"
- "but we are proceeding anyway!\n");
+ GMX_LOG(mdlog.warning).appendText(
+ "Note: The initial step is not consistent across multi simulations,\n"
+ "but we are proceeding anyway!");
}
}
{
/* PME grid + cut-off optimization with GPUs or PME nodes */
pme_loadbal_do(pme_loadbal, cr,
- (bVerbose && MASTER(cr)) ? stderr : NULL,
- fplog,
+ (bVerbose && MASTER(cr)) ? stderr : nullptr,
+ fplog, mdlog,
ir, fr, state,
wcycle,
step, step_rel,
&& (ir->bExpanded) && (step > 0) && (!startingFromCheckpoint));
}
- bDoReplEx = ((repl_ex_nst > 0) && (step > 0) && !bLastStep &&
- do_per_step(step, repl_ex_nst));
+ bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
+ do_per_step(step, replExParams.exchangeInterval));
if (bSimAnn)
{
/* Following is necessary because the graph may get out of sync
* with the coordinates if we only have every N'th coordinate set
*/
- mk_mshift(fplog, graph, fr->ePBC, state->box, state->x);
- shift_self(graph, state->box, state->x);
+ mk_mshift(fplog, graph, fr->ePBC, state->box, as_rvec_array(state->x.data()));
+ shift_self(graph, state->box, as_rvec_array(state->x.data()));
}
- construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
if (graph)
{
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
}
}
nrnb, wcycle,
do_verbose && !bPMETunePrinting);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
}
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
- wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+ wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
constr, &nullSignaller, state->box,
&totalNumberOfBondedInteractions, &bSumEkinhOld,
CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
relax_shell_flexcon(fplog, cr, bVerbose, step,
ir, bNS, force_flags, top,
constr, enerd, fcd,
- state, f, force_vir, mdatoms,
+ state, &f, force_vir, mdatoms,
nrnb, wcycle, graph, groups,
shellfc, fr, bBornRadii, t, mu_tot,
- vsite, mdoutf_get_fp_field(outf));
+ vsite,
+ ddOpenBalanceRegion, ddCloseBalanceRegion);
}
else
{
* Check comments in sim_util.c
*/
do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
- state->box, state->x, &state->hist,
- f, force_vir, mdatoms, enerd, fcd,
+ state->box, &state->x, &state->hist,
+ &f, force_vir, mdatoms, enerd, fcd,
state->lambda, graph,
- fr, vsite, mu_tot, t, mdoutf_get_fp_field(outf), ed, bBornRadii,
- (bNS ? GMX_FORCE_NS : 0) | force_flags);
+ fr, vsite, mu_tot, t, ed, bBornRadii,
+ (bNS ? GMX_FORCE_NS : 0) | force_flags,
+ ddOpenBalanceRegion, ddCloseBalanceRegion);
}
if (EI_VV(ir->eI) && !startingFromCheckpoint && !bRerunMD)
/* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
{
- rvec *vbuf = NULL;
+ rvec *vbuf = nullptr;
wallcycle_start(wcycle, ewcUPDATE);
if (ir->eI == eiVV && bInitStep)
* so that the input is actually the initial step.
*/
snew(vbuf, state->natoms);
- copy_rvecn(state->v, vbuf, 0, state->natoms); /* should make this better for parallelizing? */
+ copy_rvecn(as_rvec_array(state->v.data()), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
}
else
{
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
}
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtVELOCITY1,
cr, constr);
if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */
{
wallcycle_stop(wcycle, ewcUPDATE);
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ update_constraints(fplog, step, nullptr, ir, mdatoms,
+ state, fr->bMolPBC, graph, &f,
&top->idef, shake_vir,
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir);
{
/* Need to unshift here if a do_force has been
called in the previous step */
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
/* if VV, compute the pressure and constraints */
/* For VV2, we strictly only need this if using pressure
if (inputrecNvtTrotter(ir) && ir->eI == eiVV)
{
/* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */
- enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL, (ir->eI == eiVV), FALSE);
+ enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, nullptr, (ir->eI == eiVV), FALSE);
enerd->term[F_EKIN] = trace(ekind->ekin);
}
}
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
- wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
+ wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
constr, &nullSignaller, state->box,
- NULL, &bSumEkinhOld,
+ nullptr, &bSumEkinhOld,
CGLO_GSTAT | CGLO_TEMPERATURE);
wallcycle_start(wcycle, ewcUPDATE);
}
/* if it's the initial step, we performed this first step just to get the constraint virial */
if (ir->eI == eiVV && bInitStep)
{
- copy_rvecn(vbuf, state->v, 0, state->natoms);
+ copy_rvecn(vbuf, as_rvec_array(state->v.data()), 0, state->natoms);
sfree(vbuf);
}
wallcycle_stop(wcycle, ewcUPDATE);
/* compute the conserved quantity */
if (EI_VV(ir->eI))
{
- saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ);
+ saved_conserved_quantity = NPT_energy(ir, state, &MassQ);
if (ir->eI == eiVV)
{
last_ekin = enerd->term[F_EKIN];
statistics, but if performing simulated tempering, we
do update the velocities and the tau_t. */
- lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, &state->dfhist, step, state->v, mdatoms);
+ lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, state->dfhist, step, as_rvec_array(state->v.data()), mdatoms);
/* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
- copy_df_history(&state_global->dfhist, &state->dfhist);
+ copy_df_history(state_global->dfhist, state->dfhist);
}
/* Now we have the energies and forces corresponding to the
* the update.
*/
do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
- ir, state, state_global, top_global, fr,
- outf, mdebin, ekind, f,
+ ir, state, state_global, observablesHistory,
+ top_global, fr,
+ outf, mdebin, ekind, &f,
&nchkpt,
bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT),
bSumEkinhOld);
/* Check if IMD step and do IMD communication, if bIMD is TRUE. */
- bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle);
+ bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, as_rvec_array(state->x.data()), ir, t, wcycle);
/* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
if (startingFromCheckpoint && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
/* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
if (constr && bIfRandomize)
{
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ update_constraints(fplog, step, nullptr, ir, mdatoms,
+ state, fr->bMolPBC, graph, &f,
&top->idef, tmp_vir,
cr, nrnb, wcycle, upd, constr,
TRUE, bCalcVir);
else
{
update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
- update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep);
+ update_pcouple_before_coordinates(fplog, step, ir, state,
+ parrinellorahmanMu, M,
+ bInitStep);
}
if (EI_VV(ir->eI))
{
/* velocity half-step update */
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtVELOCITY2,
cr, constr);
}
cbuf_nalloc = state->natoms;
srenew(cbuf, cbuf_nalloc);
}
- copy_rvecn(state->x, cbuf, 0, state->natoms);
+ copy_rvecn(as_rvec_array(state->x.data()), cbuf, 0, state->natoms);
}
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
- fr->bMolPBC, graph, f,
+ fr->bMolPBC, graph, &f,
&top->idef, shake_vir,
cr, nrnb, wcycle, upd, constr,
FALSE, bCalcVir);
compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
constr, &nullSignaller, lastbox,
- NULL, &bSumEkinhOld,
+ nullptr, &bSumEkinhOld,
(bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
);
wallcycle_start(wcycle, ewcUPDATE);
trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4);
/* now we know the scaling, we can compute the positions again again */
- copy_rvecn(cbuf, state->x, 0, state->natoms);
+ copy_rvecn(cbuf, as_rvec_array(state->x.data()), 0, state->natoms);
- update_coords(fplog, step, ir, mdatoms, state, f, fcd,
+ update_coords(fplog, step, ir, mdatoms, state, &f, fcd,
ekind, M, upd, etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
- /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */
+ /* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
/* are the small terms in the shake_vir here due
* to numerical errors, or are they important
* physically? I'm thinking they are just errors, but not completely sure.
* For now, will call without actually constraining, constr=NULL*/
- update_constraints(fplog, step, NULL, ir, mdatoms,
- state, fr->bMolPBC, graph, f,
+ update_constraints(fplog, step, nullptr, ir, mdatoms,
+ state, fr->bMolPBC, graph, &f,
&top->idef, tmp_vir,
- cr, nrnb, wcycle, upd, NULL,
+ cr, nrnb, wcycle, upd, nullptr,
FALSE, bCalcVir);
}
if (EI_VV(ir->eI))
else if (graph)
{
/* Need to unshift here */
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
- if (vsite != NULL)
+ if (vsite != nullptr)
{
wallcycle_start(wcycle, ewcVSITECONSTR);
- if (graph != NULL)
+ if (graph != nullptr)
{
- shift_self(graph, state->box, state->x);
+ shift_self(graph, state->box, as_rvec_array(state->x.data()));
}
- construct_vsites(vsite, state->x, ir->delta_t, state->v,
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
top->idef.iparams, top->idef.il,
fr->ePBC, fr->bMolPBC, cr, state->box);
- if (graph != NULL)
+ if (graph != nullptr)
{
- unshift_self(graph, state->box, state->x);
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
}
wallcycle_stop(wcycle, ewcVSITECONSTR);
}
Currently done every step so that dhdl is correct in the .edr */
sum_dhdl(enerd, state->lambda, ir->fepvals);
}
- update_box(fplog, step, ir, mdatoms, state, f,
- pcoupl_mu, nrnb, upd);
+
+ update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
+ pres, force_vir, shake_vir,
+ parrinellorahmanMu,
+ state, nrnb, upd);
/* ################# END UPDATE STEP 2 ################# */
/* #### We now have r(t+dt) and v(t+dt/2) ############# */
bSumEkinhOld = TRUE;
}
- /* ######### BEGIN PREPARING EDR OUTPUT ########### */
-
- /* use the directly determined last velocity, not actually the averaged half steps */
- if (bTrotter && ir->eI == eiVV)
+ if (bCalcEner)
{
- enerd->term[F_EKIN] = last_ekin;
- }
- enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
+ /* ######### BEGIN PREPARING EDR OUTPUT ########### */
- if (EI_VV(ir->eI))
- {
- enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
- }
- else
- {
- enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir, state, &MassQ);
+ /* use the directly determined last velocity, not actually the averaged half steps */
+ if (bTrotter && ir->eI == eiVV)
+ {
+ enerd->term[F_EKIN] = last_ekin;
+ }
+ enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN];
+
+ if (integratorHasConservedEnergyQuantity(ir))
+ {
+ if (EI_VV(ir->eI))
+ {
+ enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity;
+ }
+ else
+ {
+ enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + NPT_energy(ir, state, &MassQ);
+ }
+ }
+ /* ######### END PREPARING EDR OUTPUT ########### */
}
- /* ######### END PREPARING EDR OUTPUT ########### */
/* Output stuff */
if (MASTER(cr))
if (fplog && do_log && bDoExpanded)
{
/* only needed if doing expanded ensemble */
- PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL,
- &state_global->dfhist, state->fep_state, ir->nstlog, step);
+ PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
+ state_global->dfhist, state->fep_state, ir->nstlog, step);
}
if (bCalcEner)
{
gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
gmx_bool do_or = do_per_step(step, ir->nstorireout);
- print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : NULL,
+ print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
step, t,
eprNORMAL, mdebin, fcd, groups, &(ir->opts));
do_per_step(step, ir->swap->nstswap))
{
bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle,
- bRerunMD ? rerun_fr.x : state->x,
+ bRerunMD ? rerun_fr.x : as_rvec_array(state->x.data()),
bRerunMD ? rerun_fr.box : state->box,
- top_global, MASTER(cr) && bVerbose, bRerunMD);
+ MASTER(cr) && bVerbose, bRerunMD);
if (bNeedRepartition && DOMAINDECOMP(cr))
{
vsite, constr,
nrnb, wcycle, FALSE);
shouldCheckNumberOfBondedInteractions = true;
- update_realloc(upd, state->nalloc);
+ update_realloc(upd, state->natoms);
}
bFirstStep = FALSE;
/* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
- if ( (membed != NULL) && (!bLastStep) )
+ if ( (membed != nullptr) && (!bLastStep) )
{
- rescale_membed(step_rel, membed, state_global->x);
+ rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
}
if (bRerunMD)
"resetting counters later in the run, e.g. with gmx "
"mdrun -resetstep.", step);
}
- reset_all_counters(fplog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
- use_GPU(fr->nbv) ? fr->nbv : NULL);
+ reset_all_counters(fplog, mdlog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting,
+ use_GPU(fr->nbv) ? fr->nbv : nullptr);
wcycle_set_reset_counters(wcycle, -1);
if (!(cr->duty & DUTY_PME))
{
if (bRerunMD && MASTER(cr))
{
- close_trj(status);
+ close_trx(status);
}
if (!(cr->duty & DUTY_PME))
if (bPMETune)
{
- pme_loadbal_done(pme_loadbal, cr, fplog, use_GPU(fr->nbv));
+ pme_loadbal_done(pme_loadbal, fplog, mdlog, use_GPU(fr->nbv));
}
done_shellfc(fplog, shellfc, step_rel);
- if (repl_ex_nst > 0 && MASTER(cr))
+ if (useReplicaExchange && MASTER(cr))
{
print_replica_exchange_statistics(fplog, repl_ex);
}
finish_swapcoords(ir->swap);
}
+ /* Do essential dynamics cleanup if needed. Close .edo file */
+ done_ed(&ed);
+
/* IMD cleanup, if bIMD is TRUE. */
IMD_finalize(ir->bIMD, ir->imd);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff