*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2011,2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2011-2019, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "shellfc.h"
#if GMX_FAHCORE
-#include "corewrap.h"
+# include "corewrap.h"
#endif
using gmx::SimulationSignaller;
// alias to avoid a large ripple of nearly useless changes.
// t_inputrec is being replaced by IMdpOptionsProvider, so this
// will go away eventually.
- t_inputrec *ir = inputrec;
- int64_t step, step_rel;
- double t, t0 = ir->init_t, lam0[efptNR];
- gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
- gmx_bool bNS, bNStList, bStopCM,
- bFirstStep, bInitStep, bLastStep = FALSE;
- gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
- gmx_bool do_ene, do_log, do_verbose;
- gmx_bool bMasterState;
- unsigned int force_flags;
- tensor force_vir = {{0}}, shake_vir = {{0}}, total_vir = {{0}},
- tmp_vir = {{0}}, pres = {{0}};
+ t_inputrec* ir = inputrec;
+ int64_t step, step_rel;
+ double t, t0 = ir->init_t, lam0[efptNR];
+ gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
+ gmx_bool bNS, bNStList, bStopCM, bFirstStep, bInitStep, bLastStep = FALSE;
+ gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
+ gmx_bool do_ene, do_log, do_verbose;
+ gmx_bool bMasterState;
+ unsigned int force_flags;
+ tensor force_vir = { { 0 } }, shake_vir = { { 0 } }, total_vir = { { 0 } }, tmp_vir = { { 0 } },
+ pres = { { 0 } };
int i, m;
rvec mu_tot;
matrix parrinellorahmanMu, M;
gmx_repl_ex_t repl_ex = nullptr;
gmx_localtop_t top;
- PaddedHostVector<gmx::RVec> f {};
+ PaddedHostVector<gmx::RVec> f{};
gmx_global_stat_t gstat;
- t_graph *graph = nullptr;
- gmx_shellfc_t *shellfc;
+ t_graph* graph = nullptr;
+ gmx_shellfc_t* shellfc;
gmx_bool bSumEkinhOld, bDoReplEx, bExchanged, bNeedRepartition;
gmx_bool bTemp, bPres, bTrotter;
real dvdl_constr;
std::vector<RVec> cbuf;
matrix lastbox;
- int lamnew = 0;
+ int lamnew = 0;
/* for FEP */
- int nstfep = 0;
- double cycles;
- real saved_conserved_quantity = 0;
- real last_ekin = 0;
- t_extmass MassQ;
- char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
+ int nstfep = 0;
+ double cycles;
+ real saved_conserved_quantity = 0;
+ real last_ekin = 0;
+ t_extmass MassQ;
+ char sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
/* PME load balancing data for GPU kernels */
- gmx_bool bPMETune = FALSE;
- gmx_bool bPMETunePrinting = FALSE;
+ gmx_bool bPMETune = FALSE;
+ gmx_bool bPMETunePrinting = FALSE;
- bool bInteractiveMDstep = false;
+ bool bInteractiveMDstep = false;
/* Domain decomposition could incorrectly miss a bonded
interaction, but checking for that requires a global
code. So we do that alongside the first global energy reduction
after a new DD is made. These variables handle whether the
check happens, and the result it returns. */
- bool shouldCheckNumberOfBondedInteractions = false;
- int totalNumberOfBondedInteractions = -1;
+ bool shouldCheckNumberOfBondedInteractions = false;
+ int totalNumberOfBondedInteractions = -1;
SimulationSignals signals;
// Most global communnication stages don't propagate mdrun
// turning it off is for convenience in benchmarking, which is
// something that should not show up in the general user
// interface.
- GMX_LOG(mdlog.info).asParagraph().
- appendText("The -noconfout functionality is deprecated, and may be removed in a future version.");
+ GMX_LOG(mdlog.info)
+ .asParagraph()
+ .appendText(
+ "The -noconfout functionality is deprecated, and may be removed in a "
+ "future version.");
}
/* md-vv uses averaged full step velocities for T-control
md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control)
md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */
- bTrotter = (EI_VV(ir->eI) && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
+ bTrotter = (EI_VV(ir->eI)
+ && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir) || inputrecNvtTrotter(ir)));
- const bool bRerunMD = false;
+ const bool bRerunMD = false;
- int nstglobalcomm = computeGlobalCommunicationPeriod(mdlog, ir, cr);
- bGStatEveryStep = (nstglobalcomm == 1);
+ int nstglobalcomm = computeGlobalCommunicationPeriod(mdlog, ir, cr);
+ bGStatEveryStep = (nstglobalcomm == 1);
- SimulationGroups *groups = &top_global->groups;
+ SimulationGroups* groups = &top_global->groups;
std::unique_ptr<EssentialDynamics> ed = nullptr;
if (opt2bSet("-ei", nfile, fnm))
{
/* Initialize essential dynamics sampling */
- ed = init_edsam(mdlog,
- opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm),
- top_global,
- ir, cr, constr,
- state_global, observablesHistory,
- oenv,
- startingBehavior);
+ ed = init_edsam(mdlog, opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm), top_global,
+ ir, cr, constr, state_global, observablesHistory, oenv, startingBehavior);
}
else if (observablesHistory->edsamHistory)
{
{
pleaseCiteCouplingAlgorithms(fplog, *ir);
}
- gmx_mdoutf *outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider, mdModulesNotifier, ir, top_global, oenv, wcycle,
- startingBehavior);
- gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work, mdoutf_get_fp_dhdl(outf), false, mdModulesNotifier);
+ gmx_mdoutf* outf = init_mdoutf(fplog, nfile, fnm, mdrunOptions, cr, outputProvider,
+ mdModulesNotifier, ir, top_global, oenv, wcycle, startingBehavior);
+ gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf), top_global, ir, pull_work,
+ mdoutf_get_fp_dhdl(outf), false, mdModulesNotifier);
gstat = global_stat_init(ir);
/* Check for polarizable models and flexible constraints */
- shellfc = init_shell_flexcon(fplog,
- top_global, constr ? constr->numFlexibleConstraints() : 0,
+ shellfc = init_shell_flexcon(fplog, top_global, constr ? constr->numFlexibleConstraints() : 0,
ir->nstcalcenergy, DOMAINDECOMP(cr));
{
double io = compute_io(ir, top_global->natoms, *groups, energyOutput.numEnergyTerms(), 1);
if ((io > 2000) && MASTER(cr))
{
- fprintf(stderr,
- "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
- io);
+ fprintf(stderr, "\nWARNING: This run will generate roughly %.0f Mb of data\n\n", io);
}
}
// Local state only becomes valid now.
std::unique_ptr<t_state> stateInstance;
- t_state * state;
+ t_state* state;
auto mdatoms = mdAtoms->mdatoms();
dd_init_local_state(cr->dd, state_global, state);
/* Distribute the charge groups over the nodes from the master node */
- dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
- state_global, *top_global, ir, imdSession,
- pull_work,
- state, &f, mdAtoms, &top, fr,
- vsite, constr,
+ dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1, state_global, *top_global, ir,
+ imdSession, pull_work, state, &f, mdAtoms, &top, fr, vsite, constr,
nrnb, nullptr, FALSE);
shouldCheckNumberOfBondedInteractions = true;
upd.setNumAtoms(state->natoms);
state = state_global;
/* Generate and initialize new topology */
- mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr,
- &graph, mdAtoms, constr, vsite, shellfc);
+ mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &graph, mdAtoms, constr, vsite, shellfc);
upd.setNumAtoms(state->natoms);
}
-/*****************************************************************************************/
-// TODO: The following block of code should be refactored, once:
-// 1. We have the useGpuForBufferOps variable set and available here and in do_force(...)
-// 2. The proper GPU syncronization is introduced, so that the H2D and D2H data copies can be performed in the separate
-// stream owned by the StatePropagatorDataGpu
- const auto &simulationWork = runScheduleWork->simulationWork;
+ /*****************************************************************************************/
+ // TODO: The following block of code should be refactored, once:
+ // 1. We have the useGpuForBufferOps variable set and available here and in do_force(...)
+ // 2. The proper GPU syncronization is introduced, so that the H2D and D2H data copies can be performed in the separate
+ // stream owned by the StatePropagatorDataGpu
+ const auto& simulationWork = runScheduleWork->simulationWork;
const bool useGpuForPme = simulationWork.useGpuPme;
const bool useGpuForNonbonded = simulationWork.useGpuNonbonded;
// Temporary solution to make sure that the buffer ops are offloaded when update is offloaded
- const bool useGpuForBufferOps = simulationWork.useGpuBufferOps;
- const bool useGpuForUpdate = simulationWork.useGpuUpdate;
+ const bool useGpuForBufferOps = simulationWork.useGpuBufferOps;
+ const bool useGpuForUpdate = simulationWork.useGpuUpdate;
- StatePropagatorDataGpu *stateGpu = fr->stateGpu;
+ StatePropagatorDataGpu* stateGpu = fr->stateGpu;
if (useGpuForUpdate)
{
- GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr), "Domain decomposition is not supported with the GPU update.\n");
+ GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr),
+ "Domain decomposition is not supported with the GPU update.\n");
GMX_RELEASE_ASSERT(useGpuForPme || (useGpuForNonbonded && simulationWork.useGpuBufferOps),
- "Either PME or short-ranged non-bonded interaction tasks must run on the GPU to use GPU update.\n");
- GMX_RELEASE_ASSERT(ir->eI == eiMD, "Only the md integrator is supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(ir->etc != etcNOSEHOOVER, "Nose-Hoover temperature coupling is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(ir->epc == epcNO, "Pressure coupling is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(!mdatoms->haveVsites, "Virtual sites are not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(ed == nullptr, "Essential dynamics is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(!ir->bPull && !ir->pull, "Pulling is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(fcd->orires.nr == 0, "Orientation restraints are not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(ir->efep == efepNO, "Free energy perturbations are not supported with the GPU update.");
+ "Either PME or short-ranged non-bonded interaction tasks must run on "
+ "the GPU to use GPU update.\n");
+ GMX_RELEASE_ASSERT(ir->eI == eiMD,
+ "Only the md integrator is supported with the GPU update.\n");
+ GMX_RELEASE_ASSERT(
+ ir->etc != etcNOSEHOOVER,
+ "Nose-Hoover temperature coupling is not supported with the GPU update.\n");
+ GMX_RELEASE_ASSERT(ir->epc == epcNO,
+ "Pressure coupling is not supported with the GPU update.\n");
+ GMX_RELEASE_ASSERT(!mdatoms->haveVsites,
+ "Virtual sites are not supported with the GPU update.\n");
+ GMX_RELEASE_ASSERT(ed == nullptr,
+ "Essential dynamics is not supported with the GPU update.\n");
+ GMX_RELEASE_ASSERT(!ir->bPull && !ir->pull,
+ "Pulling is not supported with the GPU update.\n");
+ GMX_RELEASE_ASSERT(fcd->orires.nr == 0,
+ "Orientation restraints are not supported with the GPU update.\n");
+ GMX_RELEASE_ASSERT(ir->efep == efepNO,
+ "Free energy perturbations are not supported with the GPU update.");
GMX_RELEASE_ASSERT(graph == nullptr, "The graph is not supported with GPU update.");
if (constr != nullptr && constr->numConstraintsTotal() > 0)
{
- GMX_LOG(mdlog.info).asParagraph().
- appendText("Updating coordinates and applying constraints on the GPU.");
+ GMX_LOG(mdlog.info)
+ .asParagraph()
+ .appendText("Updating coordinates and applying constraints on the GPU.");
}
else
{
- GMX_LOG(mdlog.info).asParagraph().
- appendText("Updating coordinates on the GPU.");
+ GMX_LOG(mdlog.info).asParagraph().appendText("Updating coordinates on the GPU.");
}
- integrator = std::make_unique<UpdateConstrainCuda>(*ir, *top_global, stateGpu->getUpdateStream(), stateGpu->xUpdatedOnDevice());
+ integrator = std::make_unique<UpdateConstrainCuda>(
+ *ir, *top_global, stateGpu->getUpdateStream(), stateGpu->xUpdatedOnDevice());
}
if (useGpuForPme || (useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
{
changePinningPolicy(&state->v, PinningPolicy::PinnedIfSupported);
}
-/*****************************************************************************************/
+ /*****************************************************************************************/
// NOTE: The global state is no longer used at this point.
// But state_global is still used as temporary storage space for writing
/* Check nstexpanded here, because the grompp check was broken */
if (ir->expandedvals->nstexpanded % ir->nstcalcenergy != 0)
{
- gmx_fatal(FARGS, "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
+ gmx_fatal(FARGS,
+ "With expanded ensemble, nstexpanded should be a multiple of nstcalcenergy");
}
- init_expanded_ensemble(startingBehavior != StartingBehavior::NewSimulation,
- ir, state->dfhist);
+ init_expanded_ensemble(startingBehavior != StartingBehavior::NewSimulation, ir, state->dfhist);
}
if (MASTER(cr))
{
- EnergyElement::initializeEnergyHistory(
- startingBehavior, observablesHistory, &energyOutput);
+ EnergyElement::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput);
}
preparePrevStepPullCom(ir, pull_work, mdatoms, state, state_global, cr,
startingBehavior != StartingBehavior::NewSimulation);
// TODO: Remove this by converting AWH into a ForceProvider
- auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms, startingBehavior != StartingBehavior::NewSimulation,
- shellfc != nullptr,
- opt2fn("-awh", nfile, fnm), pull_work);
+ auto awh = prepareAwhModule(fplog, *ir, state_global, cr, ms,
+ startingBehavior != StartingBehavior::NewSimulation,
+ shellfc != nullptr, opt2fn("-awh", nfile, fnm), pull_work);
const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
if (useReplicaExchange && MASTER(cr))
{
- repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
- replExParams);
+ repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir, replExParams);
}
/* PME tuning is only supported in the Verlet scheme, with PME for
* Coulomb. It is not supported with only LJ PME. */
- bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) &&
- !mdrunOptions.reproducible && ir->cutoff_scheme != ecutsGROUP);
+ bPMETune = (mdrunOptions.tunePme && EEL_PME(fr->ic->eeltype) && !mdrunOptions.reproducible
+ && ir->cutoff_scheme != ecutsGROUP);
- pme_load_balancing_t *pme_loadbal = nullptr;
+ pme_load_balancing_t* pme_loadbal = nullptr;
if (bPMETune)
{
- pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box,
- *fr->ic, *fr->nbv, fr->pmedata, fr->nbv->useGpu(),
- &bPMETunePrinting);
+ pme_loadbal_init(&pme_loadbal, cr, mdlog, *ir, state->box, *fr->ic, *fr->nbv, fr->pmedata,
+ fr->nbv->useGpu(), &bPMETunePrinting);
}
if (!ir->bContinuation)
/* Set the velocities of vsites, shells and frozen atoms to zero */
for (i = 0; i < mdatoms->homenr; i++)
{
- if (mdatoms->ptype[i] == eptVSite ||
- mdatoms->ptype[i] == eptShell)
+ if (mdatoms->ptype[i] == eptVSite || mdatoms->ptype[i] == eptShell)
{
clear_rvec(v[i]);
}
if (constr)
{
/* Constrain the initial coordinates and velocities */
- do_constrain_first(fplog, constr, ir, mdatoms,
- state->natoms,
- state->x.arrayRefWithPadding(),
- state->v.arrayRefWithPadding(),
- state->box, state->lambda[efptBONDED]);
+ do_constrain_first(fplog, constr, ir, mdatoms, state->natoms, state->x.arrayRefWithPadding(),
+ state->v.arrayRefWithPadding(), state->box, state->lambda[efptBONDED]);
}
if (vsite)
{
/* Construct the virtual sites for the initial configuration */
- construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr,
- top.idef.iparams, top.idef.il,
- fr->ePBC, fr->bMolPBC, cr, state->box);
+ construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, nullptr, top.idef.iparams,
+ top.idef.il, fr->ePBC, fr->bMolPBC, cr, state->box);
}
}
restore_ekinstate_from_state(cr, ekind, &state_global->ekinstate);
}
- unsigned int cglo_flags = (CGLO_TEMPERATURE | CGLO_GSTAT
- | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
- | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0)
- | (hasReadEkinState ? CGLO_READEKIN : 0));
+ unsigned int cglo_flags =
+ (CGLO_TEMPERATURE | CGLO_GSTAT | (EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
+ | (EI_VV(ir->eI) ? CGLO_CONSTRAINT : 0) | (hasReadEkinState ? CGLO_READEKIN : 0));
bSumEkinhOld = FALSE;
cglo_flags_iteration |= CGLO_STOPCM;
cglo_flags_iteration &= ~CGLO_TEMPERATURE;
}
- compute_globals(gstat, cr, ir, fr, ekind,
- state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
- mdatoms, nrnb, &vcm,
- nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
- constr, &nullSignaller, state->box,
- &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags_iteration
- | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
+ compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
+ state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, nullptr, enerd,
+ force_vir, shake_vir, total_vir, pres, mu_tot, constr, &nullSignaller,
+ state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ cglo_flags_iteration
+ | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
+ : 0));
if (cglo_flags_iteration & CGLO_STOPCM)
{
/* At initialization, do not pass x with acceleration-correction mode
* to avoid (incorrect) correction of the initial coordinates.
*/
- rvec *xPtr = nullptr;
+ rvec* xPtr = nullptr;
if (vcm.mode != ecmLINEAR_ACCELERATION_CORRECTION)
{
xPtr = state->x.rvec_array();
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
}
}
- checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
- top_global, &top, state->x.rvec_array(), state->box,
+ checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global, &top,
+ state->x.rvec_array(), state->box,
&shouldCheckNumberOfBondedInteractions);
if (ir->eI == eiVVAK)
{
kinetic energy calculation. This minimized excess variables, but
perhaps loses some logic?*/
- compute_globals(gstat, cr, ir, fr, ekind,
- state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
- mdatoms, nrnb, &vcm,
- nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
- constr, &nullSignaller, state->box,
- nullptr, &bSumEkinhOld,
- cglo_flags & ~CGLO_PRESSURE);
+ compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
+ state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, nullptr, enerd,
+ force_vir, shake_vir, total_vir, pres, mu_tot, constr, &nullSignaller,
+ state->box, nullptr, &bSumEkinhOld, cglo_flags & ~CGLO_PRESSURE);
}
/* Calculate the initial half step temperature, and save the ekinh_old */
}
}
- /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter
- temperature control */
+ /* need to make an initiation call to get the Trotter variables set, as well as other constants
+ for non-trotter temperature control */
auto trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter);
if (MASTER(cr))
{
if (constr && ir->eConstrAlg == econtLINCS)
{
- fprintf(fplog,
- "RMS relative constraint deviation after constraining: %.2e\n",
+ fprintf(fplog, "RMS relative constraint deviation after constraining: %.2e\n",
constr->rmsd());
}
if (EI_STATE_VELOCITY(ir->eI))
}
char tbuf[20];
- fprintf(stderr, "starting mdrun '%s'\n",
- *(top_global->name));
+ fprintf(stderr, "starting mdrun '%s'\n", *(top_global->name));
if (ir->nsteps >= 0)
{
- sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t);
+ sprintf(tbuf, "%8.1f", (ir->init_step + ir->nsteps) * ir->delta_t);
}
else
{
if (ir->init_step > 0)
{
fprintf(stderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n",
- gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf,
- gmx_step_str(ir->init_step, sbuf2),
- ir->init_step*ir->delta_t);
+ gmx_step_str(ir->init_step + ir->nsteps, sbuf), tbuf,
+ gmx_step_str(ir->init_step, sbuf2), ir->init_step * ir->delta_t);
}
else
{
- fprintf(stderr, "%s steps, %s ps.\n",
- gmx_step_str(ir->nsteps, sbuf), tbuf);
+ fprintf(stderr, "%s steps, %s ps.\n", gmx_step_str(ir->nsteps, sbuf), tbuf);
}
fprintf(fplog, "\n");
}
#if GMX_FAHCORE
/* safest point to do file checkpointing is here. More general point would be immediately before integrator call */
- int chkpt_ret = fcCheckPointParallel( cr->nodeid,
- NULL, 0);
+ int chkpt_ret = fcCheckPointParallel(cr->nodeid, NULL, 0);
if (chkpt_ret == 0)
{
- gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0 );
+ gmx_fatal(3, __FILE__, __LINE__, "Checkpoint error on step %d\n", 0);
}
#endif
*
************************************************************/
- bFirstStep = TRUE;
+ bFirstStep = TRUE;
/* Skip the first Nose-Hoover integration when we get the state from tpx */
bInitStep = startingBehavior == StartingBehavior::NewSimulation || EI_VV(ir->eI);
bSumEkinhOld = FALSE;
{
// TODO This implementation of ensemble orientation restraints is nasty because
// a user can't just do multi-sim with single-sim orientation restraints.
- bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
- bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
+ bool usingEnsembleRestraints =
+ (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
+ bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
// Replica exchange, ensemble restraints and AWH need all
// simulations to remain synchronized, so they need
// Inter-simulation signal communication does not need to happen
// often, so we use a minimum of 200 steps to reduce overhead.
const int c_minimumInterSimulationSignallingInterval = 200;
- nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
+ nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1) / nstglobalcomm)
+ * nstglobalcomm;
}
}
auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
- compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
- MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
- mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
+ compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
+ MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
+ mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
auto checkpointHandler = std::make_unique<CheckpointHandler>(
- compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]),
- simulationsShareState, ir->nstlist == 0, MASTER(cr),
- mdrunOptions.writeConfout, mdrunOptions.checkpointOptions.period);
+ compat::make_not_null<SimulationSignal*>(&signals[eglsCHKPT]), simulationsShareState,
+ ir->nstlist == 0, MASTER(cr), mdrunOptions.writeConfout,
+ mdrunOptions.checkpointOptions.period);
const bool resetCountersIsLocal = true;
auto resetHandler = std::make_unique<ResetHandler>(
- compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]), !resetCountersIsLocal,
- ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
- mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
+ compat::make_not_null<SimulationSignal*>(&signals[eglsRESETCOUNTERS]),
+ !resetCountersIsLocal, ir->nsteps, MASTER(cr), mdrunOptions.timingOptions.resetHalfway,
+ mdrunOptions.maximumHoursToRun, mdlog, wcycle, walltime_accounting);
const DDBalanceRegionHandler ddBalanceRegionHandler(cr);
{
if (!multisim_int_all_are_equal(ms, ir->nsteps))
{
- GMX_LOG(mdlog.warning).appendText(
- "Note: The number of steps is not consistent across multi simulations,\n"
- "but we are proceeding anyway!");
+ 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(ms, ir->init_step))
{
- GMX_LOG(mdlog.warning).appendText(
- "Note: The initial step is not consistent across multi simulations,\n"
- "but we are proceeding anyway!");
+ GMX_LOG(mdlog.warning)
+ .appendText(
+ "Note: The initial step is not consistent across multi simulations,\n"
+ "but we are proceeding anyway!");
}
}
{
/* Determine if this is a neighbor search step */
- bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
+ bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0);
if (bPMETune && bNStList)
{
stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
}
/* PME grid + cut-off optimization with GPUs or PME nodes */
- pme_loadbal_do(pme_loadbal, cr,
- (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
- fplog, mdlog,
- *ir, fr, state->box, state->x,
- wcycle,
- step, step_rel,
+ pme_loadbal_do(pme_loadbal, cr, (mdrunOptions.verbose && MASTER(cr)) ? stderr : nullptr,
+ fplog, mdlog, *ir, fr, state->box, state->x, wcycle, step, step_rel,
&bPMETunePrinting);
}
wallcycle_start(wcycle, ewcSTEP);
bLastStep = (step_rel == ir->nsteps);
- t = t0 + step*ir->delta_t;
+ t = t0 + step * ir->delta_t;
// TODO Refactor this, so that nstfep does not need a default value of zero
if (ir->efep != efepNO || ir->bSimTemp)
/* find and set the current lambdas */
setCurrentLambdasLocal(step, ir->fepvals, lam0, state->lambda, state->fep_state);
- bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
- bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
- bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded)
- && (ir->bExpanded) && (step > 0) &&
- (startingBehavior == StartingBehavior::NewSimulation));
+ bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl);
+ bDoFEP = ((ir->efep != efepNO) && do_per_step(step, nstfep));
+ bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded) && (ir->bExpanded)
+ && (step > 0) && (startingBehavior == StartingBehavior::NewSimulation));
}
- bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep &&
- do_per_step(step, replExParams.exchangeInterval));
+ bDoReplEx = (useReplicaExchange && (step > 0) && !bLastStep
+ && do_per_step(step, replExParams.exchangeInterval));
if (doSimulatedAnnealing)
{
* Note that the || bLastStep can result in non-exact continuation
* beyond the last step. But we don't consider that to be an issue.
*/
- do_log =
- (do_per_step(step, ir->nstlog) ||
- (bFirstStep && startingBehavior == StartingBehavior::NewSimulation) ||
- bLastStep);
- do_verbose = mdrunOptions.verbose &&
- (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
+ do_log = (do_per_step(step, ir->nstlog)
+ || (bFirstStep && startingBehavior == StartingBehavior::NewSimulation) || bLastStep);
+ do_verbose = mdrunOptions.verbose
+ && (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
if (useGpuForUpdate && !bFirstStep)
{
if (DOMAINDECOMP(cr))
{
/* Repartition the domain decomposition */
- dd_partition_system(fplog, mdlog, step, cr,
- bMasterState, nstglobalcomm,
- state_global, *top_global, ir, imdSession,
- pull_work,
- state, &f, mdAtoms, &top, fr,
- vsite, constr,
- nrnb, wcycle,
- do_verbose && !bPMETunePrinting);
+ dd_partition_system(fplog, mdlog, step, cr, bMasterState, nstglobalcomm, state_global,
+ *top_global, ir, imdSession, pull_work, state, &f, mdAtoms, &top,
+ fr, vsite, constr, nrnb, wcycle, do_verbose && !bPMETunePrinting);
shouldCheckNumberOfBondedInteractions = true;
upd.setNumAtoms(state->natoms);
}
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(gstat, cr, ir, fr, ekind,
- state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
- mdatoms, nrnb, &vcm,
- wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
- constr, &nullSignaller, state->box,
- &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
+ state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, wcycle, enerd,
+ nullptr, nullptr, nullptr, nullptr, mu_tot, constr, &nullSignaller,
+ state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
- checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
- top_global, &top, state->x.rvec_array(), state->box,
+ checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global,
+ &top, state->x.rvec_array(), state->box,
&shouldCheckNumberOfBondedInteractions);
}
clear_mat(force_vir);
if (EI_VV(ir->eI) && (!bInitStep))
{
bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
- bCalcVir = bCalcEnerStep ||
- (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple)));
+ bCalcVir = bCalcEnerStep
+ || (ir->epc != epcNO
+ && (do_per_step(step, ir->nstpcouple) || do_per_step(step - 1, ir->nstpcouple)));
}
else
{
bCalcEnerStep = do_per_step(step, ir->nstcalcenergy);
- bCalcVir = bCalcEnerStep ||
- (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
+ bCalcVir = bCalcEnerStep || (ir->epc != epcNO && do_per_step(step, ir->nstpcouple));
}
bCalcEner = bCalcEnerStep;
}
/* Do we need global communication ? */
- bGStat = (bCalcVir || bCalcEner || bStopCM ||
- do_per_step(step, nstglobalcomm) ||
- (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step-1, nstglobalcomm)));
-
- force_flags = (GMX_FORCE_STATECHANGED |
- ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0) |
- GMX_FORCE_ALLFORCES |
- (bCalcVir ? GMX_FORCE_VIRIAL : 0) |
- (bCalcEner ? GMX_FORCE_ENERGY : 0) |
- (bDoFEP ? GMX_FORCE_DHDL : 0)
- );
+ bGStat = (bCalcVir || bCalcEner || bStopCM || do_per_step(step, nstglobalcomm)
+ || (EI_VV(ir->eI) && inputrecNvtTrotter(ir) && do_per_step(step - 1, nstglobalcomm)));
+
+ force_flags = (GMX_FORCE_STATECHANGED | ((inputrecDynamicBox(ir)) ? GMX_FORCE_DYNAMICBOX : 0)
+ | GMX_FORCE_ALLFORCES | (bCalcVir ? GMX_FORCE_VIRIAL : 0)
+ | (bCalcEner ? GMX_FORCE_ENERGY : 0) | (bDoFEP ? GMX_FORCE_DHDL : 0));
if (shellfc)
{
/* Now is the time to relax the shells */
- relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
- enforcedRotation, step,
- ir, imdSession, pull_work, bNS, force_flags, &top,
- constr, enerd, fcd,
- state->natoms,
- state->x.arrayRefWithPadding(),
- state->v.arrayRefWithPadding(),
- state->box,
- state->lambda,
- &state->hist,
- f.arrayRefWithPadding(), force_vir, mdatoms,
- nrnb, wcycle, graph,
- shellfc, fr, runScheduleWork, t, mu_tot,
- vsite,
- ddBalanceRegionHandler);
+ relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir,
+ imdSession, pull_work, bNS, force_flags, &top, constr, enerd, fcd,
+ state->natoms, state->x.arrayRefWithPadding(),
+ state->v.arrayRefWithPadding(), state->box, state->lambda, &state->hist,
+ f.arrayRefWithPadding(), force_vir, mdatoms, nrnb, wcycle, graph,
+ shellfc, fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler);
}
else
{
- /* The AWH history need to be saved _before_ doing force calculations where the AWH bias is updated
- (or the AWH update will be performed twice for one step when continuing). It would be best to
- call this update function from do_md_trajectory_writing but that would occur after do_force.
- One would have to divide the update_awh function into one function applying the AWH force
- and one doing the AWH bias update. The update AWH bias function could then be called after
- do_md_trajectory_writing (then containing update_awh_history).
- The checkpointing will in the future probably moved to the start of the md loop which will
- rid of this issue. */
+ /* The AWH history need to be saved _before_ doing force calculations where the AWH bias
+ is updated (or the AWH update will be performed twice for one step when continuing).
+ It would be best to call this update function from do_md_trajectory_writing but that
+ would occur after do_force. One would have to divide the update_awh function into one
+ function applying the AWH force and one doing the AWH bias update. The update AWH
+ bias function could then be called after do_md_trajectory_writing (then containing
+ update_awh_history). The checkpointing will in the future probably moved to the start
+ of the md loop which will rid of this issue. */
if (awh && checkpointHandler->isCheckpointingStep() && MASTER(cr))
{
awh->updateHistory(state_global->awhHistory.get());
* This is parallellized as well, and does communication too.
* Check comments in sim_util.c
*/
- do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession,
- pull_work,
- step, nrnb, wcycle, &top,
- state->box, state->x.arrayRefWithPadding(), &state->hist,
- f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd,
- state->lambda, graph,
+ do_force(fplog, cr, ms, ir, awh.get(), enforcedRotation, imdSession, pull_work, step,
+ nrnb, wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist,
+ f.arrayRefWithPadding(), force_vir, mdatoms, enerd, fcd, state->lambda, graph,
fr, runScheduleWork, vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
- (bNS ? GMX_FORCE_NS : 0) | force_flags,
- ddBalanceRegionHandler);
+ (bNS ? GMX_FORCE_NS : 0) | force_flags, ddBalanceRegionHandler);
}
// VV integrators do not need the following velocity half step
if (EI_VV(ir->eI) && (!bFirstStep || startingBehavior == StartingBehavior::NewSimulation))
/* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */
{
- rvec *vbuf = nullptr;
+ 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.rvec_array(), vbuf, 0, state->natoms); /* should make this better for parallelizing? */
+ copy_rvecn(state->v.rvec_array(), vbuf, 0,
+ state->natoms); /* should make this better for parallelizing? */
}
else
{
/* this is for NHC in the Ekin(t+dt/2) version of vv */
- trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1);
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ,
+ trotter_seq, ettTSEQ1);
}
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
- ekind, M, &upd, etrtVELOCITY1,
- cr, constr);
+ update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
+ etrtVELOCITY1, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
- constrain_velocities(step, nullptr,
- state,
- shake_vir,
- constr,
- bCalcVir, do_log, do_ene);
+ constrain_velocities(step, nullptr, state, shake_vir, constr, bCalcVir, do_log, do_ene);
wallcycle_start(wcycle, ewcUPDATE);
/* if VV, compute the pressure and constraints */
/* For VV2, we strictly only need this if using pressure
}
/* for vv, the first half of the integration actually corresponds to the previous step.
So we need information from the last step in the first half of the integration */
- if (bGStat || do_per_step(step-1, nstglobalcomm))
+ if (bGStat || do_per_step(step - 1, nstglobalcomm))
{
wallcycle_stop(wcycle, ewcUPDATE);
- compute_globals(gstat, cr, ir, fr, ekind,
- state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
- mdatoms, nrnb, &vcm,
- wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
- constr, &nullSignaller, state->box,
- &totalNumberOfBondedInteractions, &bSumEkinhOld,
- (bGStat ? CGLO_GSTAT : 0)
- | (bCalcEner ? CGLO_ENERGY : 0)
- | (bTemp ? CGLO_TEMPERATURE : 0)
- | (bPres ? CGLO_PRESSURE : 0)
- | (bPres ? CGLO_CONSTRAINT : 0)
- | (bStopCM ? CGLO_STOPCM : 0)
- | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
- | CGLO_SCALEEKIN
- );
+ compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
+ state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, wcycle, enerd,
+ force_vir, shake_vir, total_vir, pres, mu_tot, constr, &nullSignaller,
+ state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ (bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
+ | (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
+ | (bPres ? CGLO_CONSTRAINT : 0) | (bStopCM ? CGLO_STOPCM : 0)
+ | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
+ : 0)
+ | CGLO_SCALEEKIN);
/* explanation of above:
a) We compute Ekin at the full time step
if 1) we are using the AveVel Ekin, and it's not the
&shouldCheckNumberOfBondedInteractions);
if (bStopCM)
{
- process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(), state->v.rvec_array());
+ process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(),
+ state->v.rvec_array());
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
}
wallcycle_start(wcycle, ewcUPDATE);
{
if (bTrotter)
{
- m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */
- trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2);
+ m_add(force_vir, shake_vir,
+ total_vir); /* we need the un-dispersion corrected total vir here */
+ trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ,
+ trotter_seq, ettTSEQ2);
/* TODO This is only needed when we're about to write
* a checkpoint, because we use it after the restart
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, nullptr, (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);
}
}
/* We need the kinetic energy at minus the half step for determining
* the full step kinetic energy and possibly for T-coupling.*/
/* This may not be quite working correctly yet . . . . */
- compute_globals(gstat, cr, ir, fr, ekind,
- state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
- mdatoms, nrnb, &vcm,
- wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
- constr, &nullSignaller, state->box,
- nullptr, &bSumEkinhOld,
- CGLO_GSTAT | CGLO_TEMPERATURE);
+ compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(),
+ state->v.rvec_array(), state->box, state->lambda[efptVDW],
+ mdatoms, nrnb, &vcm, wcycle, enerd, nullptr, nullptr, nullptr,
+ nullptr, mu_tot, constr, &nullSignaller, state->box, nullptr,
+ &bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
wallcycle_start(wcycle, ewcUPDATE);
}
}
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.rvec_array(), mdatoms);
+ lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state,
+ state->dfhist, step, state->v.rvec_array(), mdatoms);
/* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
if (MASTER(cr))
{
// Copy coordinate from the GPU for the output if the update is offloaded and
// coordinates have not already been copied for i) search or ii) CPU force tasks.
- if (useGpuForUpdate && !bNS && !runScheduleWork->domainWork.haveCpuLocalForceWork &&
- (do_per_step(step, ir->nstxout) || do_per_step(step, ir->nstxout_compressed)))
+ if (useGpuForUpdate && !bNS && !runScheduleWork->domainWork.haveCpuLocalForceWork
+ && (do_per_step(step, ir->nstxout) || do_per_step(step, ir->nstxout_compressed)))
{
stateGpu->copyCoordinatesFromGpu(state->x, AtomLocality::Local);
stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
// copy call in do_force(...).
// NOTE: The forces should not be copied here if the vsites are present, since they were modified
// on host after the D2H copy in do_force(...).
- if (runScheduleWork->stepWork.useGpuFBufferOps && (simulationWork.useGpuUpdate && !vsite) && do_per_step(step, ir->nstfout))
+ if (runScheduleWork->stepWork.useGpuFBufferOps && (simulationWork.useGpuUpdate && !vsite)
+ && do_per_step(step, ir->nstfout))
{
stateGpu->copyForcesFromGpu(ArrayRef<RVec>(f), AtomLocality::Local);
stateGpu->waitForcesReadyOnHost(AtomLocality::Local);
* coordinates at time t. We must output all of this before
* the update.
*/
- do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
- ir, state, state_global, observablesHistory,
- top_global, fr,
- outf, energyOutput, ekind, f,
- checkpointHandler->isCheckpointingStep(),
- bRerunMD, bLastStep,
- mdrunOptions.writeConfout,
- bSumEkinhOld);
+ do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, state_global,
+ observablesHistory, top_global, fr, outf, energyOutput, ekind, f,
+ checkpointHandler->isCheckpointingStep(), bRerunMD, bLastStep,
+ mdrunOptions.writeConfout, bSumEkinhOld);
/* Check if IMD step and do IMD communication, if bIMD is TRUE. */
bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x.rvec_array(), t);
/* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
- if (startingBehavior != StartingBehavior::NewSimulation &&
- (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
+ if (startingBehavior != StartingBehavior::NewSimulation
+ && (inputrecNptTrotter(ir) || inputrecNphTrotter(ir)))
{
copy_mat(state->svir_prev, shake_vir);
copy_mat(state->fvir_prev, force_vir);
/* ######### START SECOND UPDATE STEP ################# */
- /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled
- in preprocessing */
+ /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen
+ controlled in preprocessing */
if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
{
/* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
if (constr && bIfRandomize)
{
- constrain_velocities(step, nullptr,
- state,
- tmp_vir,
- constr,
- bCalcVir, do_log, do_ene);
+ constrain_velocities(step, nullptr, state, tmp_vir, constr, bCalcVir, do_log, do_ene);
}
}
/* Box is changed in update() when we do pressure coupling,
else
{
update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
- update_pcouple_before_coordinates(fplog, step, ir, state,
- parrinellorahmanMu, M,
- bInitStep);
+ update_pcouple_before_coordinates(fplog, step, ir, state, parrinellorahmanMu, M, bInitStep);
}
if (EI_VV(ir->eI))
{
/* velocity half-step update */
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
- ekind, M, &upd, etrtVELOCITY2,
- cr, constr);
+ update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
+ etrtVELOCITY2, cr, constr);
}
/* Above, initialize just copies ekinh into ekin,
* energies of two subsequent steps. Therefore we need to compute the
* half step kinetic energy also if we need energies at the next step.
*/
- const bool needHalfStepKineticEnergy = (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm));
+ const bool needHalfStepKineticEnergy = (!EI_VV(ir->eI) && do_per_step(step + 1, nstglobalcomm));
if (useGpuForUpdate)
{
if (bNS)
{
- integrator->set(stateGpu->getCoordinates(), stateGpu->getVelocities(), stateGpu->getForces(),
- top.idef, *mdatoms, ekind->ngtc);
+ integrator->set(stateGpu->getCoordinates(), stateGpu->getVelocities(),
+ stateGpu->getForces(), top.idef, *mdatoms, ekind->ngtc);
t_pbc pbc;
set_pbc(&pbc, epbcXYZ, state->box);
integrator->setPbc(&pbc);
stateGpu->copyForcesToGpu(ArrayRef<RVec>(f), AtomLocality::Local);
}
- bool doTempCouple = (ir->etc != etcNO && do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
- bool doParrinelloRahman = (ir->epc == epcPARRINELLORAHMAN && do_per_step(step + ir->nstpcouple - 1, ir->nstpcouple));
+ bool doTempCouple =
+ (ir->etc != etcNO && do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
+ bool doParrinelloRahman = (ir->epc == epcPARRINELLORAHMAN
+ && do_per_step(step + ir->nstpcouple - 1, ir->nstpcouple));
// This applies Leap-Frog, LINCS and SETTLE in succession
- integrator->integrate(stateGpu->getForcesReadyOnDeviceEvent(AtomLocality::Local, runScheduleWork->stepWork.useGpuFBufferOps),
- ir->delta_t, true, bCalcVir, shake_vir,
- doTempCouple, ekind->tcstat,
- doParrinelloRahman, ir->nstpcouple*ir->delta_t, M);
+ integrator->integrate(stateGpu->getForcesReadyOnDeviceEvent(
+ AtomLocality::Local, runScheduleWork->stepWork.useGpuFBufferOps),
+ ir->delta_t, true, bCalcVir, shake_vir, doTempCouple,
+ ekind->tcstat, doParrinelloRahman, ir->nstpcouple * ir->delta_t, M);
// Copy velocities D2H after update if:
// - Globals are computed this step (includes the energy output steps).
}
else
{
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
- ekind, M, &upd, etrtPOSITION, cr, constr);
+ update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
+ etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
- constrain_coordinates(step, &dvdl_constr, state,
- shake_vir,
- &upd, constr,
- bCalcVir, do_log, do_ene);
+ constrain_coordinates(step, &dvdl_constr, state, shake_vir, &upd, constr, bCalcVir,
+ do_log, do_ene);
- update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state,
- cr, nrnb, wcycle, &upd, constr, do_log, do_ene);
- finish_update(ir, mdatoms,
- state, graph,
- nrnb, wcycle, &upd, constr);
+ update_sd_second_half(step, &dvdl_constr, ir, mdatoms, state, cr, nrnb, wcycle, &upd,
+ constr, do_log, do_ene);
+ finish_update(ir, mdatoms, state, graph, nrnb, wcycle, &upd, constr);
}
if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
{
/* erase F_EKIN and F_TEMP here? */
/* just compute the kinetic energy at the half step to perform a trotter step */
- compute_globals(gstat, cr, ir, fr, ekind,
- state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
- mdatoms, nrnb, &vcm,
- wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
- constr, &nullSignaller, lastbox,
- nullptr, &bSumEkinhOld,
- (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
- );
+ compute_globals(gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
+ state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, wcycle, enerd,
+ force_vir, shake_vir, total_vir, pres, mu_tot, constr, &nullSignaller, lastbox,
+ 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 */
std::copy(cbuf.begin(), cbuf.end(), state->x.begin());
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd,
- ekind, M, &upd, etrtPOSITION, cr, constr);
+ update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
+ etrtPOSITION, cr, constr);
wallcycle_stop(wcycle, ewcUPDATE);
/* do we need an extra constraint here? just need to copy out of as_rvec_array(state->v.data()) to upd->xp? */
* 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*/
- finish_update(ir, mdatoms,
- state, graph,
- nrnb, wcycle, &upd, nullptr);
+ finish_update(ir, mdatoms, state, graph, nrnb, wcycle, &upd, nullptr);
}
if (EI_VV(ir->eI))
{
this current solution is much better than
having it completely wrong.
*/
- enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr;
+ enerd->term[F_DVDL_CONSTR] += 2 * dvdl_constr;
}
else
{
shift_self(graph, state->box, state->x.rvec_array());
}
construct_vsites(vsite, state->x.rvec_array(), ir->delta_t, state->v.rvec_array(),
- top.idef.iparams, top.idef.il,
- fr->ePBC, fr->bMolPBC, cr, state->box);
+ top.idef.iparams, top.idef.il, fr->ePBC, fr->bMolPBC, cr, state->box);
if (graph != nullptr)
{
// kinetic energy at the step before we need to write
// the full-step kinetic energy
const bool needEkinAtNextStep =
- (!EI_VV(ir->eI) && (do_per_step(step + 1, nstglobalcomm) ||
- step_rel + 1 == ir->nsteps));
+ (!EI_VV(ir->eI) && (do_per_step(step + 1, nstglobalcomm) || step_rel + 1 == ir->nsteps));
if (bGStat || needEkinAtNextStep || doInterSimSignal)
{
bool doIntraSimSignal = true;
SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
- compute_globals(gstat, cr, ir, fr, ekind,
- state->x.rvec_array(), state->v.rvec_array(), state->box, state->lambda[efptVDW],
- mdatoms, nrnb, &vcm,
- wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
- constr, &signaller,
- lastbox,
- &totalNumberOfBondedInteractions, &bSumEkinhOld,
- (bGStat ? CGLO_GSTAT : 0)
- | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
+ compute_globals(
+ gstat, cr, ir, fr, ekind, state->x.rvec_array(), state->v.rvec_array(),
+ state->box, state->lambda[efptVDW], mdatoms, nrnb, &vcm, wcycle, enerd,
+ force_vir, shake_vir, total_vir, pres, mu_tot, constr, &signaller, lastbox,
+ &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ (bGStat ? CGLO_GSTAT : 0) | (!EI_VV(ir->eI) && bCalcEner ? CGLO_ENERGY : 0)
| (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
| (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
- | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0)
- | CGLO_CONSTRAINT
- | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
- );
+ | (!EI_VV(ir->eI) ? CGLO_PRESSURE : 0) | CGLO_CONSTRAINT
+ | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
+ : 0));
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
top_global, &top, state->x.rvec_array(), state->box,
&shouldCheckNumberOfBondedInteractions);
if (!EI_VV(ir->eI) && bStopCM)
{
- process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(), state->v.rvec_array());
+ process_and_stopcm_grp(fplog, &vcm, *mdatoms, state->x.rvec_array(),
+ state->v.rvec_array());
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
// TODO: The special case of removing CM motion should be dealt more gracefully
sum_dhdl(enerd, state->lambda, *ir->fepvals);
}
- update_pcouple_after_coordinates(fplog, step, ir, mdatoms,
- pres, force_vir, shake_vir,
- parrinellorahmanMu,
- state, 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) ############# */
if (fplog && do_log && bDoExpanded)
{
/* only needed if doing expanded ensemble */
- PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : nullptr,
+ PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals,
+ ir->bSimTemp ? ir->simtempvals : nullptr,
state_global->dfhist, state->fep_state, ir->nstlog, step);
}
if (bCalcEner)
{
- energyOutput.addDataAtEnergyStep(bDoDHDL, bCalcEnerStep,
- t, mdatoms->tmass, enerd, state,
- ir->fepvals, ir->expandedvals, lastbox,
- shake_vir, force_vir, total_vir, pres,
- ekind, mu_tot, constr);
+ energyOutput.addDataAtEnergyStep(bDoDHDL, bCalcEnerStep, t, mdatoms->tmass, enerd, state,
+ ir->fepvals, ir->expandedvals, lastbox, shake_vir,
+ force_vir, total_vir, pres, ekind, mu_tot, constr);
}
else
{
energyOutput.recordNonEnergyStep();
}
- gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
- gmx_bool do_or = do_per_step(step, ir->nstorireout);
+ gmx_bool do_dr = do_per_step(step, ir->nstdisreout);
+ gmx_bool do_or = do_per_step(step, ir->nstorireout);
if (doSimulatedAnnealing)
{
if (do_log || do_ene || do_dr || do_or)
{
energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
- do_log ? fplog : nullptr,
- step, t,
- fcd, awh.get());
+ do_log ? fplog : nullptr, step, t, fcd, awh.get());
}
if (ir->bPull)
state->fep_state = lamnew;
}
/* Print the remaining wall clock time for the run */
- if (isMasterSimMasterRank(ms, MASTER(cr)) &&
- (do_verbose || gmx_got_usr_signal()) &&
- !bPMETunePrinting)
+ if (isMasterSimMasterRank(ms, MASTER(cr)) && (do_verbose || gmx_got_usr_signal()) && !bPMETunePrinting)
{
if (shellfc)
{
* Not done in last step since trajectory writing happens before this call
* in the MD loop and exchanges would be lost anyway. */
bNeedRepartition = FALSE;
- if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep &&
- do_per_step(step, ir->swap->nstswap))
+ if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep && do_per_step(step, ir->swap->nstswap))
{
- bNeedRepartition = do_swapcoords(cr, step, t, ir, swap, wcycle,
- as_rvec_array(state->x.data()),
- state->box,
- MASTER(cr) && mdrunOptions.verbose,
- bRerunMD);
+ bNeedRepartition =
+ do_swapcoords(cr, step, t, ir, swap, wcycle, as_rvec_array(state->x.data()),
+ state->box, MASTER(cr) && mdrunOptions.verbose, bRerunMD);
if (bNeedRepartition && DOMAINDECOMP(cr))
{
bExchanged = FALSE;
if (bDoReplEx)
{
- bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
- state_global, enerd,
- state, step, t);
+ bExchanged = replica_exchange(fplog, cr, ms, repl_ex, state_global, enerd, state, step, t);
}
- if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr) )
+ if ((bExchanged || bNeedRepartition) && DOMAINDECOMP(cr))
{
- dd_partition_system(fplog, mdlog, step, cr, TRUE, 1,
- state_global, *top_global, ir, imdSession,
- pull_work,
- state, &f, mdAtoms, &top, fr,
- vsite, constr,
+ dd_partition_system(fplog, mdlog, step, cr, TRUE, 1, state_global, *top_global, ir,
+ imdSession, pull_work, state, &f, mdAtoms, &top, fr, vsite, constr,
nrnb, wcycle, FALSE);
shouldCheckNumberOfBondedInteractions = true;
upd.setNumAtoms(state->natoms);
}
- bFirstStep = FALSE;
- bInitStep = FALSE;
+ bFirstStep = FALSE;
+ bInitStep = FALSE;
/* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */
/* With all integrators, except VV, we need to retain the pressure
* at the current step for coupling at the next step.
*/
- if ((state->flags & (1U<<estPRES_PREV)) &&
- (bGStatEveryStep ||
- (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
+ if ((state->flags & (1U << estPRES_PREV))
+ && (bGStatEveryStep || (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
{
/* Store the pressure in t_state for pressure coupling
* at the next MD step.
/* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
- if ( (membed != nullptr) && (!bLastStep) )
+ if ((membed != nullptr) && (!bLastStep))
{
rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
}
step++;
step_rel++;
- resetHandler->resetCounters(
- step, step_rel, mdlog, fplog, cr, fr->nbv.get(),
- nrnb, fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
+ resetHandler->resetCounters(step, step_rel, mdlog, fplog, cr, fr->nbv.get(), nrnb,
+ fr->pmedata, pme_loadbal, wcycle, walltime_accounting);
/* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
imdSession->updateEnergyRecordAndSendPositionsAndEnergies(bInteractiveMDstep, step, bCalcEner);
-
}
/* End of main MD loop */
}
global_stat_destroy(gstat);
-
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff