#include <algorithm>
#include <memory>
+#include <numeric>
-#include "gromacs/awh/awh.h"
+#include "gromacs/applied_forces/awh/awh.h"
#include "gromacs/commandline/filenm.h"
#include "gromacs/domdec/collect.h"
#include "gromacs/domdec/dlbtiming.h"
#include "gromacs/domdec/domdec.h"
#include "gromacs/domdec/domdec_network.h"
#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/domdec/gpuhaloexchange.h"
#include "gromacs/domdec/mdsetup.h"
#include "gromacs/domdec/partition.h"
#include "gromacs/essentialdynamics/edsam.h"
-#include "gromacs/ewald/pme.h"
#include "gromacs/ewald/pme_load_balancing.h"
+#include "gromacs/ewald/pme_pp.h"
#include "gromacs/fileio/trxio.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/gpu_utils/device_stream_manager.h"
#include "gromacs/gpu_utils/gpu_utils.h"
#include "gromacs/imd/imd.h"
-#include "gromacs/listed_forces/manage_threading.h"
+#include "gromacs/listed_forces/listed_forces.h"
#include "gromacs/math/functions.h"
-#include "gromacs/math/utilities.h"
+#include "gromacs/math/invertmatrix.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/mdlib/checkpointhandler.h"
#include "gromacs/mdlib/compute_io.h"
#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/coupling.h"
#include "gromacs/mdlib/ebin.h"
#include "gromacs/mdlib/enerdata_utils.h"
#include "gromacs/mdlib/energyoutput.h"
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/force_flags.h"
#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/freeenergyparameters.h"
#include "gromacs/mdlib/md_support.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
#include "gromacs/mdlib/tgroup.h"
#include "gromacs/mdlib/trajectory_writing.h"
#include "gromacs/mdlib/update.h"
-#include "gromacs/mdlib/update_constrain_cuda.h"
+#include "gromacs/mdlib/update_constrain_gpu.h"
#include "gromacs/mdlib/vcm.h"
#include "gromacs/mdlib/vsite.h"
#include "gromacs/mdrunutility/handlerestart.h"
#include "gromacs/mdtypes/df_history.h"
#include "gromacs/mdtypes/energyhistory.h"
#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forcebuffers.h"
#include "gromacs/mdtypes/forcerec.h"
#include "gromacs/mdtypes/group.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/mdrunoptions.h"
+#include "gromacs/mdtypes/multipletimestepping.h"
#include "gromacs/mdtypes/observableshistory.h"
#include "gromacs/mdtypes/pullhistory.h"
#include "gromacs/mdtypes/simulation_workload.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/mdtypes/state_propagator_data_gpu.h"
-#include "gromacs/modularsimulator/energyelement.h"
+#include "gromacs/modularsimulator/energydata.h"
#include "gromacs/nbnxm/gpu_data_mgmt.h"
#include "gromacs/nbnxm/nbnxm.h"
-#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/output.h"
#include "gromacs/pulling/pull.h"
// will go away eventually.
t_inputrec* ir = inputrec;
int64_t step, step_rel;
- double t, t0 = ir->init_t, lam0[efptNR];
+ double t, t0 = ir->init_t;
gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
- gmx_bool bNS, bNStList, bStopCM, bFirstStep, bInitStep, bLastStep = FALSE;
+ gmx_bool bNS = FALSE, 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 pressureCouplingMu, M;
- gmx_repl_ex_t repl_ex = nullptr;
- gmx_localtop_t top;
- PaddedHostVector<gmx::RVec> f{};
- gmx_global_stat_t gstat;
- 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;
+ tensor force_vir = { { 0 } }, shake_vir = { { 0 } }, total_vir = { { 0 } }, pres = { { 0 } };
+ int i, m;
+ rvec mu_tot;
+ matrix pressureCouplingMu, M;
+ gmx_repl_ex_t repl_ex = nullptr;
+ gmx_global_stat_t gstat;
+ 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;
/* for FEP */
int nstfep = 0;
double cycles;
int nstglobalcomm = computeGlobalCommunicationPeriod(mdlog, ir, cr);
bGStatEveryStep = (nstglobalcomm == 1);
- SimulationGroups* groups = &top_global->groups;
+ const SimulationGroups* groups = &top_global->groups;
std::unique_ptr<EssentialDynamics> ed = nullptr;
if (opt2bSet("-ei", nfile, fnm))
"Either specify the -ei option to mdrun, or do not use this checkpoint file.");
}
- initialize_lambdas(fplog, *ir, MASTER(cr), &state_global->fep_state, state_global->lambda, lam0);
- Update upd(ir, deform);
+ 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);
+ Update upd(*ir, deform);
const bool doSimulatedAnnealing = initSimulatedAnnealing(ir, &upd);
const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
+ const 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 =
- (fcd->disres.nsystems > 1) || ((ms != nullptr) && (fcd->orires.nr != 0));
+ (fcdata.disres->nsystems > 1) || ((ms != nullptr) && (fcdata.orires->nr != 0));
bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
// Replica exchange, ensemble restraints and AWH need all
t_state* state;
+ gmx_localtop_t top(top_global->ffparams);
+
auto mdatoms = mdAtoms->mdatoms();
- std::unique_ptr<UpdateConstrainCuda> integrator;
+ 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;
+ ForceBuffers f(fr->useMts, ((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
+ ? PinningPolicy::PinnedIfSupported
+ : PinningPolicy::CannotBePinned);
if (DOMAINDECOMP(cr))
{
- dd_init_local_top(*top_global, &top);
-
stateInstance = std::make_unique<t_state>();
state = stateInstance.get();
dd_init_local_state(cr->dd, state_global, state);
else
{
state_change_natoms(state_global, state_global->natoms);
- f.resizeWithPadding(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, &graph, mdAtoms, constr, vsite, shellfc);
+ mdAlgorithmsSetupAtomData(cr, ir, *top_global, &top, fr, &f, mdAtoms, constr, vsite, shellfc);
upd.setNumAtoms(state->natoms);
}
- 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;
+ std::unique_ptr<UpdateConstrainGpu> integrator;
StatePropagatorDataGpu* stateGpu = fr->stateGpu;
+ // TODO: the assertions below should be handled by UpdateConstraintsBuilder.
if (useGpuForUpdate)
{
GMX_RELEASE_ASSERT(!DOMAINDECOMP(cr) || ddUsesUpdateGroups(*cr->dd) || constr == nullptr
GMX_RELEASE_ASSERT(
ir->etc != etcNOSEHOOVER,
"Nose-Hoover temperature coupling is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(ir->epc == epcNO || ir->epc == epcPARRINELLORAHMAN || ir->epc == epcBERENDSEN,
- "Only Parrinello-Rahman and Berendsen pressure coupling are supported "
- "with the GPU update.\n");
+ GMX_RELEASE_ASSERT(
+ ir->epc == epcNO || ir->epc == epcPARRINELLORAHMAN || ir->epc == epcBERENDSEN
+ || ir->epc == epcCRESCALE,
+ "Only Parrinello-Rahman, Berendsen, and C-rescale pressure coupling are 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 || !pull_have_constraint(ir->pull),
"Constraints pulling is not supported with the GPU update.\n");
- GMX_RELEASE_ASSERT(fcd->orires.nr == 0,
+ GMX_RELEASE_ASSERT(fcdata.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.");
+ GMX_RELEASE_ASSERT(
+ ir->efep == efepNO
+ || (!haveFreeEnergyType(*ir, efptBONDED) && !haveFreeEnergyType(*ir, efptMASS)),
+ "Free energy perturbation of masses and constraints are not supported with the GPU "
+ "update.");
if (constr != nullptr && constr->numConstraintsTotal() > 0)
{
{
GMX_LOG(mdlog.info).asParagraph().appendText("Updating coordinates on the GPU.");
}
- integrator = std::make_unique<UpdateConstrainCuda>(
- *ir, *top_global, stateGpu->getUpdateStream(), stateGpu->xUpdatedOnDevice());
-
- t_pbc pbc;
- set_pbc(&pbc, epbcXYZ, state->box);
- integrator->setPbc(&pbc);
+ GMX_RELEASE_ASSERT(fr->deviceStreamManager != nullptr,
+ "Device stream manager should be initialized in order to use GPU "
+ "update-constraints.");
+ GMX_RELEASE_ASSERT(
+ fr->deviceStreamManager->streamIsValid(gmx::DeviceStreamType::UpdateAndConstraints),
+ "Update stream should be initialized in order to use GPU "
+ "update-constraints.");
+ integrator = std::make_unique<UpdateConstrainGpu>(
+ *ir, *top_global, fr->deviceStreamManager->context(),
+ fr->deviceStreamManager->stream(gmx::DeviceStreamType::UpdateAndConstraints),
+ stateGpu->xUpdatedOnDevice());
+
+ integrator->setPbc(PbcType::Xyz, state->box);
}
if (useGpuForPme || (useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
{
changePinningPolicy(&state->x, PinningPolicy::PinnedIfSupported);
}
- if ((useGpuForNonbonded && useGpuForBufferOps) || useGpuForUpdate)
- {
- changePinningPolicy(&f, PinningPolicy::PinnedIfSupported);
- }
if (useGpuForUpdate)
{
changePinningPolicy(&state->v, PinningPolicy::PinnedIfSupported);
if (MASTER(cr))
{
- EnergyElement::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput);
+ EnergyData::initializeEnergyHistory(startingBehavior, observablesHistory, &energyOutput);
}
- preparePrevStepPullCom(ir, pull_work, mdatoms, state, state_global, cr,
+ preparePrevStepPullCom(ir, pull_work, mdatoms->massT, state, state_global, cr,
startingBehavior != StartingBehavior::NewSimulation);
// TODO: Remove this by converting AWH into a ForceProvider
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->nr, mdatoms->homenr,
+ 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);
+ vsite->construct(state->x, ir->delta_t, {}, state->box);
}
}
nstfep = ir->fepvals->nstdhdl;
if (ir->bExpanded)
{
- nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep);
+ nstfep = std::gcd(ir->expandedvals->nstexpanded, nstfep);
}
if (useReplicaExchange)
{
- nstfep = gmx_greatest_common_divisor(replExParams.exchangeInterval, nstfep);
+ nstfep = std::gcd(replExParams.exchangeInterval, nstfep);
+ }
+ if (ir->bDoAwh)
+ {
+ nstfep = std::gcd(ir->awhParams->nstSampleCoord, nstfep);
}
}
bool hasReadEkinState = MASTER(cr) ? state_global->ekinstate.hasReadEkinState : false;
if (PAR(cr))
{
- gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr);
+ gmx_bcast(sizeof(hasReadEkinState), &hasReadEkinState, cr->mpi_comm_mygroup);
}
if (hasReadEkinState)
{
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,
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, nullptr,
+ enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
cglo_flags_iteration
| (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
/* At initialization, do not pass x with acceleration-correction mode
* to avoid (incorrect) correction of the initial coordinates.
*/
- rvec* xPtr = nullptr;
- if (vcm.mode != ecmLINEAR_ACCELERATION_CORRECTION)
- {
- xPtr = state->x.rvec_array();
- }
- process_and_stopcm_grp(fplog, &vcm, *mdatoms, xPtr, state->v.rvec_array());
+ 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);
}
}
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global, &top,
- state->x.rvec_array(), state->box,
+ makeConstArrayRef(state->x), 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,
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, nullptr,
+ enerd, force_vir, shake_vir, total_vir, pres, constr, &nullSignaller,
state->box, nullptr, &bSumEkinhOld, cglo_flags & ~CGLO_PRESSURE);
}
bExchanged = FALSE;
bNeedRepartition = FALSE;
+ step = ir->init_step;
+ step_rel = 0;
+
auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), simulationsShareState,
MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstSignalComm,
const DDBalanceRegionHandler ddBalanceRegionHandler(cr);
- step = ir->init_step;
- step_rel = 0;
-
- // TODO extract this to new multi-simulation module
if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
{
- 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!");
- }
- if (!multisim_int_all_are_equal(ms, ir->init_step))
- {
- if (simulationsShareState)
- {
- if (MASTER(cr))
- {
- gmx_fatal(FARGS,
- "The initial step is not consistent across multi simulations which "
- "share the state");
- }
- gmx_barrier(cr);
- }
- else
- {
- GMX_LOG(mdlog.warning)
- .appendText(
- "Note: The initial step is not consistent across multi "
- "simulations,\n"
- "but we are proceeding anyway!");
- }
- }
+ logInitialMultisimStatus(ms, cr, mdlog, simulationsShareState, ir->nsteps, ir->init_step);
}
/* and stop now if we should */
if (ir->efep != efepNO || ir->bSimTemp)
{
/* find and set the current lambdas */
- setCurrentLambdasLocal(step, ir->fepvals, lam0, state->lambda, state->fep_state);
+ 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));
/* Correct the new box if it is too skewed */
if (inputrecDynamicBox(ir))
{
- if (correct_box(fplog, step, state->box, graph))
+ if (correct_box(fplog, step, state->box))
{
bMasterState = TRUE;
// If update is offloaded, it should be informed about the box size change
if (useGpuForUpdate)
{
- t_pbc pbc;
- set_pbc(&pbc, epbcXYZ, state->box);
- integrator->setPbc(&pbc);
+ integrator->setPbc(PbcType::Xyz, state->box);
}
}
}
}
}
+ // Allocate or re-size GPU halo exchange object, if necessary
+ if (bNS && havePPDomainDecomposition(cr) && simulationWork.useGpuHaloExchange)
+ {
+ GMX_RELEASE_ASSERT(fr->deviceStreamManager != nullptr,
+ "GPU device manager has to be initialized to use GPU "
+ "version of halo exchange.");
+ constructGpuHaloExchange(mdlog, *cr, *fr->deviceStreamManager, wcycle);
+ }
+
if (MASTER(cr) && do_log)
{
- energyOutput.printHeader(fplog, step, t); /* can we improve the information printed here? */
+ gmx::EnergyOutput::printHeader(fplog, step,
+ t); /* can we improve the information printed here? */
}
if (ir->efep != efepNO)
/* 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,
+ compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
+ enerd, nullptr, nullptr, nullptr, nullptr, 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,
+ &top, makeConstArrayRef(state->x), state->box,
&shouldCheckNumberOfBondedInteractions);
}
clear_mat(force_vir);
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))
+ {
+ force_flags |= GMX_FORCE_DO_NOT_NEED_NORMAL_FORCE;
+ }
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,
+ imdSession, pull_work, bNS, force_flags, &top, constr, enerd,
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);
+ state->v.arrayRefWithPadding(), state->box, state->lambda,
+ &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc,
+ fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler);
}
else
{
*/
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,
+ &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork,
+ vsite, mu_tot, t, ed ? ed->getLegacyED() : nullptr,
(bNS ? GMX_FORCE_NS : 0) | force_flags, ddBalanceRegionHandler);
}
trotter_seq, ettTSEQ1);
}
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
- etrtVELOCITY1, cr, constr);
+ upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
+ M, etrtVELOCITY1, cr, constr != nullptr);
wallcycle_stop(wcycle, ewcUPDATE);
- constrain_velocities(step, nullptr, state, shake_vir, constr, bCalcVir, do_log, do_ene);
+ constrain_velocities(constr, do_log, do_ene, step, state, nullptr, bCalcVir, shake_vir);
wallcycle_start(wcycle, ewcUPDATE);
/* if VV, compute the pressure and constraints */
/* For VV2, we strictly only need this if using pressure
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,
+ 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, constr, &nullSignaller,
state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
(bGStat ? CGLO_GSTAT : 0) | (bCalcEner ? CGLO_ENERGY : 0)
| (bTemp ? CGLO_TEMPERATURE : 0) | (bPres ? CGLO_PRESSURE : 0)
b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in
EkinAveVel because it's needed for the pressure */
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
- top_global, &top, state->x.rvec_array(), state->box,
- &shouldCheckNumberOfBondedInteractions);
+ top_global, &top, makeConstArrayRef(state->x),
+ state->box, &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, makeArrayRef(state->x),
+ makeArrayRef(state->v));
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
}
wallcycle_start(wcycle, ewcUPDATE);
/* 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, makeConstArrayRef(state->x),
+ makeConstArrayRef(state->v), state->box, mdatoms, nrnb, &vcm, wcycle,
+ enerd, nullptr, nullptr, nullptr, nullptr, constr, &nullSignaller,
+ state->box, nullptr, &bSumEkinhOld, CGLO_GSTAT | CGLO_TEMPERATURE);
wallcycle_start(wcycle, ewcUPDATE);
}
}
{
saved_conserved_quantity -= enerd->term[F_DISPCORR];
}
- /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */
+ /* sum up the foreign kinetic energy and dK/dl terms for vv. currently done every step so that dhdl is correct in the .edr */
if (ir->efep != efepNO)
{
- sum_dhdl(enerd, state->lambda, *ir->fepvals);
+ accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
}
}
if (runScheduleWork->stepWork.useGpuFBufferOps && (simulationWork.useGpuUpdate && !vsite)
&& do_per_step(step, ir->nstfout))
{
- stateGpu->copyForcesFromGpu(ArrayRef<RVec>(f), AtomLocality::Local);
+ stateGpu->copyForcesFromGpu(f.view().force(), AtomLocality::Local);
stateGpu->waitForcesReadyOnHost(AtomLocality::Local);
}
/* Now we have the energies and forces corresponding to the
* the update.
*/
do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state, state_global,
- observablesHistory, top_global, fr, outf, energyOutput, ekind, f,
- checkpointHandler->isCheckpointingStep(), bRerunMD, bLastStep,
- mdrunOptions.writeConfout, bSumEkinhOld);
+ observablesHistory, top_global, fr, outf, energyOutput, ekind,
+ f.view().force(), 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);
/* 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(constr, do_log, do_ene, step, state, nullptr, false, nullptr);
}
}
/* Box is changed in update() when we do pressure coupling,
if (EI_VV(ir->eI))
{
/* velocity half-step update */
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
- etrtVELOCITY2, cr, constr);
+ upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
+ M, etrtVELOCITY2, cr, constr != nullptr);
}
/* Above, initialize just copies ekinh into ekin,
// If the buffer ops were not offloaded this step, the forces are on the host and have to be copied
if (!runScheduleWork->stepWork.useGpuFBufferOps)
{
- stateGpu->copyForcesToGpu(ArrayRef<RVec>(f), AtomLocality::Local);
+ stateGpu->copyForcesToGpu(f.view().force(), AtomLocality::Local);
}
const bool doTemperatureScaling =
}
else
{
- update_coords(step, ir, mdatoms, state, f.arrayRefWithPadding(), fcd, ekind, M, &upd,
- etrtPOSITION, cr, constr);
+ /* With multiple time stepping we need to do an additional normal
+ * update step to obtain the virial, as the actual MTS integration
+ * using an acceleration where the slow forces are multiplied by mtsFactor.
+ * 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)
+ {
+ upd.update_for_constraint_virial(*ir, *mdatoms, *state, f.view().forceWithPadding(), *ekind);
+
+ constrain_coordinates(constr, do_log, do_ene, step, state,
+ upd.xp()->arrayRefWithPadding(), &dvdl_constr, bCalcVir, shake_vir);
+ }
+
+ ArrayRefWithPadding<const RVec> forceCombined =
+ (fr->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);
wallcycle_stop(wcycle, ewcUPDATE);
- constrain_coordinates(step, &dvdl_constr, state, shake_vir, &upd, constr, bCalcVir,
- do_log, do_ene);
+ constrain_coordinates(constr, do_log, do_ene, step, state, upd.xp()->arrayRefWithPadding(),
+ &dvdl_constr, bCalcVir && !fr->useMts, shake_vir);
- 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);
+ 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);
}
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,
+ 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, 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);
+ upd.update_coords(*ir, step, mdatoms, state, f.view().forceWithPadding(), fcdata, ekind,
+ M, etrtPOSITION, cr, constr != nullptr);
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);
+ upd.finish_update(*ir, mdatoms, state, wcycle, false);
}
if (EI_VV(ir->eI))
{
/* this factor or 2 correction is necessary
because half of the constraint force is removed
in the vv step, so we have to double it. See
- the Redmine issue #1255. It is not yet clear
+ the Issue #1255. It is not yet clear
if the factor of 2 is exact, or just a very
good approximation, and this will be
investigated. The next step is to see if this
if (vsite != nullptr)
{
wallcycle_start(wcycle, ewcVSITECONSTR);
- if (graph != nullptr)
- {
- 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);
-
- if (graph != nullptr)
- {
- unshift_self(graph, state->box, state->x.rvec_array());
- }
+ vsite->construct(state->x, ir->delta_t, state->v, state->box);
wallcycle_stop(wcycle, ewcVSITECONSTR);
}
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)
- | (!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));
+ 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, 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));
checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
- top_global, &top, state->x.rvec_array(), state->box,
- &shouldCheckNumberOfBondedInteractions);
+ top_global, &top, makeConstArrayRef(state->x),
+ 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, makeArrayRef(state->x),
+ makeArrayRef(state->v));
inc_nrnb(nrnb, eNR_STOPCM, mdatoms->homenr);
// TODO: The special case of removing CM motion should be dealt more gracefully
if (ir->efep != efepNO && !EI_VV(ir->eI))
{
- /* Sum up the foreign energy and dhdl terms for md and sd.
- Currently done every step so that dhdl is correct in the .edr */
- sum_dhdl(enerd, state->lambda, *ir->fepvals);
+ /* 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, !useGpuForUpdate);
+ pressureCouplingMu, state, nrnb, upd.deform(), !useGpuForUpdate);
const bool doBerendsenPressureCoupling =
(inputrec->epc == epcBERENDSEN && do_per_step(step, inputrec->nstpcouple));
- if (useGpuForUpdate && (doBerendsenPressureCoupling || doParrinelloRahman))
+ const bool doCRescalePressureCoupling =
+ (inputrec->epc == epcCRESCALE && do_per_step(step, inputrec->nstpcouple));
+ if (useGpuForUpdate
+ && (doBerendsenPressureCoupling || doCRescalePressureCoupling || doParrinelloRahman))
{
integrator->scaleCoordinates(pressureCouplingMu);
- t_pbc pbc;
- set_pbc(&pbc, epbcXYZ, state->box);
- integrator->setPbc(&pbc);
+ if (doCRescalePressureCoupling)
+ {
+ matrix pressureCouplingInvMu;
+ gmx::invertBoxMatrix(pressureCouplingMu, pressureCouplingInvMu);
+ integrator->scaleVelocities(pressureCouplingInvMu);
+ }
+ integrator->setPbc(PbcType::Xyz, state->box);
}
/* ################# END UPDATE STEP 2 ################# */
}
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, ir->fepvals,
+ ir->expandedvals, lastbox,
+ PTCouplingArrays{ state->boxv, state->nosehoover_xi, state->nosehoover_vxi,
+ state->nhpres_xi, state->nhpres_vxi },
+ state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr);
}
else
{
if (doSimulatedAnnealing)
{
- energyOutput.printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
+ gmx::EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups,
+ &(ir->opts));
}
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,
+ fr->fcdata.get(), awh.get());
+ }
+ if (do_log && ir->bDoAwh && awh->hasFepLambdaDimension())
+ {
+ const bool isInitialOutput = false;
+ printLambdaStateToLog(fplog, state->lambda, isInitialOutput);
}
if (ir->bPull)
/* Gets written into the state at the beginning of next loop*/
state->fep_state = lamnew;
}
+ else if (ir->bDoAwh && awh->needForeignEnergyDifferences(step))
+ {
+ state->fep_state = awh->fepLambdaState();
+ }
/* Print the remaining wall clock time for the run */
if (isMasterSimMasterRank(ms, MASTER(cr)) && (do_verbose || gmx_got_usr_signal()) && !bPMETunePrinting)
{
{
if (ir->nstcalcenergy > 0)
{
- energyOutput.printAnnealingTemperatures(fplog, groups, &(ir->opts));
+ gmx::EnergyOutput::printAnnealingTemperatures(fplog, groups, &(ir->opts));
energyOutput.printAverages(fplog, groups);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff