#include "bias.h"
#include "biassharing.h"
+#include "correlationgrid.h"
#include "pointstate.h"
namespace gmx
Bias::ThisRankWillDoIO thisRankWillDoIO = (MASTER(commRecord_) ? Bias::ThisRankWillDoIO::Yes : Bias::ThisRankWillDoIO::No);
biasCoupledToSystem_.emplace_back(Bias(k, awhParams, awhParams.awhBiasParams[k], dimParams, beta, inputRecord.delta_t, numSharingSimulations, biasInitFilename, thisRankWillDoIO),
pullCoordIndex);
+
+ biasCoupledToSystem_.back().bias.printInitializationToLog(fplog);
}
/* Need to register the AWH coordinates to be allowed to apply forces to the pull coordinates. */
* sampling observables based on the input pull coordinate value,
* setting the bias force and/or updating the AWH bias state.
*/
- awh_dvec biasForce;
- double biasPotential;
- double biasPotentialJump;
+ double biasPotential;
+ double biasPotentialJump;
/* Note: In the near future this call will be split in calls
* to supports bias sharing within a single simulation.
*/
- biasCts.bias.calcForceAndUpdateBias(coordValue, biasForce,
- &biasPotential, &biasPotentialJump,
- commRecord_->ms,
- t, step, seed_, fplog);
+ gmx::ArrayRef<const double> biasForce =
+ biasCts.bias.calcForceAndUpdateBias(coordValue,
+ &biasPotential, &biasPotentialJump,
+ commRecord_->ms,
+ t, step, seed_, fplog);
awhPotential += biasPotential;
if (isOutputStep(step))
{
- /* Ensure that we output fully updated data */
+ /* We might have skipped updates for part of the grid points.
+ * Ensure all points are updated before writing out their data.
+ */
biasCts.bias.doSkippedUpdatesForAllPoints();
}
}
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/stringutil.h"
+#include "correlationgrid.h"
+#include "correlationhistory.h"
#include "pointstate.h"
namespace gmx
}
}
-void Bias::calcForceAndUpdateBias(const awh_dvec coordValue,
- awh_dvec biasForce,
- double *awhPotential,
- double *potentialJump,
- const gmx_multisim_t *ms,
- double t,
- gmx_int64_t step,
- gmx_int64_t seed,
- FILE *fplog)
+gmx::ArrayRef<const double>
+Bias::calcForceAndUpdateBias(const awh_dvec coordValue,
+ double *awhPotential,
+ double *potentialJump,
+ const gmx_multisim_t *ms,
+ double t,
+ gmx_int64_t step,
+ gmx_int64_t seed,
+ FILE *fplog)
{
if (step < 0)
{
if (sampleCoord)
{
+ updateForceCorrelationGrid(probWeightNeighbor, t);
+
state_.sampleCoordAndPmf(grid_, probWeightNeighbor, convolvedBias);
}
}
if (params_.convolveForce)
{
state_.calcConvolvedForce(dimParams_, grid_, probWeightNeighbor,
- biasForce);
+ tempForce_, biasForce_);
potential = -convolvedBias*params_.invBeta;
}
GMX_RELEASE_ASSERT(state_.points()[coordState.umbrellaGridpoint()].inTargetRegion(),
"AWH bias grid point for the umbrella reference value is outside of the target region.");
potential =
- state_.calcUmbrellaForceAndPotential(dimParams_, grid_, coordState.umbrellaGridpoint(), biasForce);
+ state_.calcUmbrellaForceAndPotential(dimParams_, grid_, coordState.umbrellaGridpoint(), biasForce_);
/* Moving the umbrella results in a force correction and
* a new potential. The umbrella center is sampled as often as
*/
if (moveUmbrella)
{
- double newPotential = state_.moveUmbrella(dimParams_, grid_, probWeightNeighbor, biasForce, step, seed, params_.biasIndex);
+ double newPotential = state_.moveUmbrella(dimParams_, grid_, probWeightNeighbor, biasForce_, step, seed, params_.biasIndex);
*potentialJump = newPotential - potential;
}
}
/* Check the sampled histograms and potentially warn user if something is suspicious */
warnForHistogramAnomalies(t, step, fplog);
+
+ return biasForce_;
}
void Bias::restoreStateFromHistory(const AwhBiasHistory *biasHistory,
{
GMX_RELEASE_ASSERT(biasHistory != nullptr, "On the master rank we need a valid history object to restore from");
state_.restoreFromHistory(*biasHistory, grid_);
+
+ if (forceCorrelationGrid_ != nullptr)
+ {
+ forceCorrelationGrid_->restoreStateFromHistory(biasHistory->forceCorrelationGrid);
+ }
}
if (PAR(cr))
GMX_RELEASE_ASSERT(biasHistory != nullptr, "Need a valid biasHistory");
state_.initHistoryFromState(biasHistory);
+
+ if (forceCorrelationGrid_ != nullptr)
+ {
+ biasHistory->forceCorrelationGrid = initCorrelationGridHistoryFromState(forceCorrelationGrid());
+ }
}
void Bias::updateHistory(AwhBiasHistory *biasHistory) const
GMX_RELEASE_ASSERT(biasHistory != nullptr, "Need a valid biasHistory");
state_.updateHistory(biasHistory, grid_);
+
+ if (forceCorrelationGrid_ != nullptr)
+ {
+ updateCorrelationGridHistory(&biasHistory->forceCorrelationGrid, forceCorrelationGrid());
+ }
}
Bias::Bias(int biasIndexInCollection,
params_(awhParams, awhBiasParams, dimParams_, beta, mdTimeStep, disableUpdateSkips, numSharingSimulations, grid_.axis(), biasIndexInCollection),
state_(awhBiasParams, params_.initialHistogramSize, dimParams_, grid_),
thisRankDoesIO_(thisRankWillDoIO == ThisRankWillDoIO::Yes),
+ biasForce_(ndim()),
alignedTempWorkSpace_(),
+ tempForce_(ndim()),
numWarningsIssued_(0)
{
/* For a global update updateList covers all points, so reserve that */
if (thisRankDoesIO_)
{
+ /* Set up the force correlation object. */
+
+ /* We let the correlation init function set its parameters
+ * to something useful for now.
+ */
+ double blockLength = 0;
+ /* Construct the force correlation object. */
+ forceCorrelationGrid_ =
+ gmx::compat::make_unique<CorrelationGrid>(state_.points().size(), ndim(),
+ blockLength, CorrelationGrid::BlockLengthMeasure::Time,
+ awhParams.nstSampleCoord*mdTimeStep);
+
writer_ = std::unique_ptr<BiasWriter>(new BiasWriter(*this));
}
}
+void Bias::printInitializationToLog(FILE *fplog) const
+{
+ if (fplog != nullptr && forceCorrelationGrid_ != nullptr)
+ {
+ std::string prefix =
+ gmx::formatString("\nawh%d:", params_.biasIndex + 1);
+
+ fprintf(fplog,
+ "%s initial force correlation block length = %g %s"
+ "%s force correlation number of blocks = %d",
+ prefix.c_str(), forceCorrelationGrid().getBlockLength(),
+ forceCorrelationGrid().blockLengthMeasure == CorrelationGrid::BlockLengthMeasure::Weight ? "" : "ps",
+ prefix.c_str(), forceCorrelationGrid().getNumBlocks());
+ }
+}
+
+void Bias::updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor,
+ double t)
+{
+ if (forceCorrelationGrid_ == nullptr)
+ {
+ return;
+ }
+
+ const std::vector<int> &neighbor = grid_.point(state_.coordState().gridpointIndex()).neighbor;
+
+ gmx::ArrayRef<double> forceFromNeighbor = tempForce_;
+ for (size_t n = 0; n < neighbor.size(); n++)
+ {
+ double weightNeighbor = probWeightNeighbor[n];
+ int indexNeighbor = neighbor[n];
+
+ /* Add the force data of this neighbor point. Note: the sum of these forces is the convolved force.
+
+ We actually add the force normalized by beta which has the units of 1/length. This means that the
+ resulting correlation time integral is directly in units of friction time/length^2 which is really what
+ we're interested in. */
+ state_.calcUmbrellaForceAndPotential(dimParams_, grid_, indexNeighbor, forceFromNeighbor);
+
+ /* Note: we might want to give a whole list of data to add instead and have this loop in the data adding function */
+ forceCorrelationGrid_->addData(indexNeighbor, weightNeighbor, forceFromNeighbor, t);
+ }
+}
+
/* Return the number of data blocks that have been prepared for writing. */
int Bias::numEnergySubblocksToWrite() const
{
#include <memory>
#include <vector>
+#include "gromacs/compat/make_unique.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/utility/alignedallocator.h"
#include "gromacs/utility/basedefinitions.h"
struct AwhHistory;
struct AwhParams;
struct AwhPointStateHistory;
+class CorrelationGrid;
class Grid;
class GridAxis;
class PointState;
ThisRankWillDoIO thisRankWillDoIO,
BiasParams::DisableUpdateSkips disableUpdateSkips = BiasParams::DisableUpdateSkips::no);
+ /*! \brief
+ * Print information about initialization to log file.
+ *
+ * Prints information about AWH variables that are set internally
+ * but might be of interest to the user.
+ *
+ * \param[in,out] fplog Log file, can be nullptr.
+ */
+ void printInitializationToLog(FILE *fplog) const;
+
/*! \brief
* Evolves the bias at every step.
*
* - reweight samples to extract the PMF.
*
* \param[in] coordValue The current coordinate value(s).
- * \param[out] biasForce The bias force.
* \param[out] awhPotential Bias potential.
* \param[out] potentialJump Change in bias potential for this bias.
* \param[in] ms Struct for multi-simulation communication.
* \param[in] step Time step.
* \param[in] seed Random seed.
* \param[in,out] fplog Log file.
+ * \returns a reference to the bias force, size \ref ndim(), valid until the next call of this method or destruction of Bias, whichever comes first.
*/
- void calcForceAndUpdateBias(const awh_dvec coordValue,
- awh_dvec biasForce,
- double *awhPotential,
- double *potentialJump,
- const gmx_multisim_t *ms,
- double t,
- gmx_int64_t step,
- gmx_int64_t seed,
- FILE *fplog);
+ gmx::ArrayRef<const double>
+ calcForceAndUpdateBias(const awh_dvec coordValue,
+ double *awhPotential,
+ double *potentialJump,
+ const gmx_multisim_t *ms,
+ double t,
+ gmx_int64_t step,
+ gmx_int64_t seed,
+ FILE *fplog);
/*! \brief
* Calculates the convolved bias for a given coordinate value.
gmx_int64_t step,
FILE *fplog);
+ /*! \brief
+ * Collect samples for the force correlation analysis on the grid.
+ *
+ * \param[in] probWeightNeighbor Probability weight of the neighboring points.
+ * \param[in] t The time.
+ */
+ void updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor,
+ double t);
+
public:
+ /*! \brief Return a const reference to the force correlation grid.
+ */
+ const CorrelationGrid &forceCorrelationGrid() const
+ {
+ GMX_RELEASE_ASSERT(forceCorrelationGrid_ != nullptr, "forceCorrelationGrid() should only be called with a valid force correlation object");
+
+ return *forceCorrelationGrid_.get();
+ }
+
/*! \brief Return the number of data blocks that have been prepared for writing.
*/
int numEnergySubblocksToWrite() const;
const bool thisRankDoesIO_; /**< Tells whether this MPI rank will do I/O (checkpointing, AWH output) */
+ std::vector<double> biasForce_; /**< Vector for returning the force to the caller. */
+
+ /* Force correlation grid */
+ std::unique_ptr<CorrelationGrid> forceCorrelationGrid_; /**< Takes care of force correlation statistics for every grid point. */
+
/* I/O */
std::unique_ptr<BiasWriter> writer_; /**< Takes care of AWH data output. */
- /* Temporary working vector used during the update.
- * This only here to avoid allocation at every MD step.
+ /* Temporary working vectors used during the update.
+ * These are only here to avoid allocation at every MD step.
*/
std::vector < double, AlignedAllocator < double>> alignedTempWorkSpace_; /**< Working vector of doubles. */
+ std::vector<double> tempForce_; /**< Bias force work buffer. */
/* Run-local counter to avoid flooding log with warnings. */
int numWarningsIssued_; /**< The number of warning issued in the current run. */
double BiasState::calcUmbrellaForceAndPotential(const std::vector<DimParams> &dimParams,
const Grid &grid,
int point,
- awh_dvec force) const
+ gmx::ArrayRef<double> force) const
{
double potential = 0;
for (size_t d = 0; d < dimParams.size(); d++)
void BiasState::calcConvolvedForce(const std::vector<DimParams> &dimParams,
const Grid &grid,
gmx::ArrayRef<const double> probWeightNeighbor,
- awh_dvec force) const
+ gmx::ArrayRef<double> forceWorkBuffer,
+ gmx::ArrayRef<double> force) const
{
for (size_t d = 0; d < dimParams.size(); d++)
{
}
/* Only neighboring points have non-negligible contribution. */
- const std::vector<int> &neighbor = grid.point(coordState_.gridpointIndex()).neighbor;
+ const std::vector<int> &neighbor = grid.point(coordState_.gridpointIndex()).neighbor;
+ gmx::ArrayRef<double> forceFromNeighbor = forceWorkBuffer;
for (size_t n = 0; n < neighbor.size(); n++)
{
double weightNeighbor = probWeightNeighbor[n];
int indexNeighbor = neighbor[n];
/* Get the umbrella force from this point. The returned potential is ignored here. */
- awh_dvec forceFromNeighbor;
calcUmbrellaForceAndPotential(dimParams, grid, indexNeighbor,
forceFromNeighbor);
double BiasState::moveUmbrella(const std::vector<DimParams> &dimParams,
const Grid &grid,
gmx::ArrayRef<const double> probWeightNeighbor,
- awh_dvec biasForce,
+ gmx::ArrayRef<double> biasForce,
gmx_int64_t step,
gmx_int64_t seed,
int indexSeed)
/* Generate and set a new coordinate reference value */
coordState_.sampleUmbrellaGridpoint(grid, coordState_.gridpointIndex(), probWeightNeighbor, step, seed, indexSeed);
- awh_dvec newForce;
- double newPotential =
+ std::vector<double> newForce(dimParams.size());
+ double newPotential =
calcUmbrellaForceAndPotential(dimParams, grid, coordState_.umbrellaGridpoint(), newForce);
/* A modification of the reference value at time t will lead to a different
at steps when the reference value has been moved. E.g. if the ref. value
is set every step to (coord dvalue +/- delta) would give zero force.
*/
- for (size_t d = 0; d < dimParams.size(); d++)
+ for (size_t d = 0; d < biasForce.size(); d++)
{
- /* clang thinks newForce[d] can be garbage */
- #ifndef __clang_analyzer__
/* Average of the current and new force */
biasForce[d] = 0.5*(biasForce[d] + newForce[d]);
- #endif
}
return newPotential;
* Allocate and initialize a bias history with the given bias state.
*
* This function will be called at the start of a new simulation.
- * Note that this only sets the correct size and does produces
- * a history object with all variables zero.
- * History data is set by \ref updateHistory.
+ * Note that this only sets the correct size and does produce
+ * a valid history object, but with all data set to zero.
+ * Actual history data is set by \ref updateHistory.
*
* \param[in,out] biasHistory AWH history to initialize.
*/
void initHistoryFromState(AwhBiasHistory *biasHistory) const;
/*! \brief
- * Update the bias history with a new state.
+ * Update the bias state history with the current state.
*
* \param[out] biasHistory Bias history struct.
* \param[in] grid The bias grid.
double calcUmbrellaForceAndPotential(const std::vector<DimParams> &dimParams,
const Grid &grid,
int point,
- awh_dvec force) const;
+ gmx::ArrayRef<double> force) const;
/*! \brief
* Calculates and sets the convolved force acting on the coordinate.
* \param[in] dimParams The bias dimensions parameters.
* \param[in] grid The grid.
* \param[in] probWeightNeighbor Probability weights of the neighbors.
+ * \param[in] forceWorkBuffer Force work buffer, values only used internally.
* \param[in,out] force Bias force vector to set.
*/
void calcConvolvedForce(const std::vector<DimParams> &dimParams,
const Grid &grid,
gmx::ArrayRef<const double> probWeightNeighbor,
- awh_dvec force) const;
+ gmx::ArrayRef<double> forceWorkBuffer,
+ gmx::ArrayRef<double> force) const;
/*! \brief
* Move the center point of the umbrella potential.
double moveUmbrella(const std::vector<DimParams> &dimParams,
const Grid &grid,
gmx::ArrayRef<const double> probWeightNeighbor,
- awh_dvec biasForce,
+ gmx::ArrayRef<double> biasForce,
gmx_int64_t step,
gmx_int64_t seed,
int indexSeed);
HistogramSize histogramSize_; /**< Global histogram size related values. */
/* Track the part of the grid sampled since the last update. */
- awh_ivec originUpdatelist_; /**< The origin of the rectangular region that has been sampled since last update. */
- awh_ivec endUpdatelist_; /**< The end of the rectangular region that has been sampled since last update. */
+ awh_ivec originUpdatelist_; /**< The origin of the rectangular region that has been sampled since last update. */
+ awh_ivec endUpdatelist_; /**< The end of the rectangular region that has been sampled since last update. */
};
//! Linewidth used for warning output
#include "gromacs/utility/smalloc.h"
#include "bias.h"
+#include "correlationgrid.h"
#include "grid.h"
#include "pointstate.h"
*/
static const std::map<AwhOutputEntryType, Normalization> outputTypeToNormalization =
{
- { AwhOutputEntryType::MetaData, Normalization::None },
- { AwhOutputEntryType::CoordValue, Normalization::Coordinate },
- { AwhOutputEntryType::Pmf, Normalization::FreeEnergy },
- { AwhOutputEntryType::Bias, Normalization::FreeEnergy },
- { AwhOutputEntryType::Visits, Normalization::Distribution },
- { AwhOutputEntryType::Weights, Normalization::Distribution },
- { AwhOutputEntryType::Target, Normalization::Distribution }
+ { AwhOutputEntryType::MetaData, Normalization::None },
+ { AwhOutputEntryType::CoordValue, Normalization::Coordinate },
+ { AwhOutputEntryType::Pmf, Normalization::FreeEnergy },
+ { AwhOutputEntryType::Bias, Normalization::FreeEnergy },
+ { AwhOutputEntryType::Visits, Normalization::Distribution },
+ { AwhOutputEntryType::Weights, Normalization::Distribution },
+ { AwhOutputEntryType::Target, Normalization::Distribution },
+ { AwhOutputEntryType::ForceCorrelationVolume, Normalization::Distribution },
+ { AwhOutputEntryType::FrictionTensor, Normalization::None }
};
/*! \brief
case AwhOutputEntryType::Target:
normalizationValue = static_cast<float>(bias.state().points().size());
break;
+ case AwhOutputEntryType::ForceCorrelationVolume:
+ normalizationValue = static_cast<double>(bias.state().points().size());
+ break;
default:
break;
}
{
outputTypeNumBlock[outputType] = bias.ndim();
}
+ else if (outputType == AwhOutputEntryType::FrictionTensor)
+ {
+ outputTypeNumBlock[outputType] = bias.forceCorrelationGrid().tensorSize();
+ }
else
{
/* Most output variable types need one block */
int blockStart = getVarStartBlock(outputType);
GMX_ASSERT(pointIndex < static_cast<int>(block_[blockStart].data().size()), "Attempt to transfer AWH data to block for point index out of range");
+ const CorrelationGrid &forceCorrelation = bias.forceCorrelationGrid();
+ int numCorrelation = forceCorrelation.tensorSize();
+
/* Transfer the point data of this variable to the right block(s) */
int b = blockStart;
switch (outputType)
case AwhOutputEntryType::Target:
block_[b].data()[pointIndex] = bias.state().points()[pointIndex].target();
break;
+ case AwhOutputEntryType::ForceCorrelationVolume:
+ block_[b].data()[pointIndex] = forceCorrelation.tensors()[pointIndex].getVolumeElement(forceCorrelation.dtSample);
+ break;
+ case AwhOutputEntryType::FrictionTensor:
+ /* Store force correlation in units of friction, i.e. time/length^2 */
+ for (int n = 0; n < numCorrelation; n++)
+ {
+ block_[b].data()[pointIndex] = forceCorrelation.tensors()[pointIndex].getTimeIntegral(n, forceCorrelation.dtSample);
+ b++;
+ }
+ break;
default:
GMX_RELEASE_ASSERT(false, "Unknown AWH output variable");
break;
//! Enum with the AWH variables to write.
enum class AwhOutputEntryType
{
- MetaData, //!< Meta data.
- CoordValue, //!< Coordinate value.
- Pmf, //!< The pmf.
- Bias, //!< The bias.
- Visits, //!< The number of visits.
- Weights, //!< The weights.
- Target //!< The target distribition.
+ MetaData, //!< Meta data.
+ CoordValue, //!< Coordinate value.
+ Pmf, //!< The pmf.
+ Bias, //!< The bias.
+ Visits, //!< The number of visits.
+ Weights, //!< The weights.
+ Target, //!< The target distribition.
+ ForceCorrelationVolume, //!< The volume of the force correlation tensor.
+ FrictionTensor //!< The full friction tensor.
};
//! Enum with the types of metadata to write.
/*! \brief Constructor
*
* \param[in] numPoints Number of points in block.
- * \param[in] normalizationType Value for normalization type enum.
- * \param[in] normalizationValue Normalization value.
+ * \param[in] normalizationType Type of normalization.
+ * \param[in] normalizationValue Value to normalization with.
*/
AwhEnergyBlock(int numPoints,
Normalization normalizationType,
*/
class BiasWriter
{
+ public:
+ /*! \brief Constructor.
+ *
+ * \param[in] bias The AWH bias.
+ */
+ BiasWriter(const Bias &bias);
+
+ /*! \brief Returns the number of data blocks.
+ *
+ * \returns the number of data blocks.
+ */
+ int numBlocks() const
+ {
+ return block_.size();
+ }
+
+ /*! \brief Collect AWH bias data in blocks and write to energy subblocks.
+ *
+ * \param[in] bias The AWH Bias.
+ * \param[in,out] subblock Energy subblocks to write to.
+ * \returns the number of blocks written.
+ */
+ int writeToEnergySubblocks(const Bias &bias, t_enxsubblock *subblock);
+
private:
/*! \brief Query if the writer has a block for the given variable.
*
const Bias &bias,
gmx::ArrayRef<const float> pmf);
- public:
- /*! \brief Constructor.
- *
- * \param[in] bias The AWH bias.
- */
- BiasWriter(const Bias &bias);
-
- /*! \brief Returns the number of data blocks.
- *
- * \returns the number of data blocks.
- */
- int numBlocks() const
- {
- return block_.size();
- }
-
- private:
/*! \brief
* Prepare the bias output data.
*
*/
void prepareBiasOutput(const Bias &bias);
- public:
- /*! \brief Collect AWH bias data in blocks and write to energy subblocks.
- *
- * \param[in] bias The AWH Bias.
- * \param[in,out] subblock Energy subblocks to write to.
- * \returns the number of blocks written.
- */
- int writeToEnergySubblocks(const Bias &bias, t_enxsubblock *subblock);
-
private:
std::vector<AwhEnergyBlock> block_; /**< The data blocks */
std::map<AwhOutputEntryType, int> outputTypeToBlock_; /**< Start block index for each variable, -1 when variable should not be written */
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015,2016,2017, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Implements the CorrelationGrid class to collect correlation statistics on a grid, using several block lengths.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_awh
+ */
+
+#include "gmxpre.h"
+
+#include "correlationgrid.h"
+
+#include "gromacs/math/utilities.h"
+#include "gromacs/utility/gmxassert.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+/*! \brief
+ * Return the number of block data structs needed for keeping a certain number of blocks.
+ *
+ * The list start with 1 block and doubles, so we need 1 + 2log(numBlocks).
+ *
+ * \param[in] numBlocks Number of blocks.
+ * \returns the number of block data structs.
+ */
+int getBlockDataListSize(int numBlocks)
+{
+ int blockDataListSize = 1;
+
+ while (numBlocks > (1 << (blockDataListSize - 1)))
+ {
+ blockDataListSize++;
+ }
+
+ GMX_RELEASE_ASSERT((1 << (blockDataListSize - 1)) == numBlocks, "numBlocks should be a power of 2");
+
+ return blockDataListSize;
+}
+
+} // namespace
+
+CorrelationGrid::CorrelationGrid(int numPoints,
+ int numDim,
+ double blockLengthInit,
+ BlockLengthMeasure blockLengthMeasure,
+ double dtSample) :
+ dtSample(dtSample),
+ blockLengthMeasure(blockLengthMeasure)
+{
+ /* Set the initial block length for the block averaging. The length doesn't really matter
+ after the block length has been doubled a few times, as long as it's set small enough */
+ if (blockLengthMeasure == BlockLengthMeasure::Weight)
+ {
+ blockLengthInit = blockLengthInit > 0 ? blockLengthInit : 1;
+ }
+ else
+ {
+ blockLengthInit = blockLengthInit > 0 ? blockLengthInit : dtSample;
+ }
+
+ /* Set the number of blocks. The number of blocks determines the current span of the data
+ and how many different block lengths (nblockdata) we need to keep track of to be able to
+ increase the block length later */
+ int numBlocks = CorrelationTensor::c_numCorrelationBlocks;
+ int BlockDataListSize = getBlockDataListSize(numBlocks);
+
+ tensors_.resize(numPoints, CorrelationTensor(numDim, BlockDataListSize, blockLengthInit));
+}
+
+int CorrelationGrid::getNumBlocks() const
+{
+ const auto &blockDataList = tensors()[0].blockDataList();
+ double maxBlockLength = blockDataList.back().blockLength();
+ double minBlockLength = blockDataList[0].blockLength();
+
+ /* If we have a finite block span we have a constant number of blocks, otherwise we are always adding more blocks (and we don't keep track of the number) */
+ if (maxBlockLength < GMX_DOUBLE_MAX)
+ {
+ return static_cast<int>(maxBlockLength/minBlockLength);
+ }
+ else
+ {
+ return -1;
+ }
+}
+
+double CorrelationGrid::getBlockLength() const
+{
+ /* Return the minimum blocklength */
+ return tensors()[0].blockDataList()[0].blockLength();
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015,2016,2017, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief
+ * Declares the CorrelationGrid class to collect correlation statistics on a grid, using several block lengths.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_awh
+ */
+
+#ifndef GMX_AWH_CORRELATIONGRID_H
+#define GMX_AWH_CORRELATIONGRID_H
+
+#include <cstddef>
+
+#include <vector>
+
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/gmxassert.h"
+
+#include "correlationtensor.h"
+
+namespace gmx
+{
+
+struct CorrelationGridHistory;
+
+/*! \internal
+ * \brief Grid of local correlation tensors.
+ *
+ * This class provides the means for a bias to interaction with the grid
+ * of correlation tensors. The grid should have the same number of points
+ * and the same dimensionality as the bias grid.
+ */
+class CorrelationGrid
+{
+ public:
+ //! Enum that sets how we measure block length.
+ enum class BlockLengthMeasure
+ {
+ Time, //!< Measure block length in time.
+ Weight //!< Measure block length in sampled weight.
+ };
+
+ /*! \brief Constructor.
+ *
+ * \param[in] numPoints Number of points in the grid.
+ * \param[in] numDims Number of dimensions of the grid.
+ * \param[in] blockLengthInit Initial length of the blocks used for block averaging.
+ * \param[in] blockLengthMeasure Sets how we measure block length.
+ * \param[in] dtSample Time step for sampling correlations.
+ */
+ CorrelationGrid(int numPoints,
+ int numDims,
+ double blockLengthInit,
+ BlockLengthMeasure blockLengthMeasure,
+ double dtSample);
+
+ /*! \brief Adds a weighted data vector to one point in the correlation grid.
+ *
+ * \param[in] pointIndex Index of the point to add data to.
+ * \param[in] weight Weight to assign to the data.
+ * \param[in] data One data point for each grid dimension.
+ * \param[in] t The time when the data was sampled.
+ */
+ void addData(int pointIndex,
+ double weight,
+ gmx::ArrayRef<const double> data,
+ double t)
+ {
+ tensors_[pointIndex].addData(weight, data, blockLengthMeasure == BlockLengthMeasure::Weight, t);
+ };
+
+ /*! \brief Restores the correlation grid state from the correlation grid history.
+ *
+ * The setup in the history should match that of this simulation.
+ * If this is not the case, an exception is thrown.
+ *
+ * \param[in] correlationGridHistory The correlation grid state history.
+ */
+ void restoreStateFromHistory(const CorrelationGridHistory &correlationGridHistory);
+
+ /*! \brief Returns the number of elements in the tensor: dim*(dim+1)/2.
+ */
+ int tensorSize() const
+ {
+ GMX_RELEASE_ASSERT(tensors_.size() > 0, "Should only call tensorSize on a valid grid");
+
+ return tensors_[0].blockDataList()[0].correlationIntegral().size();
+ }
+
+ /*! \brief Returns the size of the block data list.
+ */
+ int blockDataListSize() const
+ {
+ GMX_RELEASE_ASSERT(tensors_.size() > 0, "Should only call tensorSize on a valid grid");
+
+ return tensors_[0].blockDataList().size();
+ }
+
+ /*! \brief Get a const reference to the correlation grid data.
+ */
+ const std::vector<CorrelationTensor> &tensors() const
+ {
+ return tensors_;
+ }
+
+ /* Right now the below functions are only used for an initial log printing. */
+
+ /*! \brief Get the current blocklength.
+ */
+ double getBlockLength() const;
+
+ /*! \brief Get the current number of blocks.
+ *
+ * If we have a finite block span we have a constant number of blocks,
+ * otherwise we are always adding more blocks (and we don't keep
+ * track of the number), so we return -1.
+ */
+ int getNumBlocks() const;
+
+ public:
+ const double dtSample; /**< Time in between samples. */
+ const BlockLengthMeasure blockLengthMeasure; /**< The measure for the block length. */
+ private:
+ std::vector<CorrelationTensor> tensors_; /**< Correlation tensor grid */
+};
+
+} // namespace gmx
+
+#endif /* GMX_AWH_CORRELATIONGRID_H */
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015,2016,2017, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief
+ * Implements helper functions for checkpointing the AWH state and observables history.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_awh
+ */
+
+#include "gmxpre.h"
+
+#include "correlationhistory.h"
+
+#include <assert.h>
+
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/mdtypes/awh-correlation-history.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/gmxassert.h"
+
+#include "correlationgrid.h"
+
+namespace gmx
+{
+
+void initCorrelationGridHistory(CorrelationGridHistory *correlationGridHistory,
+ int numCorrelationTensors,
+ int tensorSize,
+ int blockDataListSize)
+{
+ correlationGridHistory->numCorrelationTensors = numCorrelationTensors;
+ correlationGridHistory->tensorSize = tensorSize;
+ correlationGridHistory->blockDataListSize = blockDataListSize;
+
+ correlationGridHistory->blockDataBuffer.resize(numCorrelationTensors*tensorSize*blockDataListSize);
+}
+
+CorrelationGridHistory initCorrelationGridHistoryFromState(const CorrelationGrid &correlationGrid)
+{
+ CorrelationGridHistory correlationGridHistory;
+
+ initCorrelationGridHistory(&correlationGridHistory, correlationGrid.tensors().size(), correlationGrid.tensorSize(), correlationGrid.blockDataListSize());
+
+ return correlationGridHistory;
+}
+
+/* Update the correlation grid history for checkpointing. */
+void updateCorrelationGridHistory(CorrelationGridHistory *correlationGridHistory,
+ const CorrelationGrid &correlationGrid)
+{
+ GMX_RELEASE_ASSERT(correlationGridHistory != nullptr, "We need a valid history object");
+
+ gmx::ArrayRef<CorrelationBlockDataHistory> blockDataBuffer = correlationGridHistory->blockDataBuffer;
+
+ /* Store the grid in a linear array */
+ size_t bufferIndex = 0;
+ for (const CorrelationTensor &tensor : correlationGrid.tensors())
+ {
+ const int numDims = tensor.blockDataList()[0].coordData().size();
+ const int tensorSize = tensor.blockDataList()[0].correlationIntegral().size();
+
+ /* Loop of the tensor elements, ignore the symmetric data */
+ int d1 = 0;
+ int d2 = 0;
+ for (int k = 0; k < tensorSize; k++)
+ {
+ /* BlockData for each correlation element */
+ for (const CorrelationBlockData &blockData : tensor.blockDataList())
+ {
+ const CorrelationBlockData::CoordData &cx = blockData.coordData()[d1];
+ const CorrelationBlockData::CoordData &cy = blockData.coordData()[d2];
+
+ CorrelationBlockDataHistory &bdh = blockDataBuffer[bufferIndex];
+
+ bdh.blockSumWeight = blockData.blockSumWeight();
+ bdh.blockSumSquareWeight = blockData.blockSumSquareWeight();
+ bdh.blockSumWeightX = cx.blockSumWeightX;
+ bdh.blockSumWeightY = cy.blockSumWeightX;
+ bdh.sumOverBlocksSquareBlockWeight = blockData.sumOverBlocksSquareBlockWeight();
+ bdh.sumOverBlocksBlockSquareWeight = blockData.sumOverBlocksBlockSquareWeight();
+ bdh.sumOverBlocksBlockWeightBlockWeightX = cx.sumOverBlocksBlockWeightBlockWeightX;
+ bdh.sumOverBlocksBlockWeightBlockWeightY = cy.sumOverBlocksBlockWeightBlockWeightX;
+ bdh.previousBlockIndex = blockData.previousBlockIndex();
+ bdh.blockLength = blockData.blockLength();
+ bdh.correlationIntegral = blockData.correlationIntegral()[k];
+
+ bufferIndex++;
+ }
+
+ d1++;
+ if (d1 == numDims)
+ {
+ d1 = 0;
+ d2++;
+ }
+ }
+ }
+
+ GMX_RELEASE_ASSERT(bufferIndex == blockDataBuffer.size(), "We should store exactly as many elements as the buffer size");
+}
+
+void CorrelationBlockData::restoreFromHistory(const CorrelationBlockDataHistory &blockHistory,
+ const std::vector<CorrelationBlockData::CoordData> &coordData,
+ const std::vector<double> &correlationIntegral)
+{
+ blockSumWeight_ = blockHistory.blockSumWeight;
+ blockSumSquareWeight_ = blockHistory.blockSumSquareWeight;
+ sumOverBlocksSquareBlockWeight_ = blockHistory.sumOverBlocksSquareBlockWeight;
+ sumOverBlocksBlockSquareWeight_ = blockHistory.sumOverBlocksBlockSquareWeight;
+ previousBlockIndex_ = blockHistory.previousBlockIndex;
+ blockLength_ = blockHistory.blockLength;
+ coordData_ = coordData;
+ correlationIntegral_ = correlationIntegral;
+}
+
+/* Restore a correlation element from history. */
+void CorrelationTensor::restoreFromHistory(const std::vector<CorrelationBlockDataHistory> &blockDataBuffer,
+ size_t *bufferIndex)
+{
+ /* Blockdata for each correlation element */
+ for (CorrelationBlockData &blockData : blockDataList_)
+ {
+ /* Correlation elements for each tensor */
+ const int numDims = blockDataList_[0].coordData().size();
+ const int tensorSize = blockDataList_[0].correlationIntegral().size();
+ int d1 = 0;
+ int d2 = 0;
+
+ /* Temporary containers to collect data */
+ std::vector<CorrelationBlockData::CoordData> coordData(numDims);
+ std::vector<double> correlationIntegral(tensorSize);
+ for (int k = 0; k < tensorSize; k++)
+ {
+ if (*bufferIndex >= blockDataBuffer.size())
+ {
+ GMX_THROW(InvalidInputError("Mismatch of the correlation tensor size for the force correlation between checkpoint and simulation. Likely you have provided a checkpoint from a different simulation."));
+ }
+ const CorrelationBlockDataHistory &blockHistory = blockDataBuffer[*bufferIndex];
+
+ /* To simplify the checkpointing, CorrelationBlockDataHistory
+ * duplicates some weight data for all tensor elements.
+ * Here we collect the coordinate and tensor data
+ * in temporary buffers.
+ */
+ coordData[d1].blockSumWeightX = blockHistory.blockSumWeightX;
+ coordData[d2].blockSumWeightX = blockHistory.blockSumWeightY;
+ coordData[d1].sumOverBlocksBlockWeightBlockWeightX = blockHistory.sumOverBlocksBlockWeightBlockWeightX;
+ coordData[d2].sumOverBlocksBlockWeightBlockWeightX = blockHistory.sumOverBlocksBlockWeightBlockWeightY;
+
+ correlationIntegral[k] = blockHistory.correlationIntegral;
+
+ /* Check if we collected all data needed for blockData */
+ if (k == tensorSize - 1)
+ {
+ blockData.restoreFromHistory(blockHistory, coordData, correlationIntegral);
+ }
+
+ (*bufferIndex)++;
+
+ d1++;
+ if (d1 == numDims)
+ {
+ d1 = 0;
+ d2++;
+ }
+ }
+ }
+}
+
+/* Restores the correlation grid state from the correlation grid history. */
+void CorrelationGrid::restoreStateFromHistory(const CorrelationGridHistory &correlationGridHistory)
+{
+ if (tensors_.size() != static_cast<size_t>(correlationGridHistory.numCorrelationTensors))
+ {
+ GMX_THROW(InvalidInputError("Mismatch of the grid size for the force correlation between checkpoint and simulation. Likely you have provided a checkpoint from a different simulation."));
+ }
+
+ /* Extract the state from the linear history array */
+ size_t bufferIndex = 0;
+ for (CorrelationTensor &tensor : tensors_)
+ {
+ tensor.restoreFromHistory(correlationGridHistory.blockDataBuffer,
+ &bufferIndex);
+ }
+
+ if (bufferIndex != correlationGridHistory.blockDataBuffer.size())
+ {
+ GMX_THROW(InvalidInputError("Mismatch of the correlation tensor size for the force correlation between checkpoint and simulation. Likely you have provided a checkpoint from a different simulation."));
+ }
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015,2016,2017, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief
+ * Contains datatypes and function declarations needed by AWH to
+ * have its force correlation data checkpointed.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_awh
+ */
+
+#ifndef GMX_AWH_CORRELATIONHISTORY_H
+#define GMX_AWH_CORRELATIONHISTORY_H
+
+struct t_commrec;
+
+namespace gmx
+{
+class CorrelationGrid;
+struct CorrelationGridHistory;
+
+/*! \brief
+ * Allocate a correlation grid history with the same structure as the given correlation grid.
+ *
+ * This function would be called at the start of a new simulation.
+ * Note that only sizes and memory are initialized here.
+ * History data is set by \ref updateCorrelationGridHistory.
+ *
+ * \param[in,out] corrGrid Correlation grid state to initialize with.
+ * \returns the correlation grid history struct.
+ */
+CorrelationGridHistory initCorrelationGridHistoryFromState(const CorrelationGrid &corrGrid);
+
+/*! \brief
+ * Restores the correlation grid state from the correlation grid history.
+ *
+ * \param[in] corrGridHist Correlation grid history to read.
+ * \param[in,out] corrGrid Correlation grid state to set.
+ */
+void restoreCorrelationGridStateFromHistory(const CorrelationGridHistory &corrGridHist, CorrelationGrid *corrGrid);
+
+/*! \brief
+ * Update the correlation grid history for checkpointing.
+ *
+ * \param[in,out] corrGridHist Correlation grid history to set.
+ * \param[in] corrGrid Correlation grid state to read.
+ */
+void updateCorrelationGridHistory(CorrelationGridHistory *corrGridHist, const CorrelationGrid &corrGrid);
+
+} // namespace gmx
+
+#endif /* GMX_AWH_CORRELATIONHISTORY_H */
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015,2016,2017, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Implements the CorrelationTensor class for correlation tensor statistics.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_awh
+ */
+
+#include "gmxpre.h"
+
+#include "correlationtensor.h"
+
+#include <assert.h>
+
+#include "gromacs/math/functions.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/utility/gmxassert.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+/*! \brief
+ * Get the block index for the current block and simulation length.
+ *
+ * This is simply how many blocks of length \p blockLength fit completely
+ * into \p totalAccumulatedLength, which is either the current time minus
+ * the start time, which time weighting, or the total weight.
+ *
+ * \param[in] blockLength Block length.
+ * \param[in] currentAccumulatedLength Sampling length of all data collected up to the current step, in time or weight.
+ * \returns the block index.
+ */
+int getBlockIndex(double blockLength,
+ double currentAccumulatedLength)
+{
+ return static_cast<int>(currentAccumulatedLength/blockLength);
+}
+
+} // namespace
+
+double CorrelationTensor::getTimeIntegral(int tensorIndex,
+ double dtSample) const
+{
+ const CorrelationBlockData &blockData = blockDataList_[0];
+ double weight = blockData.sumOverBlocksBlockSquareWeight();
+ double correlationIntegral = 0;
+ if (weight > 0)
+ {
+ correlationIntegral = blockData.correlationIntegral()[tensorIndex]/weight;
+ }
+
+ return 0.5*correlationIntegral*dtSample;
+}
+
+double CorrelationTensor::getVolumeElement(double dtSample) const
+{
+ double det;
+
+ switch (blockDataList_[0].correlationIntegral().size())
+ {
+ case 1:
+ /* 1-dimensional tensor: [a] */
+ det = getTimeIntegral(0, dtSample);
+ break;
+ case 3:
+ {
+ /* 2-dimensional tensor: [a b; b c] */
+ double a = getTimeIntegral(0, dtSample);
+ double b = getTimeIntegral(1, dtSample);
+ double c = getTimeIntegral(2, dtSample);
+
+ det = a*c - b*b;
+ }
+ break;
+ case 6:
+ {
+ /* 3-dimensional tensor: [a b d; b c e; d e f] */
+ double a = getTimeIntegral(0, dtSample);
+ double b = getTimeIntegral(1, dtSample);
+ double c = getTimeIntegral(2, dtSample);
+ double d = getTimeIntegral(3, dtSample);
+ double e = getTimeIntegral(4, dtSample);
+ double f = getTimeIntegral(5, dtSample);
+
+ det = a*c*f + 2*b*d*e - d*c*d - b*b*f - a*e*e;
+ }
+ break;
+ default:
+ det = 0;
+ /* meh */
+ break;
+ }
+
+ /* Returns 0 if no data, not supported number of dims
+ or not enough data to give a positive determinant (as it should be) */
+ return det > 0 ? std::sqrt(det) : 0;
+}
+
+void CorrelationTensor::doubleBlockLengths()
+{
+ /* We need to shift the data so that a given blockdata gets the data for double the block length.
+ The data for the shortest block length is not needed anymore. */
+
+ for (size_t i = 0; i < blockDataList_.size() - 1; i++)
+ {
+ blockDataList_[i] = blockDataList_[i + 1];
+ }
+
+ /* The blockdata which has 1 block is the same as the old one but with double the block length */
+ blockDataList_.back().doubleBlockLength();
+}
+
+void CorrelationTensor::updateBlockLengths(double samplingLength)
+{
+ /* How many times do we need to double the longest block length to fit the data? */
+ double longestLength = blockDataList_.back().blockLength();
+ int numDoublings = 0;
+ while (samplingLength > longestLength)
+ {
+ numDoublings++;
+ longestLength *= 2;
+ }
+
+ while (numDoublings > 0)
+ {
+ doubleBlockLengths();
+ numDoublings--;
+ }
+}
+
+/* Adds a filled data block to correlation time integral. */
+void CorrelationBlockData::addBlockToCorrelationIntegral()
+{
+ const bool firstBlock = (sumOverBlocksSquareBlockWeight_ == 0);
+
+ if (!firstBlock)
+ {
+ const int numDim = coordData_.size();
+ int tensorIndex = 0;
+ for (int d1 = 0; d1 < numDim; d1++)
+ {
+ const CoordData &c1 = coordData_[d1];
+
+ for (int d2 = 0; d2 <= d1; d2++)
+ {
+ const CoordData &c2 = coordData_[d2];
+
+ /* Compute change in correlaion integral due to adding
+ * the block through computing the difference of the block
+ * average with the old average for one component (we use x)
+ * and with the new component (we use y).
+ */
+ /* Need the old average, before the data of this block was added */
+ GMX_ASSERT(sumOverBlocksSquareBlockWeight_, "Denominator should be > 0 (should be guaranteed by the conditional above)");
+ double oldAverageX = c1.sumOverBlocksBlockWeightBlockWeightX/sumOverBlocksSquareBlockWeight_;
+
+ double newSumWeightBlockWeightY = c2.sumOverBlocksBlockWeightBlockWeightX + blockSumWeight_*c2.blockSumWeightX;
+ double newSumSquareBlockWeight = sumOverBlocksSquareBlockWeight_ + gmx::square(blockSumWeight_);
+
+ GMX_ASSERT(newSumSquareBlockWeight > 0, "Denominator should be > 0");
+ double newAverageY = newSumWeightBlockWeightY/newSumSquareBlockWeight;
+
+ double diffBlockWithOldAverageX = c1.blockSumWeightX - oldAverageX*blockSumWeight_;
+ double diffBlockWithNewAverageY = c2.blockSumWeightX - newAverageY*blockSumWeight_;
+
+ /* Update the correlation integral using the changes in averages. */
+ correlationIntegral_[tensorIndex] += diffBlockWithOldAverageX*diffBlockWithNewAverageY;
+ tensorIndex++;
+ }
+ }
+ }
+
+ /* Add the weights of the block to the block sums and clear the weights */
+ sumOverBlocksSquareBlockWeight_ += gmx::square(blockSumWeight_);
+ sumOverBlocksBlockSquareWeight_ += blockSumSquareWeight_;
+ for (auto &c : coordData_)
+ {
+ c.sumOverBlocksBlockWeightBlockWeightX += blockSumWeight_*c.blockSumWeightX;
+ /* Reset */
+ c.blockSumWeightX = 0;
+ }
+ /* Reset */
+ blockSumWeight_ = 0;
+ blockSumSquareWeight_ = 0;
+}
+
+void CorrelationTensor::addData(double weight,
+ gmx::ArrayRef<const double> data,
+ bool blockLengthInWeight,
+ double t)
+{
+ /* We should avoid adding data with very small weight to avoid
+ * divergence close to 0/0. The total spread weight for each sample is 1,
+ * so 1e-6 is a completely negligible amount.
+ */
+ if (weight < 1e-6)
+ {
+ /* Nothing to add */
+ return;
+ }
+
+ /* The sampling length is measured either in the total (local) weight or the current time */
+ double samplingLength = blockLengthInWeight ? getWeight() + weight : t;
+
+ /* Make sure the blocks are long enough to fit all data */
+ updateBlockLengths(samplingLength);
+
+ /* Store the data for each block length considered. First update the longest block which has all data since it's
+ used for updating the correlation function for the other block lengths. */
+ for (size_t i = 0; i < blockDataList_.size() - 1; i++)
+ {
+ CorrelationBlockData &bd = blockDataList_[i];
+
+ /* Find current block index given the block length. */
+ int blockIndex = getBlockIndex(bd.blockLength(), samplingLength);
+
+ if (bd.previousBlockIndex() >= 0 &&
+ bd.previousBlockIndex() != blockIndex)
+ {
+ /* Changed block. Update correlation with old data before adding to new block. */
+ bd.addBlockToCorrelationIntegral();
+ }
+
+ /* Keep track of which block index data was last added to */
+ bd.setPreviousBlockIndex(blockIndex);
+
+ /* Store the data */
+ bd.addData(weight, data);
+ }
+
+ /* The last blockdata has only 1 block which contains all data so far.
+ * Add the data for the largest block length.
+ */
+ blockDataList_.back().addData(weight, data);
+}
+
+CorrelationTensor::CorrelationTensor(int numDim,
+ int numBlockData,
+ double blockLengthInit)
+{
+ unsigned int scaling = 1;
+
+ GMX_RELEASE_ASSERT(numBlockData < static_cast<int>(sizeof(scaling)*8), "numBlockData should we smaller than the number of bits in scaling");
+
+ for (int n = 0; n < numBlockData; n++)
+ {
+ blockDataList_.push_back(CorrelationBlockData(numDim, scaling*blockLengthInit));
+ scaling <<= 1; /* Double block length */
+ }
+}
+
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015,2016,2017, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief
+ * Declares the CorrelationTensor class for correlation tensor statistics.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_awh
+ */
+
+#ifndef GMX_AWH_CORRELATIONTENSOR_H
+#define GMX_AWH_CORRELATIONTENSOR_H
+
+#include <cstddef>
+
+#include <vector>
+
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/gmxassert.h"
+
+namespace gmx
+{
+
+struct CorrelationBlockDataHistory;
+
+/*! \internal \brief Correlation block averaging weight-only data.
+ */
+class CorrelationBlockData
+{
+ public:
+ /*! \internal \brief Correlation block averaging data.
+ */
+ struct CoordData
+ {
+ /*! \brief Constructor.
+ */
+ CoordData() :
+ blockSumWeightX(0),
+ sumOverBlocksBlockWeightBlockWeightX(0)
+ {
+ };
+
+ double blockSumWeightX; /**< Weighted sum of x for current block. */
+ double sumOverBlocksBlockWeightBlockWeightX; /**< Sum over all blocks in the simulation of block weight times sum_wx. */
+ };
+
+ /*! \brief Constructor.
+ *
+ * \param[in] numDim The dimensionality.
+ * \param[in] blockLengthInit The initial block length.
+ */
+ CorrelationBlockData(int numDim,
+ double blockLengthInit) :
+ blockSumWeight_(0),
+ blockSumSquareWeight_(0),
+ sumOverBlocksSquareBlockWeight_(0),
+ sumOverBlocksBlockSquareWeight_(0),
+ blockLength_(blockLengthInit),
+ previousBlockIndex_(-1),
+ coordData_(numDim),
+ correlationIntegral_(numDim*(numDim + 1)/2)
+ {
+ };
+
+ /*! \brief Restore the state from history.
+ *
+ * \param[in] blockHistory The block data history containing the weight sums.
+ * \param[in] coordData The coordinate data.
+ * \param[in] correlationIntegral The correlation integral for all tensor elements.
+ */
+ void restoreFromHistory(const CorrelationBlockDataHistory &blockHistory,
+ const std::vector<CoordData> &coordData,
+ const std::vector<double> &correlationIntegral);
+
+ /*! \brief Adds a weighted data vector to one point in the correlation grid.
+ *
+ * \param[in] weight The weight of the data.
+ * \param[in] data One data point for each grid dimension.
+ */
+ void addData(double weight,
+ gmx::ArrayRef<const double> data)
+ {
+ GMX_ASSERT(data.size() == coordData_.size(), "Size of data should match the size of coordData");
+
+ blockSumWeight_ += weight;
+ blockSumSquareWeight_ += weight*weight;
+
+ for (size_t d = 0; d < coordData_.size(); d++)
+ {
+ coordData_[d].blockSumWeightX += weight*data[d];
+ }
+ };
+
+ /*! \brief Adds a filled data block to correlation time integral.
+ */
+ void addBlockToCorrelationIntegral();
+
+ /*! \brief Returns the sum weights for current block. */
+ double blockSumWeight() const
+ {
+ return blockSumWeight_;
+ };
+
+ /*! \brief Returns the sum weights^2 for current block. */
+ double blockSumSquareWeight() const
+ {
+ return blockSumSquareWeight_;
+ };
+
+ /*! \brief Returns the sum over blocks of block weight^2. */
+ double sumOverBlocksSquareBlockWeight() const
+ {
+ return sumOverBlocksSquareBlockWeight_;
+ };
+
+ /*! \brief Returns the sum over blocks of weight^2. */
+ double sumOverBlocksBlockSquareWeight() const
+ {
+ return sumOverBlocksBlockSquareWeight_;
+ };
+
+ /*! \brief Returns the length of each block used for block averaging. */
+ double blockLength() const
+ {
+ return blockLength_;
+ };
+
+ /*! \brief Double the length of each block used for block averaging. */
+ void doubleBlockLength()
+ {
+ blockLength_ *= 2;
+ };
+
+ /*! \brief Return the last block index data was added to (needed only for block length in terms of time). */
+ int previousBlockIndex() const
+ {
+ return previousBlockIndex_;
+ }
+
+ /*! \brief Set the last block index data was added to.
+ *
+ * \param[in] blockIndex The block index.
+ */
+ void setPreviousBlockIndex(int blockIndex)
+ {
+ previousBlockIndex_ = blockIndex;
+ }
+
+ /*! \brief Return sums for each coordinate dimension. */
+ const std::vector<CoordData> &coordData() const
+ {
+ return coordData_;
+ };
+
+ /*! \brief Return the correlation integral tensor. */
+ const std::vector<double> &correlationIntegral() const
+ {
+ return correlationIntegral_;
+ };
+
+ private:
+ /* Weight sum data, indentical for all dimensions */
+ double blockSumWeight_; /**< Sum weights for current block. */
+ double blockSumSquareWeight_; /**< Sum weights^2 for current block. */
+ double sumOverBlocksSquareBlockWeight_; /**< Sum over all blocks in the simulation of block weight^2. */
+ double sumOverBlocksBlockSquareWeight_; /**< Sum over all blocks in the simulation of weight^2. */
+ double blockLength_; /**< The length of each block used for block averaging */
+ int previousBlockIndex_; /**< The last block index data was added to (needed only for block length in terms of time). */
+
+ /* Sums for each coordinate dimension. */
+ std::vector<CoordData> coordData_; /**< Array with sums for each coordinate dimension. */
+
+ /* Correlation tensor. */
+ std::vector<double> correlationIntegral_; /**< Array with the correlation elements corr(x, y) in the tensor, where x, y are vector components. */
+};
+
+/*! \internal
+ * \brief Correlation data for computing the correlation tensor of one grid point.
+ *
+ * The time integrated autocorrelation of the desired quantity is computed using
+ * block averages, which is a computationally efficient and low memory method.
+ * Most of the work here goes into computing the block averages for weights
+ * and the coordinate quantity. This is done for a number of blocks in
+ * the range of \p c_numCorrelationBlocks/2 + 1 to \p c_numCorrelationBlocks,
+ * depending on the current simulation length. As the simulation time
+ * progresses, the blocks get longer. This is implemented in an efficient
+ * way by keeping track of log2(\p c_numCorrelationBlocks) \p BlockData
+ * data blocks with block length increasing progressively by a factor of 2.
+ * Once \p c_numCorrelationBlocks are reached, all block lengths are doubled.
+ */
+class CorrelationTensor
+{
+ public:
+ /*! \brief 64 blocks is a good trade-off between signal and noise */
+ static constexpr int c_numCorrelationBlocks = 64;
+
+ /*! \brief Constructor.
+ *
+ * \param[in] numDim The dimensionality.
+ * \param[in] numBlockData The number of data block data structs.
+ * \param[in] blockLengthInit The initial block length.
+ */
+ CorrelationTensor(int numDim,
+ int numBlockData,
+ double blockLengthInit);
+
+ /*! \brief Get a const reference to the list of block data.
+ */
+ const std::vector<CorrelationBlockData> &blockDataList() const
+ {
+ return blockDataList_;
+ }
+
+ /*! \brief
+ * Get the total weight of the data in the correlation matrix.
+ *
+ * \returns the weight of the added data.
+ */
+ double getWeight() const
+ {
+ /* The last blockdata has only 1 block containing all data */
+ return blockDataList().back().blockSumWeight();
+ };
+
+ /*! \brief Restore a correlation element from history.
+ *
+ * \param[in] blockDataBuffer The linear correlation grid data history buffer.
+ * \param[in,out] bufferIndex The index in \p blockDataBuffer to start reading, is increased with the number of blocks read.
+ */
+ void restoreFromHistory(const std::vector<CorrelationBlockDataHistory> &blockDataBuffer,
+ size_t *bufferIndex);
+
+ private:
+ /*! \brief Updates the block length by doubling.
+ *
+ * The length of all blocks is doubled. This is achieved by removing
+ * the shortest block, moving all other blocks and duplicating
+ * the data of longest block to a nw block of double length (but
+ * currenly only half filled with data).
+ */
+ void doubleBlockLengths();
+
+ /*! \brief Updates the block length such that data fits.
+ *
+ * \param[in] samplingLength Sampling length of all data, in time or weight.
+ */
+ void updateBlockLengths(double samplingLength);
+
+ public:
+ /*! \brief Adds a weighted data vector to one point in the correlation grid.
+ *
+ * \note To avoid rounding noise, data with weight smaller than 1e-6
+ * is ignored.
+ *
+ * \param[in] weight The weight of the data.
+ * \param[in] data One data point for each grid dimension.
+ * \param[in] blockLengthInWeight If true, a block is measured in probability weight, otherwise in time.
+ * \param[in] t The simulation time.
+ */
+ void addData(double weight,
+ gmx::ArrayRef<const double> data,
+ bool blockLengthInWeight,
+ double t);
+
+ /*! \brief Returns an element of the time integrated correlation tensor at a given point in the grid.
+ *
+ * The units of the integral are time*(units of data)^2. This will be friction units time/length^2
+ * if the data unit is 1/length.
+ *
+ * The correlation index lists the elements of the upper-triangular
+ * correlation matrix row-wise, so e.g. in 3D:
+ * 0 (0,0), 1 (1,0), 2 (1,1), 3 (2,0), 4 (2,1), 5 (2,2).
+ * (TODO: this should ideally not have to be known by the caller.)
+ *
+ * \param[in] tensorIndex Index in the tensor.
+ * \param[in] dtSample The sampling interval length.
+ * \returns the integral.
+ */
+ double getTimeIntegral(int tensorIndex,
+ double dtSample) const;
+
+ /*! \brief
+ * Returns the volume element of the correlation metric.
+ *
+ * The matrix of the metric equals the time-integrated correlation matrix. The volume element of
+ * the metric therefore equals the square-root of the absolute value of its determinant according
+ * to the standard formula for a volume element in a metric space.
+ *
+ * Since the units of the metric matrix elements are time*(units of data)^2, the volume element has
+ * units of (sqrt(time)*(units of data))^(ndim of data).
+ *
+ * \param[in] dtSample The sampling interval length.
+ * \returns the volume element.
+ */
+ double getVolumeElement(double dtSample) const;
+
+ private:
+ std::vector<CorrelationBlockData> blockDataList_; /**< The block data for different, consecutively doubling block lengths. */
+};
+
+} // namespace gmx
+
+#endif /* GMX_AWH_CORRELATIONTENSOR_H */
#include <gmock/gmock.h>
#include <gtest/gtest.h>
+#include "gromacs/awh/correlationgrid.h"
#include "gromacs/awh/pointstate.h"
#include "gromacs/mdtypes/awh-params.h"
#include "gromacs/utility/stringutil.h"
{
coordMaxValue = std::max(coordMaxValue, std::abs(coord));
- awh_dvec coordValue = { coord, 0, 0, 0 };
- awh_dvec biasForce;
- double potential = 0;
- bias.calcForceAndUpdateBias(coordValue,
- biasForce, &potential, &potentialJump,
- nullptr, step, step, seed_, nullptr);
+ awh_dvec coordValue = { coord, 0, 0, 0 };
+ double potential = 0;
+ gmx::ArrayRef<const double> biasForce =
+ bias.calcForceAndUpdateBias(coordValue,
+ &potential, &potentialJump,
+ nullptr, step, step, seed_, nullptr);
force.push_back(biasForce[0]);
pot.push_back(potential);
double coord = midPoint + halfWidth*(0.5*std::sin(t) + 0.55*std::sin(1.5*t));
awh_dvec coordValue = { coord, 0, 0, 0 };
- awh_dvec biasForce;
double potential = 0;
double potentialJump = 0;
bias.calcForceAndUpdateBias(coordValue,
- biasForce, &potential, &potentialJump,
+ &potential, &potentialJump,
nullptr, step, step, params.awhParams.seed, nullptr);
inInitialStage = bias.state().inInitialStage();
#include "gromacs/math/vec.h"
#include "gromacs/math/vecdump.h"
#include "gromacs/math/vectypes.h"
+#include "gromacs/mdtypes/awh-correlation-history.h"
#include "gromacs/mdtypes/awh-history.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/df_history.h"
eawhhUPDATELIST,
eawhhLOGSCALEDSAMPLEWEIGHT,
eawhhNUMUPDATES,
+ eawhhFORCECORRELATIONGRID,
eawhhNR
};
"awh_updatelist",
"awh_logScaledSampleWeight",
"awh_numupdates"
+ "awh_forceCorrelationGrid"
};
//! Higher level vector element type, only used for formatting checkpoint dumps
return 0;
}
+static int do_cpt_correlation_grid(XDR *xd, gmx_bool bRead, gmx_unused int fflags,
+ gmx::CorrelationGridHistory *corrGrid,
+ FILE *list, int eawhh)
+{
+ int ret = 0;
+
+ do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->numCorrelationTensors), list);
+ do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->tensorSize), list);
+ do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->blockDataListSize), list);
+
+ if (bRead)
+ {
+ initCorrelationGridHistory(corrGrid, corrGrid->numCorrelationTensors, corrGrid->tensorSize, corrGrid->blockDataListSize);
+ }
+
+ for (gmx::CorrelationBlockDataHistory &blockData : corrGrid->blockDataBuffer)
+ {
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeight), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumSquareWeight), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeightX), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeightY), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksSquareBlockWeight), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockSquareWeight), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockWeightBlockWeightX), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockWeightBlockWeightY), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockLength), list);
+ do_cpt_int_err(xd, eawhh_names[eawhh], &(blockData.previousBlockIndex), list);
+ do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.correlationIntegral), list);
+ }
+
+ return ret;
+}
+
static int do_cpt_awh_bias(XDR *xd, gmx_bool bRead,
int fflags, gmx::AwhBiasHistory *biasHistory,
FILE *list)
case eawhhNUMUPDATES:
do_cpt_step_err(xd, eawhh_names[i], &(state->numUpdates), list);
break;
+ case eawhhFORCECORRELATIONGRID:
+ ret = do_cpt_correlation_grid(xd, bRead, fflags, &biasHistory->forceCorrelationGrid, list, i);
+ break;
default:
gmx_fatal(FARGS, "Unknown awh history entry %d\n", i);
}
(1 << eawhhUMBRELLAGRIDPOINT) |
(1 << eawhhUPDATELIST) |
(1 << eawhhLOGSCALEDSAMPLEWEIGHT) |
- (1 << eawhhNUMUPDATES));
+ (1 << eawhhNUMUPDATES) |
+ (1 << eawhhFORCECORRELATIONGRID));
}
/* We can check many more things now (CPU, acceleration, etc), but
#include <string>
#include "gromacs/commandline/pargs.h"
+#include "gromacs/compat/make_unique.h"
#include "gromacs/fileio/enxio.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/oenv.h"
/*! \brief Constructor, Set the output base file name and title.
*
* Result is a valid object, but will produce empty output files.
+ *
+ * \param[in] filename The name for output files, frame time, and possibly bias number, will be added per file/frame.
+ * \param[in] baseTitle The base title of the plot, the bias number might be added.
+ * \param[in] numBias The total number of AWH biases in the system.
+ * \param[in] biasIndex The index of this bias.
*/
OutputFile(const std::string filename,
const std::string baseTitle,
int biasIndex);
/*! \brief Initializes the output file setup for the AWH output.
+ *
+ * \param[in] subBlockStart Index of the first sub-block to write in the energy frame.
+ * \param[in] numSubBlocks The total number of sub-blocks in the framw.
+ * \param[in] awhBiasParams The AWH bias parameters.
+ * \param[in] graphSelection Selects which graphs to plot.
+ * \param[in] energyUnit Requested energy unit in output.
+ * \param[in] kTValue kB*T in kJ/mol.
*/
- void initializeAwhOutputFile(int subblockStart,
- int numSubblocks,
+ void initializeAwhOutputFile(int subBlockStart,
+ int numSubBlocks,
const AwhBiasParams *awhBiasParams,
AwhGraphSelection graphSelection,
EnergyUnit energyUnit,
real kTValue);
+ /*! \brief Initializes the output file setup for the fricion output.
+ *
+ * \param[in] subBlockStart Index of the first sub-block to write in the energy frame.
+ * \param[in] numSubBlocks The total number of sub-blocks in the framw.
+ * \param[in] awhBiasParams The AWH bias parameters.
+ * \param[in] energyUnit Requested energy unit in output.
+ * \param[in] kTValue kB*T in kJ/mol.
+ */
+ void initializeFrictionOutputFile(int subBlockStart,
+ int numSubBlocks,
+ const AwhBiasParams *awhBiasParams,
+ EnergyUnit energyUnit,
+ real kTValue);
+
/*! \brief Opens a single output file for a bias, prints title and legends.
*
* \param[in] time The time for of frame to be written.
std::string baseFilename_; /**< Base of the output file name. */
std::string title_; /**< Title for the graph. */
int numDim_; /**< Number of dimensions. */
- int firstGraphSubblock_; /**< Index in the energy sub-blocks for the first graph. */
+ int firstGraphSubBlock_; /**< Index in the energy sub-blocks for the first graph. */
int numGraph_; /**< Number of actual graphs. */
bool useKTForEnergy_; /**< Whether to use kT instead of kJ/mol for energy. */
std::vector<real> scaleFactor_; /**< Scale factors from energy file data to output for each of the numGraph quantities. */
public:
//! Constructor.
BiasReader(int subblockStart,
- int numSubblocks,
- std::unique_ptr<OutputFile> awhOutputFile) :
+ int numSubBlocks,
+ std::unique_ptr<OutputFile> awhOutputFile,
+ std::unique_ptr<OutputFile> frictionOutputFile) :
subblockStart_(subblockStart),
- numSubblocks_(numSubblocks),
- awhOutputFile_(std::move(awhOutputFile))
+ numSubBlocks_(numSubBlocks),
+ awhOutputFile_(std::move(awhOutputFile)),
+ frictionOutputFile_(std::move(frictionOutputFile))
{
}
//! Return the AWH output file data.
const OutputFile &awhOutputFile() const
{
+ GMX_RELEASE_ASSERT(awhOutputFile_ != nullptr, "awhOutputFile() called without initialized AWH output file");
+
return *awhOutputFile_.get();
}
+ //! Return the a pointer to the friction output file data, can return nullptr
+ const OutputFile *frictionOutputFile() const
+ {
+ return frictionOutputFile_.get();
+ }
+
//! Return the starting subblock.
int subblockStart() const
{
}
private:
- const int subblockStart_; /**< The start index of the subblocks to read. */
- const int numSubblocks_; /**< Number of subblocks to read. */
- std::unique_ptr<OutputFile> awhOutputFile_; /**< The standard AWH output file data. */
- /* NOTE: A second OutputFile will be added soon, this will also make numSubblocks_ useful. */
+ const int subblockStart_; /**< The start index of the subblocks to read. */
+ const int numSubBlocks_; /**< Number of subblocks to read. */
+ std::unique_ptr<OutputFile> awhOutputFile_; /**< The standard AWH output file data. */
+ std::unique_ptr<OutputFile> frictionOutputFile_; /**< The friction/metric tensor output file data */
};
/*! \brief All options and meta-data needed for the AWH output */
namespace
{
-/* NOTE: A second value will be added soon. */
-enum {
- awhGraphTypeAwh
+//! Enum for selecting output file type.
+enum class OutputFileType
+{
+ Awh, //!< AWH, all data except friction tensor.
+ Friction //!< The full friction tensor.
};
/*! \brief The maximum number of observables per subblock, not including the full friction tensor.
* the output that mdrun writes. It would be better to define these
* values in a single location.
*/
-static constexpr int maxAwhGraphs = 5;
+static constexpr int maxAwhGraphs = 6;
/*! \brief Constructs a legend for a standard awh output file */
std::vector<std::string>makeLegend(const AwhBiasParams *awhBiasParams,
- int awhGraphType,
+ OutputFileType outputFileType,
size_t numLegend)
{
const std::array<std::string, maxAwhGraphs> legendBase =
"Coord bias",
"Coord distr",
"Ref value distr",
- "Target ref value distr" }
+ "Target ref value distr",
+ "Friction metric" }
};
std::vector<std::string> legend;
legend.push_back(gmx::formatString("dim%d", d));
}
- switch (awhGraphType)
+ switch (outputFileType)
{
- case awhGraphTypeAwh:
+ case OutputFileType::Awh:
{
/* Add as many legends as possible from the "base" legend list */
size_t legendBaseIndex = 0;
}
}
break;
+ case OutputFileType::Friction:
+ for (int i0 = 0; i0 < awhBiasParams->ndim; i0++)
+ {
+ for (int i1 = 0; i1 <= i0; i1++)
+ {
+ legend.push_back(gmx::formatString("%d,%d", i0, i1));
+ }
+ }
+ break;
}
GMX_RELEASE_ASSERT(legend.size() == numLegend, "The number of legends requested for printing and the number generated should be the same");
int numBias,
int biasIndex) :
numDim_(0),
- firstGraphSubblock_(0),
+ firstGraphSubBlock_(0),
numGraph_(0),
useKTForEnergy_(0),
scaleFactor_(),
}
void OutputFile::initializeAwhOutputFile(int subblockStart,
- int numSubblocks,
+ int numSubBlocks,
const AwhBiasParams *awhBiasParams,
AwhGraphSelection graphSelection,
EnergyUnit energyUnit,
{
/* The first subblock with actual graph y-values is index 1 + ndim */
numDim_ = awhBiasParams->ndim;
- firstGraphSubblock_ = subblockStart + 1 + numDim_;
+ firstGraphSubBlock_ = subblockStart + 1 + numDim_;
if (graphSelection == AwhGraphSelection::All)
{
/* There are one metadata and ndim coordinate blocks */
- numGraph_ = std::min(numSubblocks - 1 - numDim_,
+ numGraph_ = std::min(numSubBlocks - 1 - numDim_,
maxAwhGraphs);
}
else
std::fill(scaleFactor_.begin(), scaleFactor_.begin() + std::min(2, numGraph_), kTValue);
}
int numLegend = numDim_ - 1 + numGraph_;
- legend_ = makeLegend(awhBiasParams, awhGraphTypeAwh, numLegend);
+ legend_ = makeLegend(awhBiasParams, OutputFileType::Awh, numLegend);
/* We could have both length and angle coordinates in a single bias */
xLabel_ = "(nm or deg)";
yLabel_ = useKTForEnergy_ ? "(k\\sB\\NT)" : "(kJ/mol)";
}
}
+/*! \brief Initializes the output file setup for the fricion output (note that the filename is not set here). */
+void OutputFile::initializeFrictionOutputFile(int subBlockStart,
+ int numSubBlocks,
+ const AwhBiasParams *awhBiasParams,
+ EnergyUnit energyUnit,
+ real kTValue)
+{
+ /* The first subblock with actual graph y-values is index 1 + ndim */
+ numDim_ = awhBiasParams->ndim;
+ int numTensorElements = (numDim_*(numDim_ + 1))/2;
+
+ /* The friction tensor elements are always the last subblocks */
+ if (numSubBlocks < 1 + numDim_ + maxAwhGraphs + numTensorElements)
+ {
+ gmx_fatal(FARGS, "You requested friction tensor output, but the AWH data in the energy file does not contain the friction tensor");
+ }
+ GMX_ASSERT(numSubBlocks == 1 + numDim_ + maxAwhGraphs + numTensorElements, "The number of sub-blocks per bias should be 1 + ndim + maxAwhGraphs + (ndim*(ndim + 1))/2");
+
+ firstGraphSubBlock_ = subBlockStart + numSubBlocks - numTensorElements;
+ numGraph_ = numTensorElements;
+ useKTForEnergy_ = (energyUnit == EnergyUnit::KT);
+ scaleFactor_.resize(numGraph_, useKTForEnergy_ ? 1 : kTValue);
+ int numLegend = numDim_ - 1 + numGraph_;
+ legend_ = makeLegend(awhBiasParams, OutputFileType::Friction, numLegend);
+ xLabel_ = "(nm or deg)";
+ if (useKTForEnergy_)
+ {
+ yLabel_ = "friction/k\\sB\\NT (nm\\S-2\\Nps or rad\\S-2\\Nps)";
+ }
+ else
+ {
+ yLabel_ = "friction (kJ mol\\S-1\\Nnm\\S-2\\Nps or kJ mol\\S-1\\Nrad\\S-2\\Nps)";
+ }
+}
+
AwhReader::AwhReader(const AwhParams *awhParams,
int numFileOptions,
const t_filenm *filenames,
{
GMX_RELEASE_ASSERT(block != nullptr, "NULL pointer passed to initAwhReader");
+ bool outputFriction = opt2bSet("-fric", numFileOptions, filenames);
+
/* The first subblock of each AWH block has metadata about
* the number of subblocks belonging to that AWH block.
* (block->nsub gives only the total number of subblocks and not how
{
AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
- int numSubblocks = (int)block->sub[subblockStart].fval[0];
+ int numSubBlocks = (int)block->sub[subblockStart].fval[0];
- std::unique_ptr<OutputFile> outputFileAwh(new OutputFile(opt2fn("-o", numFileOptions, filenames), "AWH", awhParams->numBias, k));
+ std::unique_ptr<OutputFile> awhOutputFile(new OutputFile(opt2fn("-o", numFileOptions, filenames), "AWH", awhParams->numBias, k));
- outputFileAwh->initializeAwhOutputFile(subblockStart, numSubblocks,
+ awhOutputFile->initializeAwhOutputFile(subblockStart, numSubBlocks,
awhBiasParams, awhGraphSelection,
energyUnit, kT);
- biasReader_.emplace_back(BiasReader(subblockStart, numSubblocks, std::move(outputFileAwh)));
+ std::unique_ptr<OutputFile> frictionOutputFile;
+ if (outputFriction)
+ {
+ frictionOutputFile = gmx::compat::make_unique<OutputFile>(opt2fn("-fric", numFileOptions, filenames), "Friction tensor", awhParams->numBias, k);
+
+ frictionOutputFile->initializeFrictionOutputFile(subblockStart, numSubBlocks, awhBiasParams, energyUnit, kT);
+ }
+
+ biasReader_.emplace_back(BiasReader(subblockStart, numSubBlocks,
+ std::move(awhOutputFile),
+ std::move(frictionOutputFile)));
- subblockStart += numSubblocks;
+ subblockStart += numSubBlocks;
}
}
/* Print numGraph observables */
for (int i = 0; i < numGraph_; i++)
{
- fprintf(fp, " %g", block.sub[firstGraphSubblock_ + i].fval[j]*scaleFactor_[i]);
+ fprintf(fp, " %g", block.sub[firstGraphSubBlock_ + i].fval[j]*scaleFactor_[i]);
}
fprintf(fp, "\n");
for (auto &biasReader : biasReader_)
{
+ const int subStart = biasReader.subblockStart();
+
/* Each frame and AWH instance extracted generates one xvg file. */
{
const OutputFile &awhOutputFile = biasReader.awhOutputFile();
- const int subStart = biasReader.subblockStart();
-
FILE *fpAwh = awhOutputFile.openBiasOutputFile(time, oenv);
/* Now do the actual printing. Metadata in first subblock is treated separately. */
gmx_ffclose(fpAwh);
}
+
+ const OutputFile *frictionOutputFile = biasReader.frictionOutputFile();
+ if (frictionOutputFile != nullptr)
+ {
+ FILE *fpFriction = frictionOutputFile->openBiasOutputFile(time, oenv);
+
+ frictionOutputFile->writeData(block, subStart, fpFriction);
+
+ gmx_ffclose(fpFriction);
+ }
}
}
{
const char *desc[] = {
"[THISMODULE] extracts AWH data from an energy file.",
- "One file is written per AWH bias per time frame.",
+ "One or two files are written per AWH bias per time frame.",
"The bias index, if more than one, is appended to the file, as well as",
"the time of the frame. By default only the PMF is printed.",
"With [TT]-more[tt] the bias, target and coordinate distributions",
"are also printed.",
+ "With [TT]-more[tt] the bias, target and coordinate distributions",
+ "are also printed, as well as the metric sqrt(det(friction_tensor))",
+ "normalized such that the average is 1.",
+ "Option [TT]-fric[tt] prints all components of the friction tensor",
+ "to an additional set of files."
};
static gmx_bool moreGraphs = FALSE;
static int skip = 0;
t_filenm fnm[] = {
{ efEDR, "-f", nullptr, ffREAD },
{ efTPR, "-s", nullptr, ffREAD },
- { efXVG, "-o", "awh", ffWRITE }
+ { efXVG, "-o", "awh", ffWRITE },
+ { efXVG, "-fric", "friction", ffOPTWR }
};
const int nfile = asize(fnm);
if (!parse_common_args(&argc, argv,
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2015,2016,2017, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ *
+ * \brief
+ * Contains datatypes and function declarations needed by AWH to
+ * have its force correlation data checkpointed.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \inlibraryapi
+ * \ingroup module_mdtypes
+ */
+
+#ifndef GMX_MDTYPES_AWH_CORRELATION_HISTORY_H
+#define GMX_MDTYPES_AWH_CORRELATION_HISTORY_H
+
+#include <vector>
+
+namespace gmx
+{
+
+/*! \cond INTERNAL */
+
+//! Correlation block averaging data.
+struct CorrelationBlockDataHistory
+{
+ double blockSumWeight; /**< Sum weights for current block. */
+ double blockSumSquareWeight; /**< Sum weights^2 for current block. */
+ double blockSumWeightX; /**< Weighted sum of x for current block. */
+ double blockSumWeightY; /**< Weighted sum of y for current block. */
+ double sumOverBlocksSquareBlockWeight; /**< Sum over all blocks in the simulation of block weight^2 over the whole simulation. */
+ double sumOverBlocksBlockSquareWeight; /**< Sum over all blocks in the simulation of weight^2 over the whole simulation. */
+ double sumOverBlocksBlockWeightBlockWeightX; /**< Sum over all blocks in the simulation of block weight times blockSumWeightX over the whole simulation. */
+ double sumOverBlocksBlockWeightBlockWeightY; /**< Sum over all blocks in the simulation of block weight times blockSumWeightY over the whole simulation. */
+ double blockLength; /**< The length of each block used for block averaging. */
+ int previousBlockIndex; /**< The last block index data was added to (needed only for block length in terms of time). */
+ double correlationIntegral; /**< The time integral of the correlation function of x and y, corr(x(0), y(t)). */
+};
+
+//! Grid of local correlation matrices.
+struct CorrelationGridHistory
+{
+ /* These counts here since we curently need them for initializing the correlation grid when reading a checkpoint */
+ int numCorrelationTensors; /**< Number correlation tensors in the grid (equal to the number of points). */
+ int tensorSize; /**< The number of stored correlation matrix elements. */
+ int blockDataListSize; /**< To be able to increase the block length later on, data is saved for several block lengths for each element. */
+
+ /* We store all tensor sequentially in a buffer */
+ std::vector<CorrelationBlockDataHistory> blockDataBuffer; /**< Buffer that contains the correlation data. */
+};
+
+/*! \endcond */
+
+/*! \brief
+ * Initialize correlation grid history, sets all sizes.
+ *
+ * \param[in,out] correlationGridHistory Correlation grid history for master rank.
+ * \param[in] numCorrelationTensors Number of correlation tensors in the grid.
+ * \param[in] tensorSize Number of correlation elements in each tensor.
+ * \param[in] blockDataListSize The number of blocks in the list of each tensor element.
+ */
+void initCorrelationGridHistory(CorrelationGridHistory *correlationGridHistory,
+ int numCorrelationTensors,
+ int tensorSize,
+ int blockDataListSize);
+
+} // namespace gmx
+
+#endif /* GMX_MDTYPES_AWH_CORRELATION_HISTORY_H */
#include <vector>
+#include "gromacs/mdtypes/awh-correlation-history.h"
#include "gromacs/utility/basedefinitions.h"
namespace gmx
//! AWH bias history data. Note that this is a copy of an AWH internal struct.
struct AwhBiasHistory
{
- std::vector<AwhPointStateHistory> pointState; /**< History for grid coordinate points. */
+ std::vector<AwhPointStateHistory> pointState; /**< History for grid coordinate points. */
- AwhBiasStateHistory state; /**< The global state of the AWH bias. */
+ AwhBiasStateHistory state; /**< The global state of the AWH bias. */
+ CorrelationGridHistory forceCorrelationGrid; /**< History for force correlation statistics. */
/*! \brief Constructor. */
AwhBiasHistory() : pointState(),
- state()
+ state(),
+ forceCorrelationGrid()
{
}
};
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff