class Constraints::Impl
{
public:
- Impl(const gmx_mtop_t& mtop_p,
- const t_inputrec& ir_p,
- pull_t* pull_work,
- FILE* log_p,
- const t_commrec* cr_p,
- bool useUpdateGroups,
- const gmx_multisim_t* ms,
- t_nrnb* nrnb,
- gmx_wallcycle* wcycle_p,
- bool pbcHandlingRequired,
- int numConstraints,
- int numSettles);
+ Impl(const gmx_mtop_t& mtop_p,
+ const t_inputrec& ir_p,
+ pull_t* pull_work,
+ FILE* log_p,
+ const t_commrec* cr_p,
+ bool useUpdateGroups,
+ const gmx_multisim_t* ms,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle_p,
+ bool pbcHandlingRequired,
+ ObservablesReducerBuilder* observablesReducerBuilder,
+ int numConstraints,
+ int numSettles);
~Impl();
void setConstraints(gmx_localtop_t* top,
int numAtoms,
return bOK;
} // namespace gmx
-ArrayRef<real> Constraints::rmsdData() const
-{
- if (impl_->lincsd)
- {
- return lincs_rmsdData(impl_->lincsd);
- }
- else
- {
- return {};
- }
-}
-
real Constraints::rmsd() const
{
if (impl_->lincsd)
return mapping;
}
-Constraints::Constraints(const gmx_mtop_t& mtop,
- const t_inputrec& ir,
- pull_t* pull_work,
- FILE* log,
- const t_commrec* cr,
- const bool useUpdateGroups,
- const gmx_multisim_t* ms,
- t_nrnb* nrnb,
- gmx_wallcycle* wcycle,
- bool pbcHandlingRequired,
- int numConstraints,
- int numSettles) :
- impl_(new Impl(mtop, ir, pull_work, log, cr, useUpdateGroups, ms, nrnb, wcycle, pbcHandlingRequired, numConstraints, numSettles))
+Constraints::Constraints(const gmx_mtop_t& mtop,
+ const t_inputrec& ir,
+ pull_t* pull_work,
+ FILE* log,
+ const t_commrec* cr,
+ const bool useUpdateGroups,
+ const gmx_multisim_t* ms,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle,
+ bool pbcHandlingRequired,
+ ObservablesReducerBuilder* observablesReducerBuilder,
+ int numConstraints,
+ int numSettles) :
+ impl_(new Impl(mtop, ir, pull_work, log, cr, useUpdateGroups, ms, nrnb, wcycle, pbcHandlingRequired, observablesReducerBuilder, numConstraints, numSettles))
{
}
-Constraints::Impl::Impl(const gmx_mtop_t& mtop_p,
- const t_inputrec& ir_p,
- pull_t* pull_work,
- FILE* log_p,
- const t_commrec* cr_p,
- const bool useUpdateGroups,
- const gmx_multisim_t* ms_p,
- t_nrnb* nrnb_p,
- gmx_wallcycle* wcycle_p,
- bool pbcHandlingRequired,
- int numConstraints,
- int numSettles) :
+Constraints::Impl::Impl(const gmx_mtop_t& mtop_p,
+ const t_inputrec& ir_p,
+ pull_t* pull_work,
+ FILE* log_p,
+ const t_commrec* cr_p,
+ const bool useUpdateGroups,
+ const gmx_multisim_t* ms_p,
+ t_nrnb* nrnb_p,
+ gmx_wallcycle* wcycle_p,
+ bool pbcHandlingRequired,
+ ObservablesReducerBuilder* observablesReducerBuilder,
+ int numConstraints,
+ int numSettles) :
ncon_tot(numConstraints),
mtop(mtop_p),
pbcHandlingRequired_(pbcHandlingRequired),
if (ir.eConstrAlg == ConstraintAlgorithm::Lincs)
{
lincsd = init_lincs(
- log, mtop, nflexcon, at2con_mt, mayHaveSplitConstraints, ir.nLincsIter, ir.nProjOrder);
+ log, mtop, nflexcon, at2con_mt, mayHaveSplitConstraints, ir.nLincsIter, ir.nProjOrder, observablesReducerBuilder);
}
if (ir.eConstrAlg == ConstraintAlgorithm::Shake)
class ArrayRefWithPadding;
template<typename>
class ListOfLists;
+class ObservablesReducerBuilder;
//! Describes supported flavours of constrained updates.
enum class ConstraintVariable : int
*
* Private to enforce use of makeConstraints() factory
* function. */
- Constraints(const gmx_mtop_t& mtop,
- const t_inputrec& ir,
- pull_t* pull_work,
- FILE* log,
- const t_commrec* cr,
- bool useUpdateGroups,
- const gmx_multisim_t* ms,
- t_nrnb* nrnb,
- gmx_wallcycle* wcycle,
- bool pbcHandlingRequired,
- int numConstraints,
- int numSettles);
+ Constraints(const gmx_mtop_t& mtop,
+ const t_inputrec& ir,
+ pull_t* pull_work,
+ FILE* log,
+ const t_commrec* cr,
+ bool useUpdateGroups,
+ const gmx_multisim_t* ms,
+ t_nrnb* nrnb,
+ gmx_wallcycle* wcycle,
+ bool pbcHandlingRequired,
+ ObservablesReducerBuilder* observablesReducerBuilder,
+ int numConstraints,
+ int numSettles);
public:
/*! \brief This member type helps implement a factory
//! Getter for use by domain decomposition.
ArrayRef<const std::vector<int>> atom2settle_moltype() const;
- /*! \brief Return the data for reduction for determining
- * constraint RMS relative deviations, or an empty ArrayRef
- * when not supported for any active constraints. */
- ArrayRef<real> rmsdData() const;
/*! \brief Return the RMSD of the constraints when available. */
real rmsd() const;
#include <cstdlib>
#include <algorithm>
+#include <optional>
#include <vector>
#include "gromacs/domdec/domdec.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/observablesreducer.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pbcutil/pbc_simd.h"
#include "gromacs/simd/simd.h"
/*! @} */
//! The Lagrange multipliers times -1.
std::vector<real, AlignedAllocator<real>> mlambda;
- //! Storage for the constraint RMS relative deviation output.
- std::array<real, 2> rmsdData = { { 0 } };
+ /*! \brief Callback used after constraining to require reduction
+ * of values later used to compute the constraint RMS relative
+ * deviation, so the latter can be output. */
+ std::optional<ObservablesReducerBuilder::CallbackToRequireReduction> callbackToRequireReduction;
+ /*! \brief View used for reducing the components of the global
+ * relative RMS constraint deviation.
+ *
+ * Can be written any time, but that is only useful when followed
+ * by a call of the callbackToRequireReduction. Useful to read
+ * only from the callback that the ObservablesReducer will later
+ * make after reduction. */
+ ArrayRef<double> rmsdReductionBuffer;
+ /*! \brief The value of the constraint RMS deviation after it has
+ * been computed.
+ *
+ * When DD is active, filled by the ObservablesReducer, otherwise
+ * filled directly here. */
+ std::optional<double> constraintRmsDeviation;
};
/*! \brief Define simd_width for memory allocation used for SIMD code */
static const int simd_width = 1;
#endif
-ArrayRef<real> lincs_rmsdData(Lincs* lincsd)
-{
- return lincsd->rmsdData;
-}
-
real lincs_rmsd(const Lincs* lincsd)
{
- if (lincsd->rmsdData[0] > 0)
+ if (lincsd->constraintRmsDeviation.has_value())
{
- return std::sqrt(lincsd->rmsdData[1] / lincsd->rmsdData[0]);
+ return real(lincsd->constraintRmsDeviation.value());
}
else
{
ArrayRef<const ListOfLists<int>> atomToConstraintsPerMolType,
bool bPLINCS,
int nIter,
- int nProjOrder)
+ int nProjOrder,
+ ObservablesReducerBuilder* observablesReducerBuilder)
{
// TODO this should become a unique_ptr
Lincs* li;
}
}
+ if (observablesReducerBuilder)
+ {
+ ObservablesReducerBuilder::CallbackFromBuilder callbackFromBuilder =
+ [li](ObservablesReducerBuilder::CallbackToRequireReduction c, gmx::ArrayRef<double> v) {
+ li->callbackToRequireReduction = std::move(c);
+ li->rmsdReductionBuffer = v;
+ };
+
+ // Make the callback that runs afer reduction.
+ ObservablesReducerBuilder::CallbackAfterReduction callbackAfterReduction = [li](gmx::Step /*step*/) {
+ if (li->rmsdReductionBuffer[0] > 0)
+ {
+ li->constraintRmsDeviation =
+ std::sqrt(li->rmsdReductionBuffer[1] / li->rmsdReductionBuffer[0]);
+ }
+ };
+
+ const int reductionValuesRequired = 2;
+ observablesReducerBuilder->addSubscriber(reductionValuesRequired,
+ std::move(callbackFromBuilder),
+ std::move(callbackAfterReduction));
+ }
+
return li;
}
}
set_lincs_matrix(li, invmass, lambda);
-
- li->rmsdData[0] = 0.0;
- li->rmsdData[1] = 0.0;
}
//! Issues a warning when LINCS constraints cannot be satisfied.
if (lincsd->nc == 0 && cr->dd == nullptr)
{
- if (computeRmsd)
- {
- lincsd->rmsdData = { { 0 } };
- }
-
return bOK;
}
if (computeRmsd)
{
- // This is reduced across domains in compute_globals and
- // reported to the log file.
- lincsd->rmsdData[0] = deviations.numConstraints;
- lincsd->rmsdData[1] = deviations.sumSquaredDeviation;
- }
- else
- {
- // This is never read
- lincsd->rmsdData = { { 0 } };
+ if (lincsd->callbackToRequireReduction.has_value())
+ {
+ // This is reduced across domains in compute_globals and
+ // reported to the log file.
+ lincsd->rmsdReductionBuffer[0] = deviations.numConstraints;
+ lincsd->rmsdReductionBuffer[1] = deviations.sumSquaredDeviation;
+
+ // Call the ObservablesReducer via the callback it
+ // gave us for the purpose.
+ ObservablesReducerStatus status =
+ lincsd->callbackToRequireReduction.value()(ReductionRequirement::Soon);
+ GMX_RELEASE_ASSERT(status == ObservablesReducerStatus::ReadyToReduce,
+ "The LINCS RMSD is computed after observables have been "
+ "reduced, please reorder them.");
+ }
+ else
+ {
+ // Compute the deviation directly
+ lincsd->constraintRmsDeviation =
+ std::sqrt(deviations.sumSquaredDeviation / deviations.numConstraints);
+ }
}
if (printDebugOutput)
{
class Lincs;
template<typename>
class ListOfLists;
-
+class ObservablesReducerBuilder;
/*! \brief Return the data for determining constraint RMS relative deviations. */
ArrayRef<real> lincs_rmsdData(Lincs* lincsd);
ArrayRef<const ListOfLists<int>> atomsToConstraintsPerMolType,
bool bPLINCS,
int nIter,
- int nProjOrder);
+ int nProjOrder,
+ ObservablesReducerBuilder* observablesReducerBuilder);
/*! \brief Destructs the lincs object when it is not nullptr. */
void done_lincs(Lincs* li);
tensor shake_vir,
tensor total_vir,
tensor pres,
- gmx::ArrayRef<real> constraintsRmsdData,
gmx::SimulationSignaller* signalCoordinator,
const matrix lastbox,
gmx_bool* bSumEkinhOld,
shake_vir,
*ir,
ekind,
- constraintsRmsdData,
bStopCM ? vcm : nullptr,
signalBuffer,
*bSumEkinhOld,
{
template<typename T>
class ArrayRef;
-class Constraints;
class MDLogger;
class ObservablesReducer;
class SimulationSignaller;
tensor shake_vir,
tensor total_vir,
tensor pres,
- gmx::ArrayRef<real> constraintsRmsdData,
gmx::SimulationSignaller* signalCoordinator,
const matrix lastbox,
gmx_bool* bSumEkinhOld,
tensor svir,
const t_inputrec& inputrec,
gmx_ekindata_t* ekind,
- gmx::ArrayRef<real> constraintsRmsdData,
t_vcm* vcm,
gmx::ArrayRef<real> sig,
bool bSumEkinhOld,
gmx::ObservablesReducer* observablesReducer)
/* instead of current system, gmx_booleans for summing virial, kinetic energy, and other terms */
{
- int ie = 0, ifv = 0, isv = 0, irmsd = 0;
+ int ie = 0, ifv = 0, isv = 0;
int idedl = 0, idedlo = 0, idvdll = 0, idvdlnl = 0, iepl = 0, icm = 0, imass = 0, ica = 0;
int isig = -1;
int icj = -1, ici = -1, icx = -1;
if (bEner)
{
ie = add_binr(rb, nener, copyenerd.data());
- if (!constraintsRmsdData.empty())
- {
- irmsd = add_binr(rb, 2, constraintsRmsdData.data());
- }
for (auto key : gmx::keysOf(inn))
{
inn[key] = add_binr(rb, enerd->grpp.nener, enerd->grpp.energyGroupPairTerms[key].data());
if (bEner)
{
extract_binr(rb, ie, nener, copyenerd.data());
- if (!constraintsRmsdData.empty())
- {
- extract_binr(rb, irmsd, constraintsRmsdData);
- }
for (auto key : gmx::keysOf(inn))
{
extract_binr(rb, inn[key], enerd->grpp.nener, enerd->grpp.energyGroupPairTerms[key].data());
{
template<typename T>
class ArrayRef;
-class Constraints;
class ObservablesReducer;
} // namespace gmx
tensor svir,
const t_inputrec& inputrec,
gmx_ekindata_t* ekind,
- gmx::ArrayRef<real> constraintsRmsdData,
t_vcm* vcm,
gmx::ArrayRef<real> sig,
bool bSumEkinhOld,
at2con_mt,
false,
testData->ir_.nLincsIter,
- testData->ir_.nProjOrder);
+ testData->ir_.nProjOrder,
+ nullptr);
set_lincs(*testData->idef_,
testData->numAtoms_,
testData->invmass_,
// TODO This object will always return zero as RMSD value.
// It is more relevant to have non-zero value for testing.
constraints_ = makeConstraints(
- mtop_, inputrec_, nullptr, false, nullptr, &cr_, false, nullptr, nullptr, nullptr, false);
+ mtop_, inputrec_, nullptr, false, nullptr, &cr_, false, nullptr, nullptr, nullptr, false, nullptr);
}
/*! \brief Helper function to generate synthetic data to output
shake_vir,
total_vir,
pres,
- (bCalcEner && constr != nullptr) ? constr->rmsdData() : gmx::ArrayRef<real>{},
nullSignaller,
state->box,
bSumEkinhOld,
nullptr,
nullptr,
nullptr,
- gmx::ArrayRef<real>{},
nullSignaller,
state->box,
bSumEkinhOld,
shake_vir,
total_vir,
pres,
- gmx::ArrayRef<real>{},
nullSignaller,
lastbox,
bSumEkinhOld,
shake_vir,
total_vir,
pres,
- gmx::ArrayRef<real>{},
&nullSignaller,
state->box,
&bSumEkinhOld,
shake_vir,
total_vir,
pres,
- gmx::ArrayRef<real>{},
&nullSignaller,
state->box,
&bSumEkinhOld,
nullptr,
nullptr,
nullptr,
- gmx::ArrayRef<real>{},
&nullSignaller,
state->box,
&bSumEkinhOld,
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,
- md,
- 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,
- &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);
+ 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(
shake_vir,
total_vir,
pres,
- gmx::ArrayRef<real>{},
&nullSignaller,
state->box,
&bSumEkinhOld,
nullptr,
nullptr,
nullptr,
- constr != nullptr ? constr->rmsdData() : gmx::ArrayRef<real>{},
&signaller,
state->box,
&bSumEkinhOld,
shake_vir,
*inputrec,
nullptr,
- gmx::ArrayRef<real>{},
nullptr,
std::vector<real>(1, terminate),
FALSE,
shake_vir,
total_vir,
pres,
- gmx::ArrayRef<real>{},
&nullSignaller,
state->box,
&bSumEkinhOld,
shake_vir,
total_vir,
pres,
- constr != nullptr ? constr->rmsdData() : gmx::ArrayRef<real>{},
&signaller,
state->box,
&bSumEkinhOld,
ms,
&nrnb,
wcycle.get(),
- fr->bMolPBC);
+ fr->bMolPBC,
+ &observablesReducerBuilder);
/* Energy terms and groups */
gmx_enerdata_t enerd(mtop.groups.groups[SimulationAtomGroupType::EnergyOutput].size(),
energyData_->constraintVirial(step),
energyData_->totalVirial(step),
energyData_->pressure(step),
- (((flags & CGLO_ENERGY) != 0) && constr_ != nullptr) ? constr_->rmsdData()
- : gmx::ArrayRef<real>{},
signaller,
lastbox,
energyData_->needToSumEkinhOld(),
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff