return *pbcType_;
}
+ //! Set the simulation time step
+ void setSimulationTimeStep(double timeStep)
+ {
+ simulationTimeStep_ = timeStep;
+ }
+
+ //! Return the simulation time step
+ double simulationTimeStep()
+ {
+ return simulationTimeStep_;
+ }
+
private:
//! The reference density to fit to
std::unique_ptr<MultiDimArray<std::vector<float>, dynamicExtents3D> > referenceDensity_;
std::unique_ptr<LocalAtomSet> localAtomSet_;
//! The type of periodic boundary conditions in the simulation
std::unique_ptr<int> pbcType_;
+ //! The simulation time step
+ double simulationTimeStep_ = 1;
GMX_DISALLOW_COPY_AND_ASSIGN(DensityFittingSimulationParameterSetup);
};
};
notifier->notifier_.subscribe(setPeriodicBoundaryContionsFunction);
+ // setting the simulation time step
+ const auto setSimulationTimeStepFunction = [this](const SimulationTimeStep &simulationTimeStep) {
+ this->densityFittingSimulationParameters_.setSimulationTimeStep(simulationTimeStep.delta_t);
+ };
+ notifier->notifier_.subscribe(setSimulationTimeStepFunction);
+
// adding output to energy file
const auto requestEnergyOutput
= [this](MdModulesEnergyOutputToDensityFittingRequestChecker *energyOutputRequest) {
densityFittingSimulationParameters_.transformationToDensityLattice(),
densityFittingSimulationParameters_.localAtomSet(),
densityFittingSimulationParameters_.periodicBoundaryConditionType(),
+ densityFittingSimulationParameters_.simulationTimeStep(),
densityFittingState_);
forceProviders->addForceProvider(forceProvider_.get());
}
checkpointWriting.builder_.addValue<std::int64_t>(
DensityFittingModuleInfo::name_ + "-stepsSinceLastCalculation",
state.stepsSinceLastCalculation_);
+ checkpointWriting.builder_.addValue<real>(
+ DensityFittingModuleInfo::name_ + "-adaptiveForceConstantScale",
+ state.adaptiveForceConstantScale_);
+ KeyValueTreeObjectBuilder exponentialMovingAverageKvtEntry =
+ checkpointWriting.builder_.addObject(
+ DensityFittingModuleInfo::name_ + "-exponentialMovingAverageState");
+ exponentialMovingAverageStateAsKeyValueTree(exponentialMovingAverageKvtEntry, state.exponentialMovingAverageState_);
}
}
checkpointReading.checkpointedData_[
DensityFittingModuleInfo::name_ + "-stepsSinceLastCalculation"].cast<std::int64_t>();
}
+ if (checkpointReading.checkpointedData_.keyExists(DensityFittingModuleInfo::name_ + "-adaptiveForceConstantScale"))
+ {
+ densityFittingState_.adaptiveForceConstantScale_
+ = checkpointReading
+ .checkpointedData_[DensityFittingModuleInfo::name_ + "-adaptiveForceConstantScale"].cast<real>();
+ }
+ if (checkpointReading.checkpointedData_.keyExists(
+ DensityFittingModuleInfo::name_ + "-exponentialMovingAverageState"))
+ {
+ densityFittingState_.exponentialMovingAverageState_ =
+ exponentialMovingAverageStateFromKeyValueTree(checkpointReading
+ .checkpointedData_[DensityFittingModuleInfo::name_ + "-exponentialMovingAverageState"].asObject());
+ }
}
}
if (PAR(&(checkpointBroadcast.cr_)))
{
block_bc(&(checkpointBroadcast.cr_), densityFittingState_.stepsSinceLastCalculation_);
+ block_bc(&(checkpointBroadcast.cr_), densityFittingState_.adaptiveForceConstantScale_);
+ block_bc(&(checkpointBroadcast.cr_), densityFittingState_.exponentialMovingAverageState_);
}
}
}
#include <numeric>
+#include "gromacs/compat/optional.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/math/densityfit.h"
#include "gromacs/math/densityfittingforce.h"
+#include "gromacs/math/exponentialmovingaverage.h"
#include "gromacs/math/gausstransform.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/enerdata.h"
const TranslateAndScale &transformationToDensityLattice,
const LocalAtomSet &localAtomSet,
int pbcType,
+ double simulationTimeStep,
const DensityFittingForceProviderState &state);
~Impl();
void calculateForces(const ForceProviderInput &forceProviderInput, ForceProviderOutput *forceProviderOutput);
TranslateAndScale transformationToDensityLattice_;
RVec referenceDensityCenter_;
int pbcType_;
- real forceConstantScale_;
+ //! Optionally scale the force according to a moving average of the similarity
+ compat::optional<ExponentialMovingAverage> expAverageSimilarity_;
};
DensityFittingForceProvider::Impl::~Impl() = default;
const TranslateAndScale &transformationToDensityLattice,
const LocalAtomSet &localAtomSet,
int pbcType,
+ double simulationTimeStep,
const DensityFittingForceProviderState &state) :
parameters_(parameters),
state_(state),
amplitudeLookup_(parameters_.amplitudeLookupMethod_),
transformationToDensityLattice_(transformationToDensityLattice),
pbcType_(pbcType),
- forceConstantScale_(parameters_.calculationIntervalInSteps_)
+ expAverageSimilarity_(compat::nullopt)
{
+ if (parameters_.adaptiveForceScaling_)
+ {
+ GMX_ASSERT(simulationTimeStep > 0, "Simulation time step must be larger than zero for adaptive for scaling.");
+ expAverageSimilarity_.emplace(ExponentialMovingAverage(
+ parameters_.adaptiveForceScalingTimeConstant_
+ / (simulationTimeStep * parameters_.calculationIntervalInSteps_),
+ state.exponentialMovingAverageState_));
+ }
referenceDensityCenter_ = {
real(referenceDensity.extent(XX))/2,
real(referenceDensity.extent(YY))/2,
transformationToDensityLattice_.scaleOperationOnly().inverseIgnoringZeroScale(forces_);
- auto densityForceIterator = forces_.cbegin();
+ auto densityForceIterator = forces_.cbegin();
+ const real effectiveForceConstant = state_.adaptiveForceConstantScale_ *
+ parameters_.calculationIntervalInSteps_ * parameters_.forceConstant_;
for (const auto localAtomIndex : localAtomSet_.localIndex())
{
- forceProviderOutput->forceWithVirial_.force_[localAtomIndex] +=
- forceConstantScale_ * parameters_.forceConstant_ * *densityForceIterator;
+ forceProviderOutput->forceWithVirial_.force_[localAtomIndex]
+ += effectiveForceConstant * *densityForceIterator;
++densityForceIterator;
}
// calculate corresponding potential energy
const float similarity = measure_.similarity(gaussTransform_.constView());
- const real energy = -similarity * parameters_.forceConstant_;
+ const real energy = -similarity * parameters_.forceConstant_ * state_.adaptiveForceConstantScale_;
forceProviderOutput->enerd_.term[F_DENSITYFITTING] += energy;
+
+ if (expAverageSimilarity_.has_value())
+ {
+ expAverageSimilarity_->updateWithDataPoint(similarity);
+
+ if (expAverageSimilarity_->increasing())
+ {
+ state_.adaptiveForceConstantScale_ /= 1._real + expAverageSimilarity_->inverseTimeConstant();
+ }
+ else
+ {
+ state_.adaptiveForceConstantScale_ *= 1._real + expAverageSimilarity_->inverseTimeConstant();
+ }
+ }
}
DensityFittingForceProviderState
DensityFittingForceProvider::Impl::state()
{
+ if (expAverageSimilarity_.has_value())
+ {
+ state_.exponentialMovingAverageState_ = expAverageSimilarity_->state();
+ }
return state_;
}
const TranslateAndScale &transformationToDensityLattice,
const LocalAtomSet &localAtomSet,
int pbcType,
+ double simulationTimeStep,
const DensityFittingForceProviderState &state)
- : impl_(new Impl(parameters, referenceDensity, transformationToDensityLattice, localAtomSet, pbcType, state))
+ : impl_(new Impl(parameters, referenceDensity, transformationToDensityLattice, localAtomSet, pbcType, simulationTimeStep, state))
{}
void DensityFittingForceProvider::calculateForces(const ForceProviderInput &forceProviderInput,
#include "gromacs/domdec/localatomset.h"
#include "gromacs/math/coordinatetransformation.h"
+#include "gromacs/math/exponentialmovingaverage.h"
#include "gromacs/mdspan/extensions.h"
#include "gromacs/mdtypes/iforceprovider.h"
#include "gromacs/utility/classhelpers.h"
/*! \brief The steps since the last force calculation.
* Used if density fitting is to be calculated every N steps.
*/
- std::int64_t stepsSinceLastCalculation_ = 0;
+ std::int64_t stepsSinceLastCalculation_ = 0;
+ //! The state of the exponential moving average of the similarity measure
+ ExponentialMovingAverageState exponentialMovingAverageState_ = {};
+ //! An additional factor scaling the force for adaptive force scaling
+ real adaptiveForceConstantScale_ = 1.0_real;
};
/*! \internal \brief
const TranslateAndScale &transformationToDensityLattice,
const LocalAtomSet &localAtomSet,
int pbcType,
+ double simulationTimeStep,
const DensityFittingForceProviderState &state);
~DensityFittingForceProvider();
/*!\brief Calculate forces that maximise goodness-of-fit with a reference density map.
densityfittingMdpTransformFromString<std::string>(rules, stringIdentityTransform, c_referenceDensityFileNameTag_);
densityfittingMdpTransformFromString<std::int64_t>(rules, &fromStdString<std::int64_t>, c_everyNStepsTag_);
densityfittingMdpTransformFromString<bool>(rules, &fromStdString<bool>, c_normalizeDensitiesTag_);
+ densityfittingMdpTransformFromString<bool>(rules, &fromStdString<bool>, c_adaptiveForceScalingTag_);
+ densityfittingMdpTransformFromString<real>(rules, &fromStdString<real>, c_adaptiveForceScalingTimeConstantTag_);
}
void DensityFittingOptions::buildMdpOutput(KeyValueTreeObjectBuilder *builder) const
addDensityFittingMdpOutputValue(builder, parameters_.calculationIntervalInSteps_, c_everyNStepsTag_);
addDensityFittingMdpOutputValueComment(builder, "; Normalize the sum of density voxel values to one", c_normalizeDensitiesTag_);
addDensityFittingMdpOutputValue(builder, parameters_.normalizeDensities_, c_normalizeDensitiesTag_);
+ addDensityFittingMdpOutputValueComment(builder, "; Apply adaptive force scaling", c_adaptiveForceScalingTag_);
+ addDensityFittingMdpOutputValue(builder, parameters_.adaptiveForceScaling_, c_adaptiveForceScalingTag_);
+ addDensityFittingMdpOutputValueComment(builder, "; Time constant for adaptive force scaling in ps", c_adaptiveForceScalingTimeConstantTag_);
+ addDensityFittingMdpOutputValue(builder, parameters_.adaptiveForceScalingTimeConstant_, c_adaptiveForceScalingTimeConstantTag_);
}
}
section.addOption(StringOption(c_referenceDensityFileNameTag_.c_str()).store(&referenceDensityFileName_));
section.addOption(Int64Option(c_everyNStepsTag_.c_str()).store(¶meters_.calculationIntervalInSteps_));
section.addOption(BooleanOption(c_normalizeDensitiesTag_.c_str()).store(¶meters_.normalizeDensities_));
+ section.addOption(BooleanOption(c_adaptiveForceScalingTag_.c_str()).store(¶meters_.adaptiveForceScaling_));
+ section.addOption(RealOption(c_adaptiveForceScalingTimeConstantTag_.c_str()).store(¶meters_.adaptiveForceScalingTimeConstant_));
}
bool DensityFittingOptions::active() const
const std::string c_normalizeDensitiesTag_ = "normalize-densities";
+ const std::string c_adaptiveForceScalingTag_ = "adaptive-force-scaling";
+
+ const std::string c_adaptiveForceScalingTimeConstantTag_ = "adaptive-force-scaling-time-constant";
+
DensityFittingParameters parameters_;
};
std::int64_t calculationIntervalInSteps_ = 1;
//! Normalize reference and simulated densities
bool normalizeDensities_ = true;
+ //! Perform adaptive force scaling during the simulation
+ bool adaptiveForceScaling_ = false;
+ //! The time constant for the adaptive force scaling in ps
+ real adaptiveForceScalingTimeConstant_ = 4;
};
/*!\brief Check if two structs holding density fitting parameters are equal.
"density-guided-simulation-nst = 1\n"
"; Normalize the sum of density voxel values to one\n"
"density-guided-simulation-normalize-densities = true\n"
+ "; Apply adaptive force scaling\n"
+ "density-guided-simulation-adaptive-force-scaling = false\n"
+ "; Time constant for adaptive force scaling in ps\n"
+ "density-guided-simulation-adaptive-force-scaling-time-constant = 4\n"
};
EXPECT_EQ(expectedString, stream.toString());
#include "exponentialmovingaverage.h"
#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/keyvaluetree.h"
namespace gmx
{
-ExponentialMovingAverage::ExponentialMovingAverage(real timeConstant, const ExponentialMovingAverageState &state) :
- state_(state)
+//! Convert the exponential moving average state as key-value-tree object
+void exponentialMovingAverageStateAsKeyValueTree(KeyValueTreeObjectBuilder builder, const ExponentialMovingAverageState &state)
+{
+ builder.addValue<real>("weighted-sum", state.weightedSum_);
+ builder.addValue<real>("weighted-count", state.weightedCount_);
+ builder.addValue<bool>("increasing", state.increasing_);
+}
+
+//! Sets the exponential moving average state from a key-value-tree object
+ExponentialMovingAverageState
+exponentialMovingAverageStateFromKeyValueTree(const KeyValueTreeObject &object)
+{
+ const real weightedSum = object["weighted-sum"].cast<real>();
+ const real weightedCount = object["weighted-count"].cast<real>();
+ const bool increasing = object["increasing"].cast<bool>();
+ return {weightedSum, weightedCount, increasing};
+}
+
+ExponentialMovingAverage::ExponentialMovingAverage(real timeConstant, const ExponentialMovingAverageState &state)
+ : state_(state)
{
if (timeConstant < 1)
{
#ifndef GMX_MATH_EXPONENTIALMOVINGAVERAGE_H
#define GMX_MATH_EXPONENTIALMOVINGAVERAGE_H
+#include "gromacs/utility/keyvaluetreebuilder.h"
#include "gromacs/utility/real.h"
namespace gmx
bool increasing_ = false;
};
+//! Convert the exponential moving average state as key-value-tree object
+void exponentialMovingAverageStateAsKeyValueTree(KeyValueTreeObjectBuilder builder, const ExponentialMovingAverageState &state);
+
+//! Sets the expoential moving average state from a key-value-tree object
+ExponentialMovingAverageState
+exponentialMovingAverageStateFromKeyValueTree(const KeyValueTreeObject &object);
+
/*! \libinternal
* \brief Evaluate the exponential moving average with bias correction.
*
EXPECT_REAL_EQ(0.5, exponentialMovingAverage.inverseTimeConstant());
}
+TEST(ExponentialMovingAverage, RoundTripAsKeyValueTree)
+{
+ KeyValueTreeBuilder builder;
+ const real weightedSum = 9;
+ const real weightedCount = 1;
+ const bool increasing = true;
+ ExponentialMovingAverageState state = {weightedSum, weightedCount, increasing};
+ exponentialMovingAverageStateAsKeyValueTree(builder.rootObject(), state);
+ state = {};
+ KeyValueTreeObject result = builder.build();
+ state = exponentialMovingAverageStateFromKeyValueTree(result);
+ EXPECT_EQ(weightedSum, state.weightedSum_);
+ EXPECT_EQ(weightedCount, state.weightedCount_);
+ EXPECT_EQ(increasing, state.increasing_);
+}
+
} // namespace
} // namespace test
mdModulesNotifier.notify(*cr);
mdModulesNotifier.notify(&atomSets);
mdModulesNotifier.notify(PeriodicBoundaryConditionType {inputrec->ePBC});
+ mdModulesNotifier.notify(SimulationTimeStep { inputrec->delta_t });
/* Initiate forcerecord */
fr = new t_forcerec;
fr->forceProviders = mdModules_->initForceProviders();
std::vector<std::string> errorMessages_;
};
+struct SimulationTimeStep
+{
+ //! Time step (ps)
+ const double delta_t;
+};
+
struct MdModulesNotifier
{
//! Register callback function types for MdModule
MdModulesCheckpointReadingDataOnMaster,
MdModulesCheckpointReadingBroadcast,
MdModulesWriteCheckpointData,
- PeriodicBoundaryConditionType>::type notifier_;
+ PeriodicBoundaryConditionType,
+ const SimulationTimeStep &>::type notifier_;
};
} // namespace gmx
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff