/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2015,2016,2017,2018,2019, The GROMACS development team.
+ * Copyright (c) 2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/xvgr.h"
#include "gromacs/gmxlib/network.h"
+#include "gromacs/math/units.h"
#include "gromacs/math/utilities.h"
#include "gromacs/mdrunutility/multisim.h"
#include "gromacs/mdtypes/awh_history.h"
#include "gromacs/utility/stringutil.h"
#include "biasgrid.h"
+#include "biassharing.h"
#include "pointstate.h"
namespace gmx
namespace
{
-/*! \brief
- * Sum an array over all simulations on the master rank of each simulation.
- *
- * \param[in,out] arrayRef The data to sum.
- * \param[in] multiSimComm Struct for multi-simulation communication.
- */
-void sumOverSimulations(gmx::ArrayRef<int> arrayRef, const gmx_multisim_t* multiSimComm)
-{
- gmx_sumi_sim(arrayRef.size(), arrayRef.data(), multiSimComm);
-}
-
-/*! \brief
- * Sum an array over all simulations on the master rank of each simulation.
- *
- * \param[in,out] arrayRef The data to sum.
- * \param[in] multiSimComm Struct for multi-simulation communication.
- */
-void sumOverSimulations(gmx::ArrayRef<double> arrayRef, const gmx_multisim_t* multiSimComm)
-{
- gmx_sumd_sim(arrayRef.size(), arrayRef.data(), multiSimComm);
-}
-
-/*! \brief
- * Sum an array over all simulations on all ranks of each simulation.
- *
- * This assumes the data is identical on all ranks within each simulation.
- *
- * \param[in,out] arrayRef The data to sum.
- * \param[in] commRecord Struct for intra-simulation communication.
- * \param[in] multiSimComm Struct for multi-simulation communication.
- */
-template<typename T>
-void sumOverSimulations(gmx::ArrayRef<T> arrayRef, const t_commrec* commRecord, const gmx_multisim_t* multiSimComm)
-{
- if (MASTER(commRecord))
- {
- sumOverSimulations(arrayRef, multiSimComm);
- }
- if (commRecord->nnodes > 1)
- {
- gmx_bcast(arrayRef.size() * sizeof(T), arrayRef.data(), commRecord->mpi_comm_mygroup);
- }
-}
-
/*! \brief
* Sum PMF over multiple simulations, when requested.
*
* \param[in,out] pointState The state of the points in the bias.
* \param[in] numSharedUpdate The number of biases sharing the histogram.
- * \param[in] commRecord Struct for intra-simulation communication.
- * \param[in] multiSimComm Struct for multi-simulation communication.
+ * \param[in] biasSharing Object for sharing bias data over multiple simulations
+ * \param[in] biasIndex Index of this bias in the total list of biases in this simulation
*/
-void sumPmf(gmx::ArrayRef<PointState> pointState,
- int numSharedUpdate,
- const t_commrec* commRecord,
- const gmx_multisim_t* multiSimComm)
+void sumPmf(gmx::ArrayRef<PointState> pointState, int numSharedUpdate, const BiasSharing* biasSharing, const int biasIndex)
{
if (numSharedUpdate == 1)
{
return;
}
- GMX_ASSERT(multiSimComm != nullptr && numSharedUpdate % multiSimComm->numSimulations_ == 0,
+ GMX_ASSERT(biasSharing != nullptr
+ && numSharedUpdate % biasSharing->numSharingSimulations(biasIndex) == 0,
"numSharedUpdate should be a multiple of multiSimComm->numSimulations_");
- GMX_ASSERT(numSharedUpdate == multiSimComm->numSimulations_,
+ GMX_ASSERT(numSharedUpdate == biasSharing->numSharingSimulations(biasIndex),
"Sharing within a simulation is not implemented (yet)");
std::vector<double> buffer(pointState.size());
/* Need to temporarily exponentiate the log weights to sum over simulations */
for (size_t i = 0; i < buffer.size(); i++)
{
- buffer[i] = pointState[i].inTargetRegion() ? std::exp(-pointState[i].logPmfSum()) : 0;
+ buffer[i] = pointState[i].inTargetRegion() ? std::exp(pointState[i].logPmfSum()) : 0;
}
- sumOverSimulations(gmx::ArrayRef<double>(buffer), commRecord, multiSimComm);
+ biasSharing->sumOverSharingSimulations(gmx::ArrayRef<double>(buffer), biasIndex);
/* Take log again to get (non-normalized) PMF */
double normFac = 1.0 / numSharedUpdate;
{
if (pointState[i].inTargetRegion())
{
- pointState[i].setLogPmfSum(-std::log(buffer[i] * normFac));
+ pointState[i].setLogPmfSum(std::log(buffer[i] * normFac));
}
}
}
* \param[in] gridpointIndex The index of the current grid point.
* \returns the log of the biased probability weight.
*/
-double biasedLogWeightFromPoint(const std::vector<DimParams>& dimParams,
- const std::vector<PointState>& points,
- const BiasGrid& grid,
- int pointIndex,
- double pointBias,
- const awh_dvec value,
- gmx::ArrayRef<const double> neighborLambdaEnergies,
- int gridpointIndex)
+double biasedLogWeightFromPoint(ArrayRef<const DimParams> dimParams,
+ ArrayRef<const PointState> points,
+ const BiasGrid& grid,
+ int pointIndex,
+ double pointBias,
+ const awh_dvec value,
+ ArrayRef<const double> neighborLambdaEnergies,
+ int gridpointIndex)
{
double logWeight = detail::c_largeNegativeExponent;
* \returns The calculated marginal distribution in a 1D array with
* as many elements as there are points along the axis of interest.
*/
-std::vector<double> calculateFELambdaMarginalDistribution(const BiasGrid& grid,
- gmx::ArrayRef<const int> neighbors,
- gmx::ArrayRef<const double> probWeightNeighbor)
+std::vector<double> calculateFELambdaMarginalDistribution(const BiasGrid& grid,
+ ArrayRef<const int> neighbors,
+ ArrayRef<const double> probWeightNeighbor)
{
const std::optional<int> lambdaAxisIndex = grid.lambdaAxisIndex();
GMX_RELEASE_ASSERT(lambdaAxisIndex,
} // namespace
-void BiasState::calcConvolvedPmf(const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- std::vector<float>* convolvedPmf) const
+void BiasState::calcConvolvedPmf(ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ std::vector<float>* convolvedPmf) const
{
size_t numPoints = grid.numPoints();
{
double freeEnergyWeights = 0;
const GridPoint& point = grid.point(m);
- for (auto& neighbor : point.neighbor)
+ for (const auto& neighbor : point.neighbor)
{
/* Do not convolve the bias along a lambda axis - only use the pmf from the current point */
if (!pointsHaveDifferentLambda(grid, m, neighbor))
* \param[in,out] pointState The state of all points.
* \param[in] params The bias parameters.
*/
-void updateTargetDistribution(
gmx::ArrayRef<PointState> pointState, const BiasParams& params)
+void updateTargetDistribution(ArrayRef<PointState> pointState, const BiasParams& params)
{
double freeEnergyCutoff = 0;
- if (params.eTarget == eawhtargetCUTOFF)
+ if (params.eTarget == AwhTargetType::Cutoff)
{
freeEnergyCutoff = freeEnergyMinimumValue(pointState) + params.freeEnergyCutoffInKT;
}
/* Sum up the histograms and get their normalization */
double sumVisits = 0;
double sumWeights = 0;
- for (auto& pointState : points_)
+ for (const auto& pointState : points_)
{
if (pointState.inTargetRegion())
{
return numWarnings;
}
-double BiasState::calcUmbrellaForceAndPotential(const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- int point,
- ArrayRef<const double> neighborLambdaDhdl,
- gmx::ArrayRef<double> force) const
+double BiasState::calcUmbrellaForceAndPotential(ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ int point,
+ ArrayRef<const double> neighborLambdaDhdl,
+ ArrayRef<double> force) const
{
double potential = 0;
for (size_t d = 0; d < dimParams.size(); d++)
return potential;
}
-void BiasState::calcConvolvedForce(const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- gmx::ArrayRef<const double> probWeightNeighbor,
- ArrayRef<const double> neighborLambdaDhdl,
- gmx::ArrayRef<double> forceWorkBuffer,
- gmx::ArrayRef<double> force) const
+void BiasState::calcConvolvedForce(ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ ArrayRef<const double> probWeightNeighbor,
+ ArrayRef<const double> neighborLambdaDhdl,
+ ArrayRef<double> forceWorkBuffer,
+ ArrayRef<double> force) const
{
for (size_t d = 0; d < dimParams.size(); d++)
{
}
}
-double BiasState::moveUmbrella(const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- gmx::ArrayRef<const double> probWeightNeighbor,
- ArrayRef<const double> neighborLambdaDhdl,
- gmx::ArrayRef<double> biasForce,
- int64_t step,
- int64_t seed,
- int indexSeed,
- bool onlySampleUmbrellaGridpoint)
+double BiasState::moveUmbrella(ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ ArrayRef<const double> probWeightNeighbor,
+ ArrayRef<const double> neighborLambdaDhdl,
+ ArrayRef<double> biasForce,
+ int64_t step,
+ int64_t seed,
+ int indexSeed,
+ bool onlySampleUmbrellaGridpoint)
{
/* Generate and set a new coordinate reference value */
coordState_.sampleUmbrellaGridpoint(
/* For each neighbor point of the center point, refresh its state by adding the results of all past, skipped updates. */
const std::vector<int>& neighbors = grid.point(coordState_.gridpointIndex()).neighbor;
- for (auto& neighbor : neighbors)
+ for (const auto& neighbor : neighbors)
{
bool didUpdate = points_[neighbor].performPreviouslySkippedUpdates(
params, histogramSize_.numUpdates(), weightHistScaling, logPmfsumScaling);
*
* \param[in,out] updateList Update list for this simulation (assumed >= npoints long).
* \param[in] numPoints Total number of points.
- * \param[in] commRecord Struct for intra-simulation communication.
- * \param[in] multiSimComm Struct for multi-simulation communication.
+ * \param[in] biasSharing Object for sharing bias data over multiple simulations
+ * \param[in] biasIndex Index of this bias in the total list of biases in this simulation
*/
-void mergeSharedUpdateLists(std::vector<int>* updateList,
- int numPoints,
- const t_commrec* commRecord,
- const gmx_multisim_t* multiSimComm)
+void mergeSharedUpdateLists(std::vector<int>* updateList,
+ int numPoints,
+ const BiasSharing& biasSharing,
+ const int biasIndex)
{
std::vector<int> numUpdatesOfPoint;
}
/* Sum over the sims to get all the flagged points */
- sumOverSimulations(arrayRefFromArray(numUpdatesOfPoint.data(), numPoints), commRecord, multiSimComm);
+ biasSharing.sumOverSharingSimulations(arrayRefFromArray(numUpdatesOfPoint.data(), numPoints), biasIndex);
/* Collect the indices of the flagged points in place. The resulting array will be the merged update list.*/
updateList->clear();
* \param[in] grid The AWH bias.
* \param[in] points The point state.
* \param[in] originUpdatelist The origin of the rectangular region that has been sampled since
- * last update. \param[in] endUpdatelist The end of the rectangular that has been sampled since
- * last update. \param[in,out] updateList Local update list to set (assumed >= npoints long).
+ * last update.
+ * \param[in] endUpdatelist The end of the rectangular that has been sampled since
+ * last update.
+ * \param[in,out] updateList Local update list to set (assumed >= npoints long).
*/
-void makeLocalUpdateList(const BiasGrid& grid,
- const std::vector<PointState>& points,
- const awh_ivec originUpdatelist,
- const awh_ivec endUpdatelist,
- std::vector<int>* updateList)
+void makeLocalUpdateList(const BiasGrid& grid,
+ ArrayRef<const PointState> points,
+ const awh_ivec originUpdatelist,
+ const awh_ivec endUpdatelist,
+ std::vector<int>* updateList)
{
awh_ivec origin;
awh_ivec numPoints;
* \param[in,out] pointState The state of the points in the bias.
* \param[in,out] weightSumCovering The weights for checking covering.
* \param[in] numSharedUpdate The number of biases sharing the histrogram.
- * \param[in] commRecord Struct for intra-simulation communication.
- * \param[in] multiSimComm Struct for multi-simulation communication.
- * \param[in] localUpdateList List of points with data.
+ * \param[in] biasSharing Object for sharing bias data over multiple simulations
+ * \param[in] biasIndex Index of this bias in the total list of biases in this
+ * simulation \param[in] localUpdateList List of points with data.
*/
void sumHistograms(gmx::ArrayRef<PointState> pointState,
gmx::ArrayRef<double> weightSumCovering,
int numSharedUpdate,
- const t_commrec* commRecord,
- const gmx_multisim_t* multiSimComm,
+ const BiasSharing* biasSharing,
+ const int biasIndex,
const std::vector<int>& localUpdateList)
{
/* The covering checking histograms are added before summing over simulations, so that the
/* Sum histograms over multiple simulations if needed. */
if (numSharedUpdate > 1)
{
- GMX_ASSERT(numSharedUpdate == multiSimComm->numSimulations_,
+ GMX_ASSERT(numSharedUpdate == biasSharing->numSharingSimulations(biasIndex),
"Sharing within a simulation is not implemented (yet)");
/* Collect the weights and counts in linear arrays to be able to use gmx_sumd_sim. */
coordVisits[localIndex] = ps.numVisitsIteration();
}
- sumOverSimulations(gmx::ArrayRef<double>(weightSum), commRecord, multiSimComm);
- sumOverSimulations(gmx::ArrayRef<double>(coordVisits), commRecord, multiSimComm);
+ biasSharing->sumOverSharingSimulations(gmx::ArrayRef<double>(weightSum), biasIndex);
+ biasSharing->sumOverSharingSimulations(gmx::ArrayRef<double>(coordVisits), biasIndex);
/* Transfer back the result */
for (size_t localIndex = 0; localIndex < localUpdateList.size(); localIndex++)
} // namespace
-bool BiasState::isSamplingRegionCovered(const BiasParams& params,
- const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- const t_commrec* commRecord,
- const gmx_multisim_t* multiSimComm) const
+bool BiasState::isSamplingRegionCovered(const BiasParams& params,
+ ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid) const
{
/* Allocate and initialize arrays: one for checking visits along each dimension,
one for keeping track of which points to check and one for the covered points.
/* Set the free energy cutoff */
double maxFreeEnergy = GMX_FLOAT_MAX;
- if (params.eTarget == eawhtargetCUTOFF)
+ if (params.eTarget == AwhTargetType::Cutoff)
{
maxFreeEnergy = freeEnergyMinimumValue(points_) + params.freeEnergyCutoffInKT;
}
/* For multiple dimensions this may not be the best way to do it. */
for (int d = 0; d < grid.numDimensions(); d++)
{
- sumOverSimulations(
+ biasSharing_->sumOverSharingSimulations(
gmx::arrayRefFromArray(checkDim[d].covered.data(), grid.axis(d).numPoints()),
- commRecord,
- multiSimComm);
+ params.biasIndex);
}
}
}
}
-void BiasState::updateFreeEnergyAndAddSamplesToHistogram(const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- const BiasParams& params,
- const t_commrec* commRecord,
- const gmx_multisim_t* multiSimComm,
- double t,
- int64_t step,
- FILE* fplog,
- std::vector<int>* updateList)
+void BiasState::updateFreeEnergyAndAddSamplesToHistogram(ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ const BiasParams& params,
+ double t,
+ int64_t step,
+ FILE* fplog,
+ std::vector<int>* updateList)
{
/* Note hat updateList is only used in this scope and is always
* re-initialized. We do not use a local vector, because that would
makeLocalUpdateList(grid, points_, originUpdatelist_, endUpdatelist_, updateList);
if (params.numSharedUpdate > 1)
{
- mergeSharedUpdateLists(updateList, points_.size(), commRecord, multiSimComm);
+ mergeSharedUpdateLists(updateList, points_.size(), *biasSharing_, params.biasIndex);
}
/* Reset the range for the next update */
resetLocalUpdateRange(grid);
/* Add samples to histograms for all local points and sync simulations if needed */
- sumHistograms(points_, weightSumCovering_, params.numSharedUpdate, commRecord, multiSimComm, *updateList);
+ sumHistograms(points_, weightSumCovering_, params.numSharedUpdate, biasSharing_, params.biasIndex, *updateList);
- sumPmf(points_, params.numSharedUpdate, commRecord, multiSimComm);
+ sumPmf(points_, params.numSharedUpdate, biasSharing_, params.biasIndex);
/* Renormalize the free energy if values are too large. */
bool needToNormalizeFreeEnergy = false;
/* Update target distribution? */
bool needToUpdateTargetDistribution =
- (params.eTarget != eawhtargetCONSTANT && params.isUpdateTargetStep(step));
+ (params.eTarget != AwhTargetType::Constant && params.isUpdateTargetStep(step));
/* In the initial stage, the histogram grows dynamically as a function of the number of coverings. */
bool detectedCovering = false;
if (inInitialStage())
{
detectedCovering =
- (params.isCheckCoveringStep(step)
- && isSamplingRegionCovered(params, dimParams, grid, commRecord, multiSimComm));
+ (params.isCheckCoveringStep(step) && isSamplingRegionCovered(params, dimParams, grid));
}
/* The weighthistogram size after this update. */
histogramSize_.incrementNumUpdates();
}
-double BiasState::updateProbabilityWeightsAndConvolvedBias(const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- gmx::ArrayRef<const double> neighborLambdaEnergies,
+double BiasState::updateProbabilityWeightsAndConvolvedBias(ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ ArrayRef<const double> neighborLambdaEnergies,
std::vector<double, AlignedAllocator<double>>* weight) const
{
/* Only neighbors of the current coordinate value will have a non-negligible chance of getting sampled */
return std::log(weightSum);
}
-double BiasState::calcConvolvedBias(const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- const awh_dvec& coordValue) const
+double BiasState::calcConvolvedBias(ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ const awh_dvec& coordValue) const
{
int point = grid.nearestIndex(coordValue);
const GridPoint& gridPoint = grid.point(point);
gmx_bcast(sizeof(histogramSize_), &histogramSize_, commRecord->mpi_comm_mygroup);
}
-void BiasState::setFreeEnergyToConvolvedPmf(const std::vector<DimParams>& dimParams, const BiasGrid& grid)
+void BiasState::setFreeEnergyToConvolvedPmf(ArrayRef<const DimParams> dimParams, const BiasGrid& grid)
{
std::vector<float> convolvedPmf;
* \param[in] biasIndex The index of the bias.
* \param[in,out] pointState The state of the points in this bias.
*/
-static void readUserPmfAndTargetDistribution(const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- const std::string& filename,
- int numBias,
- int biasIndex,
- std::vector<PointState>* pointState)
+static void readUserPmfAndTargetDistribution(ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ const std::string& filename,
+ int numBias,
+ int biasIndex,
+ std::vector<PointState>* pointState)
{
/* Read the PMF and target distribution.
From the PMF, the convolved PMF, or the reference value free energy, can be calculated
}
}
-void BiasState::initGridPointState(const AwhBiasParams& awhBiasParams,
- const std::vector<DimParams>& dimParams,
- const BiasGrid& grid,
- const BiasParams& params,
- const std::string& filename,
- int numBias)
+void BiasState::initGridPointState(const AwhBiasParams& awhBiasParams,
+ ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ const BiasParams& params,
+ const std::string& filename,
+ int numBias)
{
/* Modify PMF, free energy and the constant target distribution factor
* to user input values if there is data given.
*/
- if (awhBiasParams.bUserData)
+ if (awhBiasParams.userPMFEstimate())
{
readUserPmfAndTargetDistribution(dimParams, grid, filename, numBias, params.biasIndex, &points_);
setFreeEnergyToConvolvedPmf(dimParams, grid);
}
/* The local Boltzmann distribution is special because the target distribution is updated as a function of the reference weighthistogram. */
- GMX_RELEASE_ASSERT(params.eTarget != eawhtargetLOCALBOLTZMANN || points_[0].weightSumRef() != 0,
+ GMX_RELEASE_ASSERT(params.eTarget != AwhTargetType::LocalBoltzmann || points_[0].weightSumRef() != 0,
"AWH reference weight histogram not initialized properly with local "
"Boltzmann target distribution.");
normalizePmf(params.numSharedUpdate);
}
-BiasState::BiasState(const AwhBiasParams& awhBiasParams,
- double histogramSizeInitial,
- const std::vector<DimParams>& dimParams,
- const BiasGrid& grid) :
+BiasState::BiasState(const AwhBiasParams& awhBiasParams,
+ double histogramSizeInitial,
+ ArrayRef<const DimParams> dimParams,
+ const BiasGrid& grid,
+ const BiasSharing* biasSharing) :
coordState_(awhBiasParams, dimParams, grid),
points_(grid.numPoints()),
weightSumCovering_(grid.numPoints()),
- histogramSize_(awhBiasParams, histogramSizeInitial)
+ histogramSize_(awhBiasParams, histogramSizeInitial),
+ biasSharing_(biasSharing)
{
/* The minimum and maximum multidimensional point indices that are affected by the next update */
for (size_t d = 0; d < dimParams.size(); d++)
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff