#include <numeric>
#include "gromacs/applied_forces/awh/awh.h"
+#include "gromacs/applied_forces/awh/read_params.h"
#include "gromacs/commandline/filenm.h"
#include "gromacs/domdec/collect.h"
#include "gromacs/domdec/dlbtiming.h"
#include "gromacs/domdec/domdec_network.h"
#include "gromacs/domdec/domdec_struct.h"
#include "gromacs/domdec/gpuhaloexchange.h"
+#include "gromacs/domdec/localtopologychecker.h"
#include "gromacs/domdec/mdsetup.h"
#include "gromacs/domdec/partition.h"
#include "gromacs/essentialdynamics/edsam.h"
#include "gromacs/mdlib/update_vv.h"
#include "gromacs/mdlib/vcm.h"
#include "gromacs/mdlib/vsite.h"
+#include "gromacs/mdrunutility/freeenergy.h"
#include "gromacs/mdrunutility/handlerestart.h"
#include "gromacs/mdrunutility/multisim.h"
#include "gromacs/mdrunutility/printtime.h"
#include "gromacs/mdtypes/mdrunoptions.h"
#include "gromacs/mdtypes/multipletimestepping.h"
#include "gromacs/mdtypes/observableshistory.h"
+#include "gromacs/mdtypes/observablesreducer.h"
#include "gromacs/mdtypes/pullhistory.h"
#include "gromacs/mdtypes/simulation_workload.h"
#include "gromacs/mdtypes/state.h"
// will go away eventually.
const t_inputrec* ir = inputrec;
- int64_t step, step_rel;
double t, t0 = ir->init_t;
gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
gmx_bool bNS = FALSE, bNStList, bStopCM, bFirstStep, bInitStep, bLastStep = FALSE;
- gmx_bool bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
+ gmx_bool bDoExpanded = FALSE;
gmx_bool do_ene, do_log, do_verbose;
gmx_bool bMasterState;
unsigned int force_flags;
matrix lastbox;
int lamnew = 0;
/* for FEP */
- int nstfep = 0;
double cycles;
real saved_conserved_quantity = 0;
real last_ekin = 0;
bool bInteractiveMDstep = false;
- /* Domain decomposition could incorrectly miss a bonded
- interaction, but checking for that requires a global
- communication stage, which does not otherwise happen in DD
- 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;
-
SimulationSignals signals;
// Most global communnication stages don't propagate mdrun
// signals, and will use this object to achieve that.
int nstglobalcomm = computeGlobalCommunicationPeriod(mdlog, ir, cr);
bGStatEveryStep = (nstglobalcomm == 1);
- const SimulationGroups* groups = &top_global->groups;
+ const SimulationGroups* groups = &top_global.groups;
std::unique_ptr<EssentialDynamics> ed = nullptr;
if (opt2bSet("-ei", nfile, fnm))
ed = init_edsam(mdlog,
opt2fn_null("-ei", nfile, fnm),
opt2fn("-eo", nfile, fnm),
- *top_global,
+ top_global,
*ir,
cr,
constr,
int* fep_state = MASTER(cr) ? &state_global->fep_state : nullptr;
gmx::ArrayRef<real> lambda = MASTER(cr) ? state_global->lambda : gmx::ArrayRef<real>();
- initialize_lambdas(fplog, *ir, MASTER(cr), fep_state, lambda);
+ initialize_lambdas(fplog,
+ ir->efep,
+ ir->bSimTemp,
+ *ir->fepvals,
+ ir->simtempvals->temperatures,
+ gmx::arrayRefFromArray(ir->opts.ref_t, ir->opts.ngtc),
+ MASTER(cr),
+ fep_state,
+ lambda);
Update upd(*ir, deform);
bool doSimulatedAnnealing = false;
{
}
const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
- const t_fcdata& fcdata = *fr->fcdata;
+ t_fcdata& fcdata = *fr->fcdata;
bool simulationsShareState = false;
int nstSignalComm = nstglobalcomm;
{
// 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 =
- (fcdata.disres->nsystems > 1) || ((ms != nullptr) && (fcdata.orires->nr != 0));
- bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
+ bool usingEnsembleRestraints = (fcdata.disres->nsystems > 1) || ((ms != nullptr) && fcdata.orires);
+ bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim() && (ms != nullptr));
// Replica exchange, ensemble restraints and AWH need all
// simulations to remain synchronized, so they need
mdrunOptions,
cr,
outputProvider,
- mdModulesNotifier,
+ mdModulesNotifiers,
ir,
top_global,
oenv,
ms);
gmx::EnergyOutput energyOutput(mdoutf_get_fp_ene(outf),
top_global,
- ir,
+ *ir,
pull_work,
mdoutf_get_fp_dhdl(outf),
false,
startingBehavior,
simulationsShareState,
- mdModulesNotifier);
+ mdModulesNotifiers);
gstat = global_stat_init(ir);
const auto& simulationWork = runScheduleWork->simulationWork;
const bool useGpuForPme = simulationWork.useGpuPme;
const bool useGpuForNonbonded = simulationWork.useGpuNonbonded;
- const bool useGpuForBufferOps = simulationWork.useGpuBufferOps;
const bool useGpuForUpdate = simulationWork.useGpuUpdate;
/* Check for polarizable models and flexible constraints */
top_global,
constr ? constr->numFlexibleConstraints() : 0,
ir->nstcalcenergy,
- DOMAINDECOMP(cr),
+ haveDDAtomOrdering(*cr),
useGpuForPme);
{
- double io = compute_io(ir, top_global->natoms, *groups, energyOutput.numEnergyTerms(), 1);
+ 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);
}
}
- // Local state only becomes valid now.
- std::unique_ptr<t_state> stateInstance;
- t_state* state;
-
- gmx_localtop_t top(top_global->ffparams);
-
- auto mdatoms = mdAtoms->mdatoms();
+ ObservablesReducer observablesReducer = observablesReducerBuilder->build();
- ForceBuffers f(fr->useMts,
- ((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
- ? PinningPolicy::PinnedIfSupported
- : PinningPolicy::CannotBePinned);
- if (DOMAINDECOMP(cr))
+ ForceBuffers f(simulationWork.useMts,
+ (simulationWork.useGpuFBufferOps || useGpuForUpdate) ? PinningPolicy::PinnedIfSupported
+ : PinningPolicy::CannotBePinned);
+ const t_mdatoms* md = mdAtoms->mdatoms();
+ if (haveDDAtomOrdering(*cr))
{
- stateInstance = std::make_unique<t_state>();
- state = stateInstance.get();
- dd_init_local_state(cr->dd, state_global, state);
+ // Local state only becomes valid now.
+ dd_init_local_state(*cr->dd, state_global, state);
/* Distribute the charge groups over the nodes from the master node */
dd_partition_system(fplog,
TRUE,
1,
state_global,
- *top_global,
- ir,
+ top_global,
+ *ir,
imdSession,
pull_work,
state,
&f,
mdAtoms,
- &top,
+ top,
fr,
vsite,
constr,
nrnb,
nullptr,
FALSE);
- shouldCheckNumberOfBondedInteractions = true;
- upd.setNumAtoms(state->natoms);
+ upd.updateAfterPartition(state->natoms,
+ md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
+ fr->longRangeNonbondeds->updateAfterPartition(*md);
}
else
{
state_change_natoms(state_global, state_global->natoms);
- /* Copy the pointer to the global state */
- state = state_global;
/* Generate and initialize new topology */
- mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &f, mdAtoms, constr, vsite, shellfc);
-
- upd.setNumAtoms(state->natoms);
+ mdAlgorithmsSetupAtomData(cr, *ir, top_global, top, fr, &f, mdAtoms, constr, vsite, shellfc);
+
+ upd.updateAfterPartition(state->natoms,
+ md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
+ fr->longRangeNonbondeds->updateAfterPartition(*md);
}
std::unique_ptr<UpdateConstrainGpu> integrator;
// TODO: the assertions below should be handled by UpdateConstraintsBuilder.
if (useGpuForUpdate)
{
- GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr) || ddUsesUpdateGroups(*cr->dd) || constr == nullptr
- || constr->numConstraintsTotal() == 0,
+ GMX_RELEASE_ASSERT(!haveDDAtomOrdering(*cr) || ddUsesUpdateGroups(*cr->dd)
+ || constr == nullptr || constr->numConstraintsTotal() == 0,
"Constraints in domain decomposition are only supported with update "
"groups if using GPU update.\n");
- GMX_RELEASE_ASSERT(ir->eConstrAlg != econtSHAKE || constr == nullptr
+ GMX_RELEASE_ASSERT(ir->eConstrAlg != ConstraintAlgorithm::Shake || constr == nullptr
|| constr->numConstraintsTotal() == 0,
"SHAKE is not supported with GPU update.");
- GMX_RELEASE_ASSERT(useGpuForPme || (useGpuForNonbonded && simulationWork.useGpuBufferOps),
+ GMX_RELEASE_ASSERT(useGpuForPme || (useGpuForNonbonded && simulationWork.useGpuXBufferOps),
"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,
+ GMX_RELEASE_ASSERT(ir->eI == IntegrationAlgorithm::MD,
"Only the md integrator is supported with the GPU update.\n");
GMX_RELEASE_ASSERT(
ir->etc != TemperatureCoupling::NoseHoover,
|| ir->epc == PressureCoupling::Berendsen || ir->epc == PressureCoupling::CRescale,
"Only Parrinello-Rahman, Berendsen, and C-rescale pressure coupling are supported "
"with the GPU update.\n");
- GMX_RELEASE_ASSERT(!mdatoms->haveVsites,
+ GMX_RELEASE_ASSERT(!md->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 || !pull_have_constraint(*ir->pull),
"Constraints pulling is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(fcdata.orires->nr == 0,
+ GMX_RELEASE_ASSERT(fcdata.orires == nullptr,
"Orientation restraints are not supported with the GPU update.\n");
GMX_RELEASE_ASSERT(
- ir->efep == efepNO
- || (!haveFepPerturbedMasses(*top_global) && !havePerturbedConstraints(*top_global)),
+ ir->efep == FreeEnergyPerturbationType::No
+ || (!haveFepPerturbedMasses(top_global) && !havePerturbedConstraints(top_global)),
"Free energy perturbation of masses and constraints are not supported with the GPU "
"update.");
"update-constraints.");
integrator = std::make_unique<UpdateConstrainGpu>(
*ir,
- *top_global,
+ top_global,
ekind->ngtc,
fr->deviceStreamManager->context(),
fr->deviceStreamManager->stream(gmx::DeviceStreamType::UpdateAndConstraints),
- stateGpu->xUpdatedOnDevice(),
wcycle);
+ stateGpu->setXUpdatedOnDeviceEvent(integrator->xUpdatedOnDeviceEvent());
+
integrator->setPbc(PbcType::Xyz, state->box);
}
- if (useGpuForPme || (useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
+ if (useGpuForPme || simulationWork.useGpuXBufferOps || useGpuForUpdate)
{
changePinningPolicy(&state->x, PinningPolicy::PinnedIfSupported);
}
// the global state to file and potentially for replica exchange.
// (Global topology should persist.)
- update_mdatoms(mdatoms, state->lambda[efptMASS]);
+ update_mdatoms(mdAtoms->mdatoms(), state->lambda[FreeEnergyPerturbationCouplingType::Mass]);
if (ir->bExpanded)
{
EnergyData::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput);
}
- preparePrevStepPullCom(
- ir, pull_work, mdatoms->massT, state, state_global, cr, startingBehavior != StartingBehavior::NewSimulation);
+ preparePrevStepPullCom(ir,
+ pull_work,
+ gmx::arrayRefFromArray(md->massT, md->nr),
+ state,
+ state_global,
+ cr,
+ startingBehavior != StartingBehavior::NewSimulation);
// TODO: Remove this by converting AWH into a ForceProvider
auto awh = prepareAwhModule(fplog,
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);
+ && ir->cutoff_scheme != CutoffScheme::Group);
pme_load_balancing_t* pme_loadbal = nullptr;
if (bPMETune)
if (!ir->bContinuation)
{
- if (state->flags & (1U << estV))
+ if (state->flags & enumValueToBitMask(StateEntry::V))
{
auto v = makeArrayRef(state->v);
/* Set the velocities of vsites, shells and frozen atoms to zero */
- for (i = 0; i < mdatoms->homenr; i++)
+ for (i = 0; i < md->homenr; i++)
{
- if (mdatoms->ptype[i] == eptShell)
+ if (md->ptype[i] == ParticleType::Shell)
{
clear_rvec(v[i]);
}
- else if (mdatoms->cFREEZE)
+ else if (md->cFREEZE)
{
for (m = 0; m < DIM; m++)
{
- if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m])
+ if (ir->opts.nFreeze[md->cFREEZE[i]][m])
{
v[i][m] = 0;
}
do_constrain_first(fplog,
constr,
ir,
- mdatoms->nr,
- mdatoms->homenr,
+ md->nr,
+ md->homenr,
state->x.arrayRefWithPadding(),
state->v.arrayRefWithPadding(),
state->box,
- state->lambda[efptBONDED]);
+ state->lambda[FreeEnergyPerturbationCouplingType::Bonded]);
}
}
- if (ir->efep != efepNO)
- {
- /* Set free energy calculation frequency as the greatest common
- * denominator of nstdhdl and repl_ex_nst. */
- nstfep = ir->fepvals->nstdhdl;
- if (ir->bExpanded)
- {
- nstfep = std::gcd(ir->expandedvals->nstexpanded, nstfep);
- }
- if (useReplicaExchange)
- {
- nstfep = std::gcd(replExParams.exchangeInterval, nstfep);
- }
- if (ir->bDoAwh)
- {
- nstfep = std::gcd(ir->awhParams->nstSampleCoord, nstfep);
- }
- }
+ const int nstfep = computeFepPeriod(*ir, replExParams);
/* Be REALLY careful about what flags you set here. You CANNOT assume
* this is the first step, since we might be restarting from a checkpoint,
* and in that case we should not do any modifications to the state.
*/
- bStopCM = (ir->comm_mode != ecmNO && !ir->bContinuation);
+ bStopCM = (ir->comm_mode != ComRemovalAlgorithm::No && !ir->bContinuation);
// When restarting from a checkpoint, it can be appropriate to
// initialize ekind from quantities in the checkpoint. Otherwise,
bSumEkinhOld = FALSE;
- t_vcm vcm(top_global->groups, *ir);
+ t_vcm vcm(top_global.groups, *ir);
reportComRemovalInfo(fplog, vcm);
+ int64_t step = ir->init_step;
+ int64_t step_rel = 0;
+
/* To minimize communication, compute_globals computes the COM velocity
* and the kinetic energy for the velocities without COM motion removed.
* Thus to get the kinetic energy without the COM contribution, we need
makeConstArrayRef(state->x),
makeConstArrayRef(state->v),
state->box,
- mdatoms,
+ md,
nrnb,
&vcm,
nullptr,
shake_vir,
total_vir,
pres,
- gmx::ArrayRef<real>{},
&nullSignaller,
state->box,
- &totalNumberOfBondedInteractions,
&bSumEkinhOld,
- cglo_flags_iteration
- | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
- : 0));
+ cglo_flags_iteration,
+ step,
+ &observablesReducer);
+ // Clean up after pre-step use of compute_globals()
+ observablesReducer.markAsReadyToReduce();
+
if (cglo_flags_iteration & CGLO_STOPCM)
{
/* At initialization, do not pass x with acceleration-correction mode
* to avoid (incorrect) correction of the initial coordinates.
*/
- auto x = (vcm.mode == ecmLINEAR_ACCELERATION_CORRECTION) ? ArrayRef<RVec>()
- : makeArrayRef(state->x);
- process_and_stopcm_grp(fplog, &vcm, *mdatoms, x, makeArrayRef(state->v));
- inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
+ auto x = (vcm.mode == ComRemovalAlgorithm::LinearAccelerationCorrection)
+ ? ArrayRef<RVec>()
+ : makeArrayRef(state->x);
+ process_and_stopcm_grp(fplog, &vcm, *md, x, makeArrayRef(state->v));
+ inc_nrnb(nrnb, eNR_STOPCM, md->homenr);
}
}
- checkNumberOfBondedInteractions(mdlog,
- cr,
- totalNumberOfBondedInteractions,
- top_global,
- &top,
- makeConstArrayRef(state->x),
- state->box,
- &shouldCheckNumberOfBondedInteractions);
- if (ir->eI == eiVVAK)
+ if (ir->eI == IntegrationAlgorithm::VVAK)
{
/* a second call to get the half step temperature initialized as well */
/* we do the same call as above, but turn the pressure off -- internally to
makeConstArrayRef(state->x),
makeConstArrayRef(state->v),
state->box,
- mdatoms,
+ md,
nrnb,
&vcm,
nullptr,
shake_vir,
total_vir,
pres,
- gmx::ArrayRef<real>{},
&nullSignaller,
state->box,
- nullptr,
&bSumEkinhOld,
- cglo_flags & ~CGLO_PRESSURE);
+ cglo_flags & ~CGLO_PRESSURE,
+ step,
+ &observablesReducer);
+ // Clean up after pre-step use of compute_globals()
+ observablesReducer.markAsReadyToReduce();
}
/* Calculate the initial half step temperature, and save the ekinh_old */
{
if (!ir->bContinuation)
{
- if (constr && ir->eConstrAlg == econtLINCS)
+ if (constr && ir->eConstrAlg == ConstraintAlgorithm::Lincs)
{
fprintf(fplog,
"RMS relative constraint deviation after constraining: %.2e\n",
if (EI_STATE_VELOCITY(ir->eI))
{
real temp = enerd->term[F_TEMP];
- if (ir->eI != eiVV)
+ if (ir->eI != IntegrationAlgorithm::VV)
{
/* Result of Ekin averaged over velocities of -half
* and +half step, while we only have -half step here.
}
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);
}
walltime_accounting_start_time(walltime_accounting);
- wallcycle_start(wcycle, ewcRUN);
+ wallcycle_start(wcycle, WallCycleCounter::Run);
print_start(fplog, cr, walltime_accounting, "mdrun");
/***********************************************************
bExchanged = FALSE;
bNeedRepartition = FALSE;
- step = ir->init_step;
- step_rel = 0;
-
auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]),
simulationsShareState,
simulationWork.useGpuPmePpCommunication);
}
- wallcycle_start(wcycle, ewcSTEP);
+ wallcycle_start(wcycle, WallCycleCounter::Step);
bLastStep = (step_rel == ir->nsteps);
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)
+ if (ir->efep != FreeEnergyPerturbationType::No || ir->bSimTemp)
{
/* find and set the current lambdas */
state->lambda = currentLambdas(step, *(ir->fepvals), 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)
&& (!bFirstStep));
}
}
/* Stop Center of Mass motion */
- bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm));
+ bStopCM = (ir->comm_mode != ComRemovalAlgorithm::No && do_per_step(step, ir->nstcomm));
/* Determine whether or not to do Neighbour Searching */
bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
do_verbose = mdrunOptions.verbose
&& (step % mdrunOptions.verboseStepPrintInterval == 0 || bFirstStep || bLastStep);
- if (useGpuForUpdate && !bFirstStep && bNS)
+ // On search steps, when doing the update on the GPU, copy
+ // the coordinates and velocities to the host unless they are
+ // already there (ie on the first step and after replica
+ // exchange).
+ if (useGpuForUpdate && bNS && !bFirstStep && !bExchanged)
{
- // Copy velocities from the GPU on search steps to keep a copy on host (device buffers are reinitialized).
stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
- stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
- // Copy coordinate from the GPU when needed at the search step.
- // NOTE: The cases when coordinates needed on CPU for force evaluation are handled in sim_utils.
- // NOTE: If the coordinates are to be written into output file they are also copied separately before the output.
stateGpu->copyCoordinatesFromGpu(state->x, AtomLocality::Local);
+ stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
}
if (vsite != nullptr)
{
// Virtual sites need to be updated before domain decomposition and forces are calculated
- wallcycle_start(wcycle, ewcVSITECONSTR);
+ wallcycle_start(wcycle, WallCycleCounter::VsiteConstr);
// md-vv calculates virtual velocities once it has full-step real velocities
vsite->construct(state->x,
state->v,
(!EI_VV(inputrec->eI) && needVirtualVelocitiesThisStep)
? VSiteOperation::PositionsAndVelocities
: VSiteOperation::Positions);
- wallcycle_stop(wcycle, ewcVSITECONSTR);
+ wallcycle_stop(wcycle, WallCycleCounter::VsiteConstr);
}
if (bNS && !(bFirstStep && ir->bContinuation))
if (correct_box(fplog, step, state->box))
{
bMasterState = TRUE;
- // If update is offloaded, it should be informed about the box size change
- if (useGpuForUpdate)
- {
- integrator->setPbc(PbcType::Xyz, state->box);
- }
}
}
- if (DOMAINDECOMP(cr) && bMasterState)
+ // If update is offloaded, and the box was changed either
+ // above or in a replica exchange on the previous step,
+ // the GPU Update object should be informed
+ if (useGpuForUpdate && (bMasterState || bExchanged))
+ {
+ integrator->setPbc(PbcType::Xyz, state->box);
+ }
+ if (haveDDAtomOrdering(*cr) && bMasterState)
{
dd_collect_state(cr->dd, state, state_global);
}
- if (DOMAINDECOMP(cr))
+ if (haveDDAtomOrdering(*cr))
{
/* Repartition the domain decomposition */
dd_partition_system(fplog,
bMasterState,
nstglobalcomm,
state_global,
- *top_global,
- ir,
+ top_global,
+ *ir,
imdSession,
pull_work,
state,
&f,
mdAtoms,
- &top,
+ top,
fr,
vsite,
constr,
nrnb,
wcycle,
do_verbose && !bPMETunePrinting);
- shouldCheckNumberOfBondedInteractions = true;
- upd.setNumAtoms(state->natoms);
+ upd.updateAfterPartition(state->natoms,
+ md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
+ fr->longRangeNonbondeds->updateAfterPartition(*md);
}
}
// Allocate or re-size GPU halo exchange object, if necessary
- if (bNS && havePPDomainDecomposition(cr) && simulationWork.useGpuHaloExchange)
+ if (bNS && simulationWork.havePpDomainDecomposition && simulationWork.useGpuHaloExchange)
{
GMX_RELEASE_ASSERT(fr->deviceStreamManager != nullptr,
"GPU device manager has to be initialized to use GPU "
fplog, step, t); /* can we improve the information printed here? */
}
- if (ir->efep != efepNO)
+ if (ir->efep != FreeEnergyPerturbationType::No)
{
- update_mdatoms(mdatoms, state->lambda[efptMASS]);
+ update_mdatoms(mdAtoms->mdatoms(), state->lambda[FreeEnergyPerturbationCouplingType::Mass]);
}
if (bExchanged)
/* 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 . . . . */
+ int cglo_flags = CGLO_GSTAT | CGLO_TEMPERATURE;
compute_globals(gstat,
cr,
ir,
makeConstArrayRef(state->x),
makeConstArrayRef(state->v),
state->box,
- mdatoms,
+ md,
nrnb,
&vcm,
wcycle,
nullptr,
nullptr,
nullptr,
- gmx::ArrayRef<real>{},
&nullSignaller,
state->box,
- &totalNumberOfBondedInteractions,
&bSumEkinhOld,
- CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
- checkNumberOfBondedInteractions(mdlog,
- cr,
- totalNumberOfBondedInteractions,
- top_global,
- &top,
- makeConstArrayRef(state->x),
- state->box,
- &shouldCheckNumberOfBondedInteractions);
+ cglo_flags,
+ step,
+ &observablesReducer);
}
clear_mat(force_vir);
bCalcEner = TRUE;
}
+ // bCalcEner is only here for when the last step is not a mulitple of nstfep
+ const bool computeDHDL = ((ir->efep != FreeEnergyPerturbationType::No || ir->bSimTemp)
+ && (do_per_step(step, nstfep) || bCalcEner));
+
/* 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));
- if (fr->useMts && !do_per_step(step, ir->nstfout))
+ | (bCalcEner ? GMX_FORCE_ENERGY : 0) | (computeDHDL ? GMX_FORCE_DHDL : 0));
+ if (simulationWork.useMts && !do_per_step(step, ir->nstfout))
{
+ // TODO: merge this with stepWork.useOnlyMtsCombinedForceBuffer
force_flags |= GMX_FORCE_DO_NOT_NEED_NORMAL_FORCE;
}
pull_work,
bNS,
force_flags,
- &top,
+ top,
constr,
enerd,
state->natoms,
&state->hist,
&f.view(),
force_vir,
- mdatoms,
+ *md,
+ fr->longRangeNonbondeds.get(),
nrnb,
wcycle,
shellfc,
step,
nrnb,
wcycle,
- &top,
+ top,
state->box,
state->x.arrayRefWithPadding(),
&state->hist,
&f.view(),
force_vir,
- mdatoms,
+ md,
enerd,
state->lambda,
fr,
mu_tot,
t,
ed ? ed->getLegacyED() : nullptr,
+ fr->longRangeNonbondeds.get(),
(bNS ? GMX_FORCE_NS : 0) | force_flags,
ddBalanceRegionHandler);
}
fr,
cr,
state,
- mdatoms,
- fcdata,
+ mdAtoms->mdatoms(),
+ &fcdata,
&MassQ,
&vcm,
- top_global,
- top,
enerd,
+ &observablesReducer,
ekind,
gstat,
&last_ekin,
bTrotter,
bExchanged,
&bSumEkinhOld,
- &shouldCheckNumberOfBondedInteractions,
&saved_conserved_quantity,
&f,
&upd,
&nullSignaller,
trotter_seq,
nrnb,
- mdlog,
fplog,
wcycle);
if (vsite != nullptr && needVirtualVelocitiesThisStep)
{
// Positions were calculated earlier
- wallcycle_start(wcycle, ewcVSITECONSTR);
+ wallcycle_start(wcycle, WallCycleCounter::VsiteConstr);
vsite->construct(state->x, state->v, state->box, VSiteOperation::Velocities);
- wallcycle_stop(wcycle, ewcVSITECONSTR);
+ wallcycle_stop(wcycle, WallCycleCounter::VsiteConstr);
}
}
state->dfhist,
step,
state->v.rvec_array(),
- mdatoms);
+ md->homenr,
+ md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
/* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */
if (MASTER(cr))
{
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);
+ bInteractiveMDstep = imdSession->run(step, bNS, state->box, state->x, t);
/* kludge -- virial is lost with restart for MTTK NPT control. Must reload (saved earlier). */
if (startingBehavior != StartingBehavior::NewSimulation && bFirstStep
if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */
{
gmx_bool bIfRandomize;
- bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state->v, &upd, constr);
+ bIfRandomize = update_randomize_velocities(ir,
+ step,
+ cr,
+ md->homenr,
+ md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ gmx::arrayRefFromArray(md->invmass, md->nr),
+ state->v,
+ &upd,
+ constr);
/* if we have constraints, we have to remove the kinetic energy parallel to the bonds */
if (constr && bIfRandomize)
{
if (!useGpuForUpdate)
{
- wallcycle_start(wcycle, ewcUPDATE);
+ wallcycle_start(wcycle, WallCycleCounter::Update);
}
/* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */
if (bTrotter)
{
- trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3);
+ trotter_update(ir,
+ step,
+ ekind,
+ enerd,
+ state,
+ total_vir,
+ md->homenr,
+ md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ gmx::arrayRefFromArray(md->invmass, md->nr),
+ &MassQ,
+ trotter_seq,
+ TrotterSequence::Three);
/* We can only do Berendsen coupling after we have summed
* the kinetic energy or virial. Since the happens
* in global_state after update, we should only do it at
}
else
{
- update_tcouple(step, ir, state, ekind, &MassQ, mdatoms);
+ update_tcouple(step,
+ ir,
+ state,
+ ekind,
+ &MassQ,
+ md->homenr,
+ md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
update_pcouple_before_coordinates(fplog, step, ir, state, pressureCouplingMu, M, bInitStep);
}
fr,
cr,
state,
- mdatoms,
- fcdata,
+ mdAtoms->mdatoms(),
+ &fcdata,
&MassQ,
&vcm,
pull_work,
enerd,
+ &observablesReducer,
ekind,
gstat,
&dvdl_constr,
{
if (useGpuForUpdate)
{
-
- wallcycle_stop(wcycle, ewcUPDATE);
-
- if (bNS && (bFirstStep || DOMAINDECOMP(cr)))
+ // On search steps, update handles to device vectors
+ if (bNS && (bFirstStep || haveDDAtomOrdering(*cr) || bExchanged))
{
integrator->set(stateGpu->getCoordinates(),
stateGpu->getVelocities(),
stateGpu->getForces(),
- top.idef,
- *mdatoms);
+ top->idef,
+ *md);
// Copy data to the GPU after buffers might have being reinitialized
+ /* The velocity copy is redundant if we had Center-of-Mass motion removed on
+ * the previous step. We don't check that now. */
stateGpu->copyVelocitiesToGpu(state->v, AtomLocality::Local);
- stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
+ if (bExchanged
+ || (!runScheduleWork->stepWork.haveGpuPmeOnThisRank
+ && !runScheduleWork->stepWork.useGpuXBufferOps))
+ {
+ stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
+ }
}
- if (simulationWork.useGpuPme && !runScheduleWork->simulationWork.useGpuPmePpCommunication
- && !thisRankHasDuty(cr, DUTY_PME))
- {
- // The PME forces were recieved to the host, so have to be copied
- stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::All);
- }
- else if (!runScheduleWork->stepWork.useGpuFBufferOps)
+ if ((simulationWork.useGpuPme && simulationWork.useCpuPmePpCommunication)
+ || (!runScheduleWork->stepWork.useGpuFBufferOps))
{
- // The buffer ops were not offloaded this step, so the forces are on the
- // host and have to be copied
+ // The PME forces were recieved to the host, and reduced on the CPU with the
+ // rest of the forces computed on the GPU, so the final forces have to be copied
+ // back to the GPU. Or the buffer ops were not offloaded this step, so the
+ // forces are on the host and have to be copied
stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::Local);
}
-
const bool doTemperatureScaling =
(ir->etc != TemperatureCoupling::No
&& do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
// 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,
- doTemperatureScaling,
- 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).
- // - Temperature is needed for the next step.
- if (bGStat || needHalfStepKineticEnergy)
- {
- stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
- stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
- }
+ integrator->integrate(stateGpu->getLocalForcesReadyOnDeviceEvent(
+ runScheduleWork->stepWork, runScheduleWork->simulationWork),
+ ir->delta_t,
+ true,
+ bCalcVir,
+ shake_vir,
+ doTemperatureScaling,
+ ekind->tcstat,
+ doParrinelloRahman,
+ ir->nstpcouple * ir->delta_t,
+ M);
}
else
{
* Using that acceleration would result in a virial with the slow
* force contribution would be a factor mtsFactor too large.
*/
- if (fr->useMts && bCalcVir && constr != nullptr)
+ if (simulationWork.useMts && bCalcVir && constr != nullptr)
{
- upd.update_for_constraint_virial(
- *ir, *mdatoms, *state, f.view().forceWithPadding(), *ekind);
+ upd.update_for_constraint_virial(*ir,
+ md->homenr,
+ md->havePartiallyFrozenAtoms,
+ gmx::arrayRefFromArray(md->invmass, md->nr),
+ gmx::arrayRefFromArray(md->invMassPerDim, md->nr),
+ *state,
+ f.view().forceWithPadding(),
+ *ekind);
constrain_coordinates(constr,
do_log,
}
ArrayRefWithPadding<const RVec> forceCombined =
- (fr->useMts && step % ir->mtsLevels[1].stepFactor == 0)
+ (simulationWork.useMts && step % ir->mtsLevels[1].stepFactor == 0)
? f.view().forceMtsCombinedWithPadding()
: f.view().forceWithPadding();
- upd.update_coords(
- *ir, step, mdatoms, state, forceCombined, fcdata, ekind, M, etrtPOSITION, cr, constr != nullptr);
+ upd.update_coords(*ir,
+ step,
+ md->homenr,
+ md->havePartiallyFrozenAtoms,
+ gmx::arrayRefFromArray(md->ptype, md->nr),
+ gmx::arrayRefFromArray(md->invmass, md->nr),
+ gmx::arrayRefFromArray(md->invMassPerDim, md->nr),
+ state,
+ forceCombined,
+ &fcdata,
+ ekind,
+ M,
+ etrtPOSITION,
+ cr,
+ constr != nullptr);
- wallcycle_stop(wcycle, ewcUPDATE);
+ wallcycle_stop(wcycle, WallCycleCounter::Update);
constrain_coordinates(constr,
do_log,
state,
upd.xp()->arrayRefWithPadding(),
&dvdl_constr,
- bCalcVir && !fr->useMts,
+ bCalcVir && !simulationWork.useMts,
shake_vir);
- upd.update_sd_second_half(
- *ir, step, &dvdl_constr, mdatoms, state, cr, nrnb, wcycle, constr, do_log, do_ene);
- upd.finish_update(*ir, mdatoms, state, wcycle, constr != nullptr);
+ upd.update_sd_second_half(*ir,
+ step,
+ &dvdl_constr,
+ md->homenr,
+ gmx::arrayRefFromArray(md->ptype, md->nr),
+ gmx::arrayRefFromArray(md->invmass, md->nr),
+ state,
+ cr,
+ nrnb,
+ wcycle,
+ constr,
+ do_log,
+ do_ene);
+ upd.finish_update(
+ *ir, md->havePartiallyFrozenAtoms, md->homenr, state, wcycle, constr != nullptr);
}
if (ir->bPull && ir->pull->bSetPbcRefToPrevStepCOM)
{
- updatePrevStepPullCom(pull_work, state);
+ updatePrevStepPullCom(pull_work, state->pull_com_prev_step);
}
enerd->term[F_DVDL_CONSTR] += dvdl_constr;
// and when algorithms require it.
const bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
- if (bGStat || needHalfStepKineticEnergy || doInterSimSignal)
+ if (useGpuForUpdate)
{
- // Copy coordinates when needed to stop the CM motion.
- if (useGpuForUpdate && !EI_VV(ir->eI) && bStopCM)
+ const bool coordinatesRequiredForStopCM =
+ bStopCM && (bGStat || needHalfStepKineticEnergy || doInterSimSignal)
+ && !EI_VV(ir->eI);
+
+ // Copy coordinates when needed to stop the CM motion or for replica exchange
+ if (coordinatesRequiredForStopCM || bDoReplEx)
{
stateGpu->copyCoordinatesFromGpu(state->x, AtomLocality::Local);
stateGpu->waitCoordinatesReadyOnHost(AtomLocality::Local);
}
+
+ // Copy velocities back to the host if:
+ // - Globals are computed this step (includes the energy output steps).
+ // - Temperature is needed for the next step.
+ // - This is a replica exchange step (even though we will only need
+ // the velocities if an exchange succeeds)
+ if (bGStat || needHalfStepKineticEnergy || bDoReplEx)
+ {
+ stateGpu->copyVelocitiesFromGpu(state->v, AtomLocality::Local);
+ stateGpu->waitVelocitiesReadyOnHost(AtomLocality::Local);
+ }
+ }
+
+ if (bGStat || needHalfStepKineticEnergy || doInterSimSignal)
+ {
// Since we're already communicating at this step, we
// can propagate intra-simulation signals. Note that
// check_nstglobalcomm has the responsibility for
bool doIntraSimSignal = true;
SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
- compute_globals(
- gstat,
- cr,
- ir,
- fr,
- ekind,
- makeConstArrayRef(state->x),
- makeConstArrayRef(state->v),
- state->box,
- mdatoms,
- nrnb,
- &vcm,
- wcycle,
- enerd,
- force_vir,
- shake_vir,
- total_vir,
- pres,
- (!EI_VV(ir->eI) && bCalcEner && constr != nullptr) ? constr->rmsdData()
- : gmx::ArrayRef<real>{},
- &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));
- checkNumberOfBondedInteractions(mdlog,
- cr,
- totalNumberOfBondedInteractions,
- top_global,
- &top,
- makeConstArrayRef(state->x),
- state->box,
- &shouldCheckNumberOfBondedInteractions);
+ compute_globals(gstat,
+ cr,
+ ir,
+ fr,
+ ekind,
+ makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v),
+ state->box,
+ md,
+ nrnb,
+ &vcm,
+ wcycle,
+ enerd,
+ force_vir,
+ shake_vir,
+ total_vir,
+ pres,
+ &signaller,
+ lastbox,
+ &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,
+ step,
+ &observablesReducer);
if (!EI_VV(ir->eI) && bStopCM)
{
process_and_stopcm_grp(
- fplog, &vcm, *mdatoms, makeArrayRef(state->x), makeArrayRef(state->v));
- inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
+ fplog, &vcm, *md, makeArrayRef(state->x), makeArrayRef(state->v));
+ inc_nrnb(nrnb, eNR_STOPCM, md->homenr);
// TODO: The special case of removing CM motion should be dealt more gracefully
if (useGpuForUpdate)
{
+ // Issue #3988, #4106.
+ stateGpu->resetCoordinatesCopiedToDeviceEvent(AtomLocality::Local);
stateGpu->copyCoordinatesToGpu(state->x, AtomLocality::Local);
// Here we block until the H2D copy completes because event sync with the
// force kernels that use the coordinates on the next steps is not implemented
// (not because of a race on state->x being modified on the CPU while H2D is in progress).
stateGpu->waitCoordinatesCopiedToDevice(AtomLocality::Local);
// If the COM removal changed the velocities on the CPU, this has to be accounted for.
- if (vcm.mode != ecmNO)
+ if (vcm.mode != ComRemovalAlgorithm::No)
{
stateGpu->copyVelocitiesToGpu(state->v, AtomLocality::Local);
}
but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
- if (ir->efep != efepNO && !EI_VV(ir->eI))
+ if (ir->efep != FreeEnergyPerturbationType::No && !EI_VV(ir->eI))
{
/* Sum up the foreign energy and dK/dl terms for md and sd.
Currently done every step so that dH/dl is correct in the .edr */
accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
}
- update_pcouple_after_coordinates(
- fplog, step, ir, mdatoms, pres, force_vir, shake_vir, pressureCouplingMu, state, nrnb, upd.deform(), !useGpuForUpdate);
+ bool scaleCoordinates = !useGpuForUpdate || bDoReplEx;
+ update_pcouple_after_coordinates(fplog,
+ step,
+ ir,
+ md->homenr,
+ md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ pres,
+ force_vir,
+ shake_vir,
+ pressureCouplingMu,
+ state,
+ nrnb,
+ upd.deform(),
+ scaleCoordinates);
const bool doBerendsenPressureCoupling = (inputrec->epc == PressureCoupling::Berendsen
&& do_per_step(step, inputrec->nstpcouple));
/* ######### BEGIN PREPARING EDR OUTPUT ########### */
/* use the directly determined last velocity, not actually the averaged half steps */
- if (bTrotter && ir->eI == eiVV)
+ if (bTrotter && ir->eI == IntegrationAlgorithm::VV)
{
enerd->term[F_EKIN] = last_ekin;
}
{
/* only needed if doing expanded ensemble */
PrintFreeEnergyInfoToFile(fplog,
- ir->fepvals,
- ir->expandedvals,
- ir->bSimTemp ? ir->simtempvals : nullptr,
+ ir->fepvals.get(),
+ ir->expandedvals.get(),
+ ir->bSimTemp ? ir->simtempvals.get() : nullptr,
state_global->dfhist,
state->fep_state,
ir->nstlog,
}
if (bCalcEner)
{
- energyOutput.addDataAtEnergyStep(bDoDHDL,
+ const bool outputDHDL = (computeDHDL && do_per_step(step, ir->fepvals->nstdhdl));
+
+ energyOutput.addDataAtEnergyStep(outputDHDL,
bCalcEnerStep,
t,
- mdatoms->tmass,
+ md->tmass,
enerd,
- ir->fepvals,
- ir->expandedvals,
+ ir->fepvals.get(),
lastbox,
PTCouplingArrays{ state->boxv,
state->nosehoover_xi,
state->nhpres_xi,
state->nhpres_vxi },
state->fep_state,
- shake_vir,
- force_vir,
total_vir,
pres,
ekind,
* 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 != SwapType::No) && (step > 0) && !bLastStep
+ && do_per_step(step, ir->swap->nstswap))
{
bNeedRepartition = do_swapcoords(cr,
step,
MASTER(cr) && mdrunOptions.verbose,
bRerunMD);
- if (bNeedRepartition && DOMAINDECOMP(cr))
+ if (bNeedRepartition && haveDDAtomOrdering(*cr))
{
dd_collect_state(cr->dd, state, state_global);
}
bExchanged = replica_exchange(fplog, cr, ms, repl_ex, state_global, enerd, state, step, t);
}
- if ((bExchanged || bNeedRepartition) && DOMAINDECOMP(cr))
+ if ((bExchanged || bNeedRepartition) && haveDDAtomOrdering(*cr))
{
dd_partition_system(fplog,
mdlog,
TRUE,
1,
state_global,
- *top_global,
- ir,
+ top_global,
+ *ir,
imdSession,
pull_work,
state,
&f,
mdAtoms,
- &top,
+ top,
fr,
vsite,
constr,
nrnb,
wcycle,
FALSE);
- shouldCheckNumberOfBondedInteractions = true;
- upd.setNumAtoms(state->natoms);
+ upd.updateAfterPartition(state->natoms,
+ md->cFREEZE ? gmx::arrayRefFromArray(md->cFREEZE, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cTC ? gmx::arrayRefFromArray(md->cTC, md->nr)
+ : gmx::ArrayRef<const unsigned short>(),
+ md->cACC ? gmx::arrayRefFromArray(md->cACC, md->nr)
+ : gmx::ArrayRef<const unsigned short>());
+ fr->longRangeNonbondeds->updateAfterPartition(*md);
}
bFirstStep = FALSE;
/* 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))
+ if ((state->flags & enumValueToBitMask(StateEntry::PressurePrevious))
&& (bGStatEveryStep || (ir->nstpcouple > 0 && step % ir->nstpcouple == 0)))
{
/* Store the pressure in t_state for pressure coupling
rescale_membed(step_rel, membed, as_rvec_array(state_global->x.data()));
}
- cycles = wallcycle_stop(wcycle, ewcSTEP);
- if (DOMAINDECOMP(cr) && wcycle)
+ cycles = wallcycle_stop(wcycle, WallCycleCounter::Step);
+ if (haveDDAtomOrdering(*cr) && wcycle)
{
dd_cycles_add(cr->dd, cycles, ddCyclStep);
}
/* increase the MD step number */
step++;
step_rel++;
+ observablesReducer.markAsReadyToReduce();
#if GMX_FAHCORE
if (MASTER(cr))
/* Stop measuring walltime */
walltime_accounting_end_time(walltime_accounting);
- if (!thisRankHasDuty(cr, DUTY_PME))
+ if (simulationWork.haveSeparatePmeRank)
{
/* Tell the PME only node to finish */
gmx_pme_send_finish(cr);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff