each input file should contain the coordinate values, such that each row defines a point in
coordinate space. Column :mdp:`awh1-ndim` + 1 should contain the PMF value for each point.
The target distribution column can either follow the PMF (column :mdp:`awh1-ndim` + 2) or
- be in the same column as written by ``gmx awh``.
+ be in the same column as written by :ref:`gmx awh`.
.. mdp:: awh1-share-group
#include <algorithm>
+#include "gromacs/fileio/enxio.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/math/units.h"
#include "gromacs/mdtypes/awh-history.h"
BiasCoupledToSystem::BiasCoupledToSystem(Bias bias,
const std::vector<int> &pullCoordIndex) :
- bias(bias),
+ bias(std::move(bias)),
pullCoordIndex(pullCoordIndex)
{
/* We already checked for this in grompp, but check again here. */
const std::string &biasInitFilename,
pull_t *pull_work) :
seed_(awhParams.seed),
+ nstout_(awhParams.nstOut),
commRecord_(commRecord),
pull_(pull_work),
potentialOffset_(0)
Awh::~Awh() = default;
+bool Awh::isOutputStep(gmx_int64_t step) const
+{
+ return (nstout_ > 0 && step % nstout_ == 0);
+}
+
real Awh::applyBiasForcesAndUpdateBias(int ePBC,
const t_mdatoms &mdatoms,
const matrix box,
biasForce[d], &mdatoms,
forceWithVirial);
}
+
+ if (isOutputStep(step))
+ {
+ /* Ensure that we output fully updated data */
+ biasCts.bias.doSkippedUpdatesForAllPoints();
+ }
}
wallcycle_stop(wallcycle, ewcAWH);
}
}
+/* Fill the AWH data block of an energy frame with data (if there is any). */
+void Awh::writeToEnergyFrame(gmx_int64_t step,
+ t_enxframe *frame) const
+{
+ GMX_ASSERT(MASTER(commRecord_), "writeToEnergyFrame should only be called on the master rank");
+ GMX_ASSERT(frame != nullptr, "Need a valid energy frame");
+
+ if (!isOutputStep(step))
+ {
+ /* This is not an AWH output step, don't write any AWH data */
+ return;
+ }
+
+ /* Get the total number of energy subblocks that AWH needs */
+ int numSubblocks = 0;
+ for (auto &biasCoupledToSystem : biasCoupledToSystem_)
+ {
+ numSubblocks += biasCoupledToSystem.bias.numEnergySubblocksToWrite();
+ }
+ GMX_ASSERT(numSubblocks > 0, "We should always have data to write");
+
+ /* Add 1 energy block */
+ add_blocks_enxframe(frame, frame->nblock + 1);
+
+ /* Take the block that was just added and set the number of subblocks. */
+ t_enxblock *awhEnergyBlock = &(frame->block[frame->nblock - 1]);
+ add_subblocks_enxblock(awhEnergyBlock, numSubblocks);
+
+ /* Claim it as an AWH block. */
+ awhEnergyBlock->id = enxAWH;
+
+ /* Transfer AWH data blocks to energy sub blocks */
+ int energySubblockCount = 0;
+ for (auto &biasCoupledToSystem : biasCoupledToSystem_)
+ {
+ energySubblockCount += biasCoupledToSystem.bias.writeToEnergySubblocks(&(awhEnergyBlock->sub[energySubblockCount]));
+ }
+}
+
} // namespace gmx
struct pull_t;
class t_state;
struct t_commrec;
+struct t_enxframe;
struct t_inputrec;
struct t_mdatoms;
* \param[in] box Box vectors.
* \param[in,out] forceWithVirial Force and virial buffers, should cover at least the local atoms.
* \param[in] t Time.
- * \param[in] step Time step.
+ * \param[in] step The current MD step.
* \param[in,out] wallcycle Wallcycle counter, can be nullptr.
* \param[in,out] fplog General output file, normally md.log, can be nullptr.
* \returns the potential energy for the bias.
*/
void restoreStateFromHistory(const AwhHistory *awhHistory);
+ /*! \brief Fills the AWH data block of an energy frame with data at certain steps.
+ *
+ * \param[in] step The current MD step.
+ * \param[in,out] fr Energy data frame.
+ */
+ void writeToEnergyFrame(gmx_int64_t step,
+ t_enxframe *fr) const;
+
/*! \brief Returns string "AWH" for registering AWH as an external potential provider with the pull module.
*/
static const char *externalPotentialString();
static void registerAwhWithPull(const AwhParams &awhParams,
pull_t *pull_work);
+ private:
+ /*! \brief Returns whether we need to write output at the current step.
+ *
+ * \param[in] step The current MD step.
+ */
+ bool isOutputStep(gmx_int64_t step) const;
+
private:
std::vector<BiasCoupledToSystem> biasCoupledToSystem_; /**< AWH biases and definitions of their coupling to the system. */
const gmx_int64_t seed_; /**< Random seed for MC jumping with umbrella type bias potential. */
+ const int nstout_; /**< Interval in steps for writing to energy file. */
const t_commrec *commRecord_; /**< Pointer to the communication record. */
pull_t *pull_; /**< Pointer to the pull working data. */
double potentialOffset_; /**< The offset of the bias potential which changes due to bias updates. */
void Bias::doSkippedUpdatesForAllPoints()
{
- state_.doSkippedUpdatesForAllPoints(params_);
+ if (params_.skipUpdates())
+ {
+ state_.doSkippedUpdatesForAllPoints(params_);
+ }
}
void Bias::calcForceAndUpdateBias(const awh_dvec coordValue,
updateList_.reserve(grid_.numPoints());
state_.initGridPointState(awhBiasParams, dimParams_, grid_, params_, biasInitFilename, awhParams.numBias);
+
+ if (thisRankDoesIO_)
+ {
+ writer_ = std::unique_ptr<BiasWriter>(new BiasWriter(*this));
+ }
+}
+
+/* Return the number of data blocks that have been prepared for writing. */
+int Bias::numEnergySubblocksToWrite() const
+{
+ GMX_RELEASE_ASSERT(writer_ != nullptr, "Should only request data from an initialized writer");
+
+ return writer_->numBlocks();
+}
+
+/* Write bias data blocks to energy subblocks. */
+int Bias::writeToEnergySubblocks(t_enxsubblock *subblock) const
+{
+ GMX_RELEASE_ASSERT(writer_ != nullptr, "Should only request data from an initialized writer");
+
+ return writer_->writeToEnergySubblocks(*this, subblock);
}
} // namespace gmx
#include "biasparams.h"
#include "biasstate.h"
+#include "biaswriter.h"
#include "dimparams.h"
#include "grid.h"
struct gmx_multisim_t;
struct t_commrec;
+struct t_enxsubblock;
namespace gmx
{
* sense as a relative, to other coordinate values, measure of the bias.
*
* \param[in] coordValue The coordinate value.
- * \returns the convolved bias >= -GMX_DOUBLE_MAX.
+ * \returns the convolved bias >= -GMX_FLOAT_MAX.
*/
double calcConvolvedBias(const awh_dvec &coordValue) const
{
return params_.biasIndex;
}
+ /*! \brief Return the coordinate value for a grid point.
+ *
+ * \param[in] gridPointIndex The index of the grid point.
+ */
+ inline const awh_dvec &getGridCoordValue(size_t gridPointIndex) const
+ {
+ GMX_ASSERT(gridPointIndex < grid_.numPoints(), "gridPointIndex should be in the range of the grid");
+
+ return grid_.point(gridPointIndex).coordValue;
+ }
+
private:
/*! \brief
* Performs statistical checks on the collected histograms and warns if issues are detected.
gmx_int64_t step,
FILE *fplog);
+ public:
+ /*! \brief Return the number of data blocks that have been prepared for writing.
+ */
+ int numEnergySubblocksToWrite() const;
+
+ /*! \brief Write bias data blocks to energy subblocks.
+ *
+ * \param[in,out] subblock Energy subblocks to write to.
+ * \returns the number of subblocks written.
+ */
+ int writeToEnergySubblocks(t_enxsubblock *subblock) const;
+
/* Data members. */
private:
const std::vector<DimParams> dimParams_; /**< Parameters for each dimension. */
const bool thisRankDoesIO_; /**< Tells whether this MPI rank will do I/O (checkpointing, AWH output) */
+ /* 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.
*/
numSharedUpdate(getNumSharedUpdate(awhBiasParams, numSharingSimulations)),
updateWeight(numSamplesUpdateFreeEnergy_*numSharedUpdate),
localWeightScaling(eTarget == eawhtargetLOCALBOLTZMANN ? temperatureScaleFactor : 1),
+ initialErrorInKT(beta*awhBiasParams.errorInitial),
initialHistogramSize(getInitialHistogramSizeEstimate(dimParams, awhBiasParams, gridAxis, beta, numStepsSampleCoord_*mdTimeStep)),
convolveForce(awhParams.ePotential == eawhpotentialCONVOLVED),
biasIndex(biasIndex),
const int numSharedUpdate; /**< The number of (multi-)simulations sharing the bias update */
const double updateWeight; /**< The probability weight accumulated for each update. */
const double localWeightScaling; /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */
+ const double initialErrorInKT; /**< Estimated initial free energy error in kT. */
const double initialHistogramSize; /**< Initial reference weight histogram size. */
private:
awh_ivec coverRadius_; /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */
namespace gmx
{
-void BiasState::getPmf(std::vector<float> *pmf) const
+void BiasState::getPmf(gmx::ArrayRef<float> pmf) const
{
+ GMX_ASSERT(pmf.size() == points_.size(), "pmf should have the size of the bias grid");
+
/* The PMF is just the negative of the log of the sampled PMF histogram.
* Points with zero target weight are ignored, they will mostly contain noise.
*/
-
- pmf->resize(points_.size());
-
for (size_t i = 0; i < points_.size(); i++)
{
- (*pmf)[i] = points_[i].inTargetRegion() ? -points_[i].logPmfSum() : GMX_FLOAT_MAX;
+ pmf[i] = points_[i].inTargetRegion() ? -points_[i].logPmfSum() : GMX_FLOAT_MAX;
}
}
*/
double freeEnergyMinimumValue(gmx::ArrayRef<const PointState> pointState)
{
- double fMin = GMX_DOUBLE_MAX;
+ double fMin = GMX_FLOAT_MAX;
for (auto const &ps : pointState)
{
convolvedPmf->resize(numPoints);
/* Get the PMF to convolve. */
- std::vector<float> pmf;
- getPmf(&pmf);
+ std::vector<float> pmf(numPoints);
+ getPmf(pmf);
for (size_t m = 0; m < numPoints; m++)
{
FILE *fplog,
int maxNumWarnings) const
{
- const double maxHistogramRatio = 0.5; /* Tolerance for printing a warning about the histogram ratios */
+ const double maxHistogramRatio = 0.5; /* Tolerance for printing a warning about the histogram ratios */
/* Sum up the histograms and get their normalization */
double sumVisits = 0;
for (size_t d = 0; d < dimParams.size(); d++)
{
/* clang thinks newForce[d] can be garbage */
-#ifndef __clang_analyzer__
+ #ifndef __clang_analyzer__
/* Average of the current and new force */
biasForce[d] = 0.5*(biasForce[d] + newForce[d]);
-#endif
+ #endif
}
return newPotential;
of the target region are ignored. */
/* Set the free energy cutoff */
- double maxFreeEnergy = GMX_DOUBLE_MAX;
+ double maxFreeEnergy = GMX_FLOAT_MAX;
if (params.eTarget == eawhtargetCUTOFF)
{
weightSum += std::exp(logWeight);
}
- /* Returns -GMX_DOUBLE_MAX if no neighboring points were in the target region. */
- return (weightSum > 0) ? std::log(weightSum) : -GMX_DOUBLE_MAX;
+ /* Returns -GMX_FLOAT_MAX if no neighboring points were in the target region. */
+ return (weightSum > 0) ? std::log(weightSum) : -GMX_FLOAT_MAX;
}
void BiasState::sampleProbabilityWeights(const Grid &grid,
* \param[in] dimParams The bias dimensions parameters
* \param[in] grid The grid.
* \param[in] coordValue Coordinate value.
- * \returns the convolved bias >= -GMX_DOUBLE_MAX.
+ * \returns the convolved bias >= -GMX_FLOAT_MAX.
*/
double calcConvolvedBias(const std::vector<DimParams> &dimParams,
const Grid &grid,
const awh_dvec &coordValue) const;
/*! \brief
- * Fills the given array with PMF values, resizes if necessary.
+ * Fills the given array with PMF values.
*
* Points outside of the biasing target region will get PMF = GMX_FLOAT_MAX.
* \note: The PMF is in single precision, because it is a statistical
* quantity and therefore never reaches full float precision.
*
- * \param[in,out] pmf Array returned will be of the same length as the AWH grid to store the PMF in.
+ * \param[out] pmf Array(ref) to be filled with the PMF values, should have the same size as the bias grid.
*/
- void getPmf(std::vector<float> *pmf) const;
+ void getPmf(gmx::ArrayRef<float>) const;
/*! \brief Returns the current coordinate state.
*/
return histogramSize_.inInitialStage();
};
+ /*! \brief Returns the current histogram size.
+ */
+ inline HistogramSize histogramSize() const
+ {
+ return histogramSize_;
+ };
+
/* Data members */
private:
CoordState coordState_; /**< The Current coordinate state */
--- /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.
+ */
+
+#include "gmxpre.h"
+
+#include "biaswriter.h"
+
+#include <assert.h>
+
+#include <cmath>
+
+#include "gromacs/awh/awh.h"
+#include "gromacs/fileio/enxio.h"
+#include "gromacs/mdtypes/awh-params.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+
+#include "bias.h"
+#include "grid.h"
+#include "pointstate.h"
+
+namespace gmx
+{
+
+namespace
+{
+
+/*! \brief
+ * Map the output entry type to a normalization type.
+ *
+ * The data is written to energy file blocks in the order given by
+ * the iterator of this map, which is based on the enum value
+ * (and matches the order of the lines below).
+ */
+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 }
+};
+
+/*! \brief
+ * Gets the coordinate normalization value for the given dimension.
+ *
+ * \param[in] bias The AWH bias.
+ * \param[in] dimIndex Dimensional index.
+ * \returns the coordinate normalization value.
+ */
+float getCoordNormalizationValue(const Bias &bias,
+ int dimIndex)
+{
+ /* AWH may use different units internally but here we convert to user units */
+ return bias.dimParams()[dimIndex].scaleInternalToUserInput(1);
+}
+
+/*! \brief
+ * Gets the normalization value for the given output entry type.
+ *
+ * \param[in] outputType Output entry type.
+ * \param[in] bias The AWH bias.
+ * \param[in] numBlocks The number of blocks for this output type.
+ * \returns the normalization value.
+ */
+float getNormalizationValue(AwhOutputEntryType outputType,
+ const Bias &bias,
+ int numBlocks)
+{
+ float normalizationValue = 0;
+
+ switch (outputType)
+ {
+ case AwhOutputEntryType::CoordValue:
+ normalizationValue = getCoordNormalizationValue(bias, numBlocks);
+ break;
+ case AwhOutputEntryType::Visits:
+ case AwhOutputEntryType::Weights:
+ case AwhOutputEntryType::Target:
+ normalizationValue = static_cast<float>(bias.state().points().size());
+ break;
+ default:
+ break;
+ }
+
+ return normalizationValue;
+}
+
+} // namespace
+
+AwhEnergyBlock::AwhEnergyBlock(int numPoints,
+ Normalization normalizationType,
+ float normalizationValue) :
+ normalizationType(normalizationType),
+ normalizationValue(normalizationValue),
+ data_(numPoints)
+{
+}
+
+BiasWriter::BiasWriter(const Bias &bias)
+{
+ std::map<AwhOutputEntryType, int> outputTypeNumBlock; /* Number of blocks per output type */
+
+ /* Different output variable types need different number of blocks.
+ * We keep track of the starting block for each variable.
+ */
+ int blockCount = 0;
+ for (const auto &pair : outputTypeToNormalization)
+ {
+ const AwhOutputEntryType outputType = pair.first;
+ {
+ outputTypeToBlock_[outputType] = blockCount;
+
+ if (outputType == AwhOutputEntryType::CoordValue)
+ {
+ outputTypeNumBlock[outputType] = bias.ndim();
+ }
+ else
+ {
+ /* Most output variable types need one block */
+ outputTypeNumBlock[outputType] = 1;
+ }
+ }
+ blockCount += outputTypeNumBlock[outputType];
+ }
+
+ /* Initialize the data blocks for each variable */
+ for (const auto &pair : outputTypeToNormalization)
+ {
+ const AwhOutputEntryType outputType = pair.first;
+ int numPoints;
+ if (outputType == AwhOutputEntryType::MetaData)
+ {
+ numPoints = static_cast<int>(AwhOutputMetaData::Count);
+ }
+ else
+ {
+ numPoints = bias.state().points().size();
+ }
+ for (int b = 0; b < outputTypeNumBlock[outputType]; b++)
+ {
+ block_.push_back(AwhEnergyBlock(numPoints,
+ pair.second,
+ getNormalizationValue(outputType, bias, b)));
+ }
+ }
+}
+
+/*! \brief
+ * Normalizes block data for output.
+ *
+ * \param[in,out] block The block to normalize.
+ * \param[in] bias The AWH bias.
+ */
+static void normalizeBlock(AwhEnergyBlock *block, const Bias &bias)
+{
+ gmx::ArrayRef<float> data = block->data();
+
+ /* Here we operate on float data (which is accurate enough, since it
+ * is statistical data that will never reach full float precision).
+ * But since we can have very many data points, we sum into a double.
+ */
+ double sum = 0;
+ float minValue = GMX_FLOAT_MAX;
+ float recipNorm = 0;
+
+ switch (block->normalizationType)
+ {
+ case Normalization::None:
+ break;
+ case Normalization::Coordinate:
+ /* Normalize coordinate values by a scale factor */
+ for (float &point : data)
+ {
+ point *= block->normalizationValue;
+ }
+ break;
+ case Normalization::FreeEnergy:
+ /* Normalize free energy values by subtracting the minimum value */
+ for (size_t index = 0; index < data.size(); index++)
+ {
+ if (bias.state().points()[index].inTargetRegion() && data[index] < minValue)
+ {
+ minValue = data[index];
+ }
+ }
+ for (size_t index = 0; index < data.size(); index++)
+ {
+ if (bias.state().points()[index].inTargetRegion())
+ {
+ data[index] -= minValue;
+ }
+ }
+
+ break;
+ case Normalization::Distribution:
+ /* Normalize distribution values by normalizing their sum */
+ for (float &point : data)
+ {
+ sum += point;
+ }
+ if (sum > 0)
+ {
+ recipNorm = block->normalizationValue/static_cast<float>(sum);
+ }
+ for (float &point : data)
+ {
+ point *= recipNorm;
+ }
+ break;
+ default:
+ GMX_RELEASE_ASSERT(false, "Unknown AWH normalization type");
+ break;
+ }
+}
+
+void BiasWriter::transferMetaDataToWriter(size_t metaDataIndex,
+ AwhOutputMetaData metaDataType,
+ const Bias &bias)
+{
+ gmx::ArrayRef<float> data = block_[getVarStartBlock(AwhOutputEntryType::MetaData)].data();
+ GMX_ASSERT(metaDataIndex < data.size(), "Attempt to transfer AWH meta data to block for index out of range");
+
+ /* Transfer the point data of this variable to the right block(s) */
+ switch (metaDataType)
+ {
+ case AwhOutputMetaData::NumBlock:
+ /* The number of subblocks per awh (needed by gmx_energy) */
+ data[metaDataIndex] = static_cast<double>(block_.size());
+ /* Note: a single subblock takes only a single type and we need doubles. */
+ break;
+ case AwhOutputMetaData::TargetError:
+ /* The theoretical target error */
+ data[metaDataIndex] = bias.params().initialErrorInKT*std::sqrt(bias.params().initialHistogramSize/bias.state().histogramSize().histogramSize());
+ break;
+ case AwhOutputMetaData::ScaledSampleWeight:
+ /* The logarithm of the sample weight relative to a sample weight of 1 at the initial time.
+ In the normal case: this will increase in the initial stage and then stay at a constant value. */
+ data[metaDataIndex] = bias.state().histogramSize().logScaledSampleWeight();
+ break;
+ case AwhOutputMetaData::Count:
+ break;
+ }
+}
+
+void
+BiasWriter::transferPointDataToWriter(AwhOutputEntryType outputType,
+ int pointIndex,
+ const Bias &bias,
+ gmx::ArrayRef<const float> pmf)
+{
+ /* The starting block index of this output type.
+ * Note that some variables need several (contiguous) blocks.
+ */
+ 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");
+
+ /* Transfer the point data of this variable to the right block(s) */
+ int b = blockStart;
+ switch (outputType)
+ {
+ case AwhOutputEntryType::MetaData:
+ GMX_RELEASE_ASSERT(false, "MetaData is handled by a different function");
+ break;
+ case AwhOutputEntryType::CoordValue:
+ {
+ const awh_dvec &coordValue = bias.getGridCoordValue(pointIndex);
+ for (int d = 0; d < bias.ndim(); d++)
+ {
+ block_[b].data()[pointIndex] = coordValue[d];
+ b++;
+ }
+ }
+ break;
+ case AwhOutputEntryType::Pmf:
+ block_[b].data()[pointIndex] = bias.state().points()[pointIndex].inTargetRegion() ? pmf[pointIndex] : 0;
+ break;
+ case AwhOutputEntryType::Bias:
+ {
+ const awh_dvec &coordValue = bias.getGridCoordValue(pointIndex);
+ block_[b].data()[pointIndex] = bias.state().points()[pointIndex].inTargetRegion() ? bias.calcConvolvedBias(coordValue) : 0;
+ }
+ break;
+ case AwhOutputEntryType::Visits:
+ block_[b].data()[pointIndex] = bias.state().points()[pointIndex].numVisitsTot();
+ break;
+ case AwhOutputEntryType::Weights:
+ block_[b].data()[pointIndex] = bias.state().points()[pointIndex].weightSumTot();
+ break;
+ case AwhOutputEntryType::Target:
+ block_[b].data()[pointIndex] = bias.state().points()[pointIndex].target();
+ break;
+ default:
+ GMX_RELEASE_ASSERT(false, "Unknown AWH output variable");
+ break;
+ }
+}
+
+void BiasWriter::prepareBiasOutput(const Bias &bias)
+{
+ /* Pack the AWH data into the writer data. */
+
+ /* Evaluate the PMF for all points */
+ gmx::ArrayRef<float> pmf = block_[getVarStartBlock(AwhOutputEntryType::Pmf)].data();
+ bias.state().getPmf(pmf);
+
+ /* Pack the the data point by point.
+ * Unfortunately we can not loop over a class enum, so we cast to int.
+ * \todo Use strings instead of enum when we port the output to TNG.
+ */
+ for (int i = 0; i < static_cast<int>(AwhOutputMetaData::Count); i++)
+ {
+ transferMetaDataToWriter(i, static_cast<AwhOutputMetaData>(i), bias);
+ }
+ for (const auto &pair : outputTypeToNormalization)
+ {
+ const AwhOutputEntryType outputType = pair.first;
+ /* Skip metadata (transfered above) and unused blocks */
+ if (outputType == AwhOutputEntryType::MetaData || !hasVarBlock(outputType))
+ {
+ continue;
+ }
+ for (size_t m = 0; m < bias.state().points().size(); m++)
+ {
+ transferPointDataToWriter(outputType, m, bias, pmf);
+ }
+ }
+
+ /* For looks of the output, normalize it */
+ for (AwhEnergyBlock &block : block_)
+ {
+ normalizeBlock(&block, bias);
+ }
+}
+
+int BiasWriter::writeToEnergySubblocks(const Bias &bias,
+ t_enxsubblock *sub)
+{
+ prepareBiasOutput(bias);
+
+ for (size_t b = 0; b < block_.size(); b++)
+ {
+ sub[b].type = xdr_datatype_float;
+ sub[b].nr = block_[b].data().size();
+ sub[b].fval = block_[b].data().data();
+ }
+
+ return block_.size();
+}
+
+} // namepace 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
+ * This file contains the BiasWriter class that prepares and writes data of a Bias to an energy file.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess <hess@kth.se>
+ * \ingroup module_awh
+ */
+
+#ifndef GMX_AWH_BIASWRITER_H
+#define GMX_AWH_BIASWRITER_H
+
+#include <map>
+#include <vector>
+
+#include "gromacs/fileio/enxio.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/basedefinitions.h"
+
+struct gmx_multisim_t;
+struct t_enxsubblock;
+
+namespace gmx
+{
+class Bias;
+
+/* TODO: the post-simulations AWH reader and this AWH writer are totally
+ * disconnected although they read/write the same data. I'm not sure how
+ * to handle that or if it should be left as it is until the writing is done
+ * in a differen format (i.e. TNG) than the current energy file.
+ */
+
+//! 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.
+};
+
+//! Enum with the types of metadata to write.
+enum class AwhOutputMetaData
+{
+ NumBlock, //!< The number of blocks.
+ TargetError, //!< The target error.
+ ScaledSampleWeight, //!< The logarithm of the sample weight relative to a sample weight of 1 at the initial time.
+ Count //!< The number of enum values, not including Count.
+};
+
+//! Enum with different ways of normalizing the output.
+enum class Normalization
+{
+ None, //!< No normalization.
+ Coordinate, //!< Scale using the internal/user input coordinate scaling factor.
+ FreeEnergy, //!< Normalize free energy values by subtracting the minimum value.
+ Distribution //!< Normalize the distribution to 1.
+};
+
+/*! \internal \brief AWH output data block that can be written to an energy file block.
+ */
+class AwhEnergyBlock
+{
+ public:
+ /*! \brief Constructor
+ *
+ * \param[in] numPoints Number of points in block.
+ * \param[in] normalizationType Value for normalization type enum.
+ * \param[in] normalizationValue Normalization value.
+ */
+ AwhEnergyBlock(int numPoints,
+ Normalization normalizationType,
+ float normalizationValue);
+
+ /*! \brief Returns an ArrarRef to the data in the block.
+ */
+ gmx::ArrayRef<float> data()
+ {
+ return data_;
+ }
+
+ const Normalization normalizationType; /**< How to normalize the output data */
+ const float normalizationValue; /**< The normalization value */
+ private:
+ std::vector<float> data_; /**< The data, always float which is enough since this is statistical data */
+};
+
+/*! \internal \brief Class organizing the output data storing and writing of an AWH bias.
+ */
+class BiasWriter
+{
+ private:
+ /*! \brief Query if the writer has a block for the given variable.
+ *
+ * \param[in] outputType Output entry type.
+ */
+ bool hasVarBlock(AwhOutputEntryType outputType) const
+ {
+ return outputTypeToBlock_.find(outputType)->second >= 0;
+ }
+
+ /*! \brief* Find the first block containing the given variable.
+ *
+ * \param[in] outputType Output entry type.
+ * \returns the first block index for the variable, or -1 there is no block.
+ */
+ int getVarStartBlock(AwhOutputEntryType outputType) const
+ {
+ return outputTypeToBlock_.find(outputType)->second;
+ }
+
+ /*! \brief Transfer AWH point data to writer data blocks.
+ *
+ * \param[in] metaDataIndex Meta data block index.
+ * \param[in] metaDataType The type of meta data to write.
+ * \param[in] bias The AWH Bias.
+ */
+ void transferMetaDataToWriter(size_t metaDataIndex,
+ AwhOutputMetaData metaDataType,
+ const Bias &bias);
+
+ /*! \brief Transfer AWH point data to writer data blocks.
+ *
+ * \param[in] outputType Output entry type.
+ * \param[in] pointIndex The point index.
+ * \param[in] bias The AWH Bias.
+ * \param[in] pmf PMF values.
+ */
+ void transferPointDataToWriter(AwhOutputEntryType outputType,
+ int pointIndex,
+ 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.
+ *
+ * \param[in] bias The AWH Bias.
+ */
+ 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 */
+};
+
+} // namespace gmx
+
+#endif /* GMX_AWH_BIASWRITER_H */
/*! \brief Returns true if we are in the initial stage of the AWH method.
*/
- inline bool inInitialStage() const
+ bool inInitialStage() const
{
return inInitialStage_;
};
+ /*! \brief Returns The log of the current sample weight, scaled because of the histogram rescaling.
+ */
+ double logScaledSampleWeight() const
+ {
+ return logScaledSampleWeight_;
+ };
+
private:
gmx_int64_t numUpdates_; /**< The number of updates performed since the start of the simulation. */
opt, awhParams->nstOut);
warning_error(wi, buf);
}
+ /* This restriction can be removed by changing a flag of print_ebin() */
+ if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
+ {
+ char buf[STRLEN];
+ sprintf(buf, "%s (%d) should be a multiple of nstenergy (%d)",
+ opt, awhParams->nstOut, ir->nstenergy);
+ warning_error(wi, buf);
+ }
CTYPE("Coordinate sampling interval in number of steps");
sprintf(opt, "%s-nstsample", prefix);
"Distance restraints",
"Free energy data",
"BAR histogram",
- "Delta H raw data"
+ "Delta H raw data",
+ "AWH data"
};
enxDHHIST, /* BAR histogram */
enxDH, /* BAR raw delta H data */
+ enxAWH, /* AWH data */
+
enxNR /* Total number of extra blocks in the current code,
* note that the enxio code can read files written by
* future code which contain more blocks.
/* the subblocks that are contained in energy file blocks. Each of these
has a number of values of a single data type in a .edr file. */
-typedef struct
+struct t_enxsubblock
{
int nr; /* number of items in subblock */
xdr_datatype type; /* the block type */
int lval_alloc;
int cval_alloc;
int sval_alloc;
-} t_enxsubblock;
+};
/* the energy file blocks. Each block contains a number of sub-blocks
/* The frames that are read/written */
-typedef struct {
+struct t_enxframe {
double t; /* Timestamp of this frame */
gmx_int64_t step; /* MD step */
gmx_int64_t nsteps; /* The number of steps between frames */
int nblock; /* Number of following energy blocks */
t_enxblock *block; /* The blocks */
int nblock_alloc; /* The number of blocks allocated */
-} t_enxframe;
+};
/* file handle */
typedef struct ener_file *ener_file_t;
}
}
-void xvgr_legend(FILE *out, int nsets, const char** setname,
- const gmx_output_env_t *oenv)
+static bool stringIsEmpty(const std::string s)
+{
+ return s.empty();
+}
+
+static bool stringIsEmpty(const char *s)
+{
+ return (s == nullptr || s[0] == '\0');
+}
+
+template <typename T>
+static void xvgr_legend(FILE *out, int nsets, const T * setname,
+ const gmx_output_env_t *oenv)
{
int i;
char buf[STRLEN];
fprintf(out, "@ legend length %d\n", 2);
for (i = 0; (i < nsets); i++)
{
- if (setname[i])
+ if (!stringIsEmpty(setname[i]))
{
if (output_env_get_xvg_format(oenv) == exvgXMGR)
{
}
}
+void xvgrLegend(FILE *out,
+ const std::vector<std::string> &setNames,
+ const struct gmx_output_env_t *oenv)
+{
+ xvgr_legend(out, setNames.size(), setNames.data(), oenv);
+}
+void xvgr_legend(FILE *out, int nsets, const char** setnames,
+ const struct gmx_output_env_t *oenv)
+{
+ xvgr_legend<char const *>(out, nsets, setnames, oenv);
+}
+
void xvgr_new_dataset(FILE *out, int nr_first, int nsets,
const char **setname,
const gmx_output_env_t *oenv)
#include <stdio.h>
#include <string>
+#include <vector>
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/real.h"
const struct gmx_output_env_t *oenv);
/* Set the world in xvgr */
+void xvgrLegend(FILE *out,
+ const std::vector<std::string> &setNames,
+ const struct gmx_output_env_t *oenv);
+/* Make a legend box, and also modifies the view to make room for the legend */
+
void xvgr_legend(FILE *out, int nsets, const char** setnames,
const struct gmx_output_env_t *oenv);
/* Make a legend box, and also modifies the view to make room for the legend */
-
void xvgr_new_dataset(FILE *out,
int nr_first, int nsets, const char **setnames,
const struct gmx_output_env_t *oenv);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2008, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,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.
int
gmx_anaeig(int argc, char *argv[]);
+int
+gmx_awh(int argc, char *argv[]);
+
int
gmx_g_angle(int argc, char *argv[]);
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2013,2014,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
+ * Tool for extracting AWH (Accelerated Weight Histogram) data from energy files.
+ *
+ * \author Viveca Lindahl
+ * \author Berk Hess
+ */
+
+#include "gmxpre.h"
+
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
+#include <array>
+#include <memory>
+#include <string>
+
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/fileio/enxio.h"
+#include "gromacs/fileio/gmxfio.h"
+#include "gromacs/fileio/oenv.h"
+#include "gromacs/fileio/tpxio.h"
+#include "gromacs/fileio/trxio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/math/units.h"
+#include "gromacs/mdtypes/awh-params.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/arraysize.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/stringutil.h"
+
+using gmx::AwhParams;
+using gmx::AwhBiasParams;
+
+namespace
+{
+
+//! Enum for choosing AWH graph output
+enum class AwhGraphSelection
+{
+ Pmf, //!< Only the PMF
+ All //!< All possible AWH graphs
+};
+
+//! Energy unit options
+enum class EnergyUnit
+{
+ KJPerMol, //!< kJ/mol
+ KT //!< kT
+};
+
+} // namespace
+
+/*! \brief All meta-data that is shared for one output file type for one bias */
+class OutputFile
+{
+ public:
+ /*! \brief Constructor, Set the output base file name and title.
+ *
+ * Result is a valid object, but will produce empty output files.
+ */
+ OutputFile(const std::string filename,
+ const std::string baseTitle,
+ int numBias,
+ int biasIndex);
+
+ /*! \brief Initializes the output file setup for the AWH output.
+ */
+ void initializeAwhOutputFile(int subblockStart,
+ int numSubblocks,
+ const AwhBiasParams *awhBiasParams,
+ AwhGraphSelection graphSelection,
+ 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.
+ * \param[in] oenv The output environment.
+ * \returns the file pointer.
+ */
+ FILE * openBiasOutputFile(double time,
+ const gmx_output_env_t *oenv) const;
+
+ /*! \brief Writes data selected from \p block to file.
+ *
+ * \param[in] block The energy block with the data to write.
+ * \param[in] subBlockStart The sub-block start index.
+ * \param[in] fp The file pointer.
+ */
+ void writeData(const t_enxblock &block,
+ int subBlockStart,
+ FILE *fp) const;
+
+ private:
+ 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 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. */
+
+ std::vector<std::string> legend_; /**< Legends for the output */
+ std::string xLabel_; /**< Label for the x-axis. */
+ std::string yLabel_; /**< Label for the y-axis. */
+};
+
+/*! \brief All meta-data that is shared for all output files for one bias */
+class BiasReader
+{
+ public:
+ //! Constructor.
+ BiasReader(int subblockStart,
+ int numSubblocks,
+ std::unique_ptr<OutputFile> awhOutputFile) :
+ subblockStart_(subblockStart),
+ numSubblocks_(numSubblocks),
+ awhOutputFile_(std::move(awhOutputFile))
+ {
+ }
+
+ //! Return the AWH output file data.
+ const OutputFile &awhOutputFile() const
+ {
+ return *awhOutputFile_.get();
+ }
+
+ //! Return the starting subblock.
+ int subblockStart() const
+ {
+ return subblockStart_;
+ }
+
+ 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. */
+};
+
+/*! \brief All options and meta-data needed for the AWH output */
+class AwhReader
+{
+ public:
+ //! Constructor
+ AwhReader(const AwhParams *awhParams,
+ int numFileOptions,
+ const t_filenm *filenames,
+ AwhGraphSelection awhGraphSelection,
+ EnergyUnit energyUnit,
+ real kT,
+ const t_enxblock *block);
+
+ //! Extract and print AWH data for an AWH block of one energy frame
+ void processAwhFrame(const t_enxblock &block,
+ double time,
+ const gmx_output_env_t *oenv) const;
+
+ private:
+ std::vector<BiasReader> biasReader_; /**< The readers, one for each AWH bias. */
+ public:
+ const real kT_; /**< kB*T in kJ/mol. */
+};
+
+namespace
+{
+
+/* NOTE: A second value will be added soon. */
+enum {
+ awhGraphTypeAwh
+};
+
+/*! \brief The maximum number of observables per subblock, not including the full friction tensor.
+ *
+ * This number, as well as the legends in make_legend() should match
+ * the output that mdrun writes. It would be better to define these
+ * values in a single location.
+ */
+static constexpr int maxAwhGraphs = 5;
+
+/*! \brief Constructs a legend for a standard awh output file */
+std::vector<std::string>makeLegend(const AwhBiasParams *awhBiasParams,
+ int awhGraphType,
+ size_t numLegend)
+{
+ const std::array<std::string, maxAwhGraphs> legendBase =
+ {
+ { "PMF",
+ "Coord bias",
+ "Coord distr",
+ "Ref value distr",
+ "Target ref value distr" }
+ };
+
+ std::vector<std::string> legend;
+ /* Give legends to dimensions higher than the first */
+ for (int d = 1; d < awhBiasParams->ndim; d++)
+ {
+ legend.push_back(gmx::formatString("dim%d", d));
+ }
+
+ switch (awhGraphType)
+ {
+ case awhGraphTypeAwh:
+ {
+ /* Add as many legends as possible from the "base" legend list */
+ size_t legendBaseIndex = 0;
+ while (legend.size() < numLegend && legendBaseIndex < legendBase.size())
+ {
+ legend.push_back(legendBase[legendBaseIndex]);
+ legendBaseIndex++;
+ }
+ }
+ break;
+ }
+
+ GMX_RELEASE_ASSERT(legend.size() == numLegend, "The number of legends requested for printing and the number generated should be the same");
+
+ return legend;
+}
+
+} // namespace
+
+OutputFile::OutputFile(const std::string filename,
+ const std::string baseTitle,
+ int numBias,
+ int biasIndex) :
+ numDim_(0),
+ firstGraphSubblock_(0),
+ numGraph_(0),
+ useKTForEnergy_(0),
+ scaleFactor_(),
+ legend_(),
+ xLabel_(),
+ yLabel_()
+{
+ // cppcheck-suppress useInitializationList
+ baseFilename_ = filename.substr(0, filename.find('.'));
+ title_ = baseTitle;
+ if (numBias > 1)
+ {
+ baseFilename_ += gmx::formatString("%d", biasIndex + 1);
+ title_ += gmx::formatString(" %d", biasIndex + 1);
+ }
+}
+
+void OutputFile::initializeAwhOutputFile(int subblockStart,
+ int numSubblocks,
+ const AwhBiasParams *awhBiasParams,
+ AwhGraphSelection graphSelection,
+ EnergyUnit energyUnit,
+ real kTValue)
+{
+ /* The first subblock with actual graph y-values is index 1 + ndim */
+ numDim_ = awhBiasParams->ndim;
+ firstGraphSubblock_ = subblockStart + 1 + numDim_;
+ if (graphSelection == AwhGraphSelection::All)
+ {
+ /* There are one metadata and ndim coordinate blocks */
+ numGraph_ = std::min(numSubblocks - 1 - numDim_,
+ maxAwhGraphs);
+ }
+ else
+ {
+ numGraph_ = 1;
+ }
+ useKTForEnergy_ = (energyUnit == EnergyUnit::KT);
+ scaleFactor_.resize(numGraph_, 1);
+ if (!useKTForEnergy_)
+ {
+ /* The first two graphs are in units of energy, multiply by kT */
+ std::fill(scaleFactor_.begin(), scaleFactor_.begin() + std::min(2, numGraph_), kTValue);
+ }
+ int numLegend = numDim_ - 1 + numGraph_;
+ legend_ = makeLegend(awhBiasParams, awhGraphTypeAwh, numLegend);
+ /* We could have both length and angle coordinates in a single bias */
+ xLabel_ = "(nm or deg)";
+ yLabel_ = useKTForEnergy_ ? "(k\\sB\\NT)" : "(kJ/mol)";
+ if (graphSelection == AwhGraphSelection::All)
+ {
+ yLabel_ += gmx::formatString(", (nm\\S-%d\\N or rad\\S-%d\\N), (-)",
+ awhBiasParams->ndim, awhBiasParams->ndim);
+ }
+}
+
+AwhReader::AwhReader(const AwhParams *awhParams,
+ int numFileOptions,
+ const t_filenm *filenames,
+ AwhGraphSelection awhGraphSelection,
+ EnergyUnit energyUnit,
+ real kT,
+ const t_enxblock *block) :
+ kT_(kT)
+{
+ GMX_RELEASE_ASSERT(block != nullptr, "NULL pointer passed to initAwhReader");
+
+ /* 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
+ * they are distributed between the AWH biases).
+ */
+
+ /* Keep track of the first subblock of this AWH */
+ int subblockStart = 0;
+ for (int k = 0; k < awhParams->numBias; k++)
+ {
+ AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
+
+ int numSubblocks = (int)block->sub[subblockStart].fval[0];
+
+ std::unique_ptr<OutputFile> outputFileAwh(new OutputFile(opt2fn("-o", numFileOptions, filenames), "AWH", awhParams->numBias, k));
+
+ outputFileAwh->initializeAwhOutputFile(subblockStart, numSubblocks,
+ awhBiasParams, awhGraphSelection,
+ energyUnit, kT);
+
+ biasReader_.emplace_back(BiasReader(subblockStart, numSubblocks, std::move(outputFileAwh)));
+
+ subblockStart += numSubblocks;
+ }
+}
+
+FILE * OutputFile::openBiasOutputFile(double time,
+ const gmx_output_env_t *oenv) const
+{
+ std::string filename = baseFilename_ + gmx::formatString("_t%g.xvg", time);
+
+ FILE *fp = xvgropen(filename.c_str(), title_.c_str(), xLabel_, yLabel_, oenv);
+ xvgrLegend(fp, legend_, oenv);
+
+ return fp;
+}
+
+/*! \brief Prints data selected by \p outputFile from \p block to \p fp */
+void OutputFile::writeData(const t_enxblock &block,
+ int subBlockStart,
+ FILE *fp) const
+{
+ int numPoints = block.sub[subBlockStart + 1].nr;
+ for (int j = 0; j < numPoints; j++)
+ {
+ /* Print the coordinates for numDim dimensions */
+ for (int d = 0; d < numDim_; d++)
+ {
+ fprintf(fp, " %8.4f", block.sub[subBlockStart + 1 + d].fval[j]);
+ }
+
+ /* Print numGraph observables */
+ for (int i = 0; i < numGraph_; i++)
+ {
+ fprintf(fp, " %g", block.sub[firstGraphSubblock_ + i].fval[j]*scaleFactor_[i]);
+ }
+
+ fprintf(fp, "\n");
+ }
+}
+
+void AwhReader::processAwhFrame(const t_enxblock &block,
+ double time,
+ const gmx_output_env_t *oenv) const
+{
+ /* We look for AWH data every energy frame and count the no of AWH frames found. We only extract every 'skip' AWH frame. */
+
+ for (auto &biasReader : biasReader_)
+ {
+ /* 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. */
+ fprintf(fpAwh, "# AWH metadata: target error = %.2f kT = %.2f kJ/mol\n",
+ block.sub[subStart].fval[1],
+ block.sub[subStart].fval[1]*kT_);
+
+ fprintf(fpAwh, "# AWH metadata: log sample weight = %4.2f\n",
+ block.sub[subStart].fval[2]);
+
+ awhOutputFile.writeData(block, subStart, fpAwh);
+
+ gmx_ffclose(fpAwh);
+ }
+ }
+}
+
+/*! \brief The main function for the AWH tool */
+int gmx_awh(int argc, char *argv[])
+{
+ const char *desc[] = {
+ "[THISMODULE] extracts AWH data from an energy file.",
+ "One file is 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.",
+ };
+ static gmx_bool moreGraphs = FALSE;
+ static int skip = 0;
+ static gmx_bool kTUnit = FALSE;
+ t_pargs pa[] = {
+ { "-skip", FALSE, etINT, {&skip},
+ "Skip number of frames between data points" },
+ { "-more", FALSE, etBOOL, {&moreGraphs},
+ "Print more output" },
+ { "-kt", FALSE, etBOOL, {&kTUnit},
+ "Print free energy output in units of kT instead of kJ/mol" }
+ };
+
+ ener_file_t fp;
+ t_inputrec ir;
+ gmx_enxnm_t *enm = nullptr;
+ t_enxframe *frame;
+ int nre;
+ gmx_output_env_t *oenv;
+
+ t_filenm fnm[] = {
+ { efEDR, "-f", nullptr, ffREAD },
+ { efTPR, "-s", nullptr, ffREAD },
+ { efXVG, "-o", "awh", ffWRITE }
+ };
+ const int nfile = asize(fnm);
+ if (!parse_common_args(&argc, argv,
+ PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END,
+ nfile, fnm, asize(pa), pa, asize(desc), desc, 0, nullptr, &oenv))
+ {
+ return 0;
+ }
+
+ snew(frame, 1);
+ fp = open_enx(ftp2fn(efEDR, nfile, fnm), "r");
+ do_enxnms(fp, &nre, &enm);
+
+ /* We just need the AWH parameters from inputrec. These are used to initialize
+ the AWH reader when we have a frame to read later on. */
+ gmx_mtop_t mtop;
+ int natoms;
+ matrix box;
+ read_tpx(ftp2fn(efTPR, nfile, fnm), &ir, box, &natoms, nullptr, nullptr, &mtop);
+
+ if (!ir.bDoAwh)
+ {
+ gmx_fatal(FARGS, "No AWH data in %s\n", opt2fn("-f", nfile, fnm));
+ }
+
+ std::unique_ptr<AwhReader> awhReader;
+
+ /* Initiate counters */
+ gmx_bool haveFrame;
+ int awhFrameCounter = 0;
+ int timeCheck = 0;
+ do
+ {
+ haveFrame = do_enx(fp, frame);
+
+ bool useFrame = false;
+
+ t_enxblock *block = nullptr;
+
+ if (haveFrame)
+ {
+ timeCheck = check_times(frame->t);
+
+ if (timeCheck == 0)
+ {
+ block = find_block_id_enxframe(frame, enxAWH, nullptr);
+
+ if (block != nullptr)
+ {
+ if ((skip == 0) || (awhFrameCounter % skip == 0))
+ {
+ useFrame = true;
+ }
+ awhFrameCounter++;
+ }
+ }
+ }
+
+ if (useFrame)
+ {
+ /* We read a valid frame with an AWH block, so we can use it */
+
+ /* Currently we have the number of subblocks per AWH stored
+ * in the first subblock (we cannot get this directly from the tpr),
+ * so we have to read an AWH block before making the legends.
+ */
+ if (awhReader == nullptr)
+ {
+ AwhGraphSelection awhGraphSelection = (moreGraphs ? AwhGraphSelection::All : AwhGraphSelection::Pmf);
+ EnergyUnit energyUnit = (kTUnit ? EnergyUnit::KT : EnergyUnit::KJPerMol);
+ awhReader =
+ std::unique_ptr<AwhReader>(new AwhReader(ir.awhParams,
+ nfile, fnm,
+ awhGraphSelection,
+ energyUnit, BOLTZ*ir.opts.ref_t[0],
+ block));
+ }
+
+ awhReader->processAwhFrame(*block, frame->t, oenv);
+ }
+ }
+ while (haveFrame && (timeCheck <= 0));
+
+ fprintf(stderr, "\n");
+
+ close_enx(fp);
+
+ return 0;
+}
#include <stdlib.h>
#include <string.h>
+#include "gromacs/awh/awh.h"
#include "gromacs/fileio/enxio.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/xvgr.h"
gmx_int64_t step, double time,
int mode,
t_mdebin *md, t_fcdata *fcd,
- gmx_groups_t *groups, t_grpopts *opts)
+ gmx_groups_t *groups, t_grpopts *opts,
+ gmx::Awh *awh)
{
/*static char **grpnms=NULL;*/
char buf[246];
mde_delta_h_coll_reset(md->dhc);
}
+ /* AWH bias blocks. */
+ if (awh != nullptr) // TODO: add boolean in t_mdebin. See in mdebin.h.
+ {
+ awh->writeToEnergyFrame(step, &fr);
+ }
+
/* do the actual I/O */
do_enx(fp_ene, &fr);
if (fr.nre)
struct t_lambda;
class t_state;
+namespace gmx
+{
+class Awh;
+}
+
/* The functions & data structures here determine the content for outputting
the .edr file; the file format and actual writing is done with functions
defined in enxio.h */
gmx_int64_t step, double time,
int mode,
t_mdebin *md, t_fcdata *fcd,
- gmx_groups_t *groups, t_grpopts *opts);
+ gmx_groups_t *groups, t_grpopts *opts,
+ gmx::Awh *awh);
if (ir->eI == eiNM)
{
- GMX_ASSERT(shellfc != NULL, "With NM we always support shells");
+ GMX_ASSERT(shellfc != nullptr, "With NM we always support shells");
*shellfc = init_shell_flexcon(stdout,
top_global,
print_ebin_header(fplog, step, step);
print_ebin(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, step, eprNORMAL,
- mdebin, fcd, &(top_global->groups), &(inputrec->opts));
+ mdebin, fcd, &(top_global->groups), &(inputrec->opts), nullptr);
}
where();
}
print_ebin(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
do_log ? fplog : nullptr, step, step, eprNORMAL,
- mdebin, fcd, &(top_global->groups), &(inputrec->opts));
+ mdebin, fcd, &(top_global->groups), &(inputrec->opts), nullptr);
}
/* Send energies and positions to the IMD client if bIMD is TRUE. */
/* Write final energy file entries */
print_ebin(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
!do_log ? fplog : nullptr, step, step, eprNORMAL,
- mdebin, fcd, &(top_global->groups), &(inputrec->opts));
+ mdebin, fcd, &(top_global->groups), &(inputrec->opts), nullptr);
}
}
print_ebin_header(fplog, step, step);
print_ebin(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, step, eprNORMAL,
- mdebin, fcd, &(top_global->groups), &(inputrec->opts));
+ mdebin, fcd, &(top_global->groups), &(inputrec->opts), nullptr);
}
where();
}
print_ebin(mdoutf_get_fp_ene(outf), do_ene, FALSE, FALSE,
do_log ? fplog : nullptr, step, step, eprNORMAL,
- mdebin, fcd, &(top_global->groups), &(inputrec->opts));
+ mdebin, fcd, &(top_global->groups), &(inputrec->opts), nullptr);
}
/* Send x and E to IMD client, if bIMD is TRUE. */
{
print_ebin(mdoutf_get_fp_ene(outf), !do_ene, FALSE, FALSE,
!do_log ? fplog : nullptr, step, step, eprNORMAL,
- mdebin, fcd, &(top_global->groups), &(inputrec->opts));
+ mdebin, fcd, &(top_global->groups), &(inputrec->opts), nullptr);
}
/* Print some stuff... */
do_per_step(steps_accepted, inputrec->nstdisreout),
do_per_step(steps_accepted, inputrec->nstorireout),
fplog, count, count, eprNORMAL,
- mdebin, fcd, &(top_global->groups), &(inputrec->opts));
+ mdebin, fcd, &(top_global->groups), &(inputrec->opts), nullptr);
fflush(fplog);
}
}
"Analyze data sets");
registerModule(manager, &gmx_g_angle, "angle",
"Calculate distributions and correlations for angles and dihedrals");
+ registerModule(manager, &gmx_awh, "awh",
+ "Extract data from an accelerated weight histogram (AWH) run");
registerModule(manager, &gmx_bar, "bar",
"Calculate free energy difference estimates through Bennett's acceptance ratio");
registerObsoleteTool(manager, "bond");
gmx::CommandLineModuleGroup group =
manager->addModuleGroup("Tools");
group.addModule("analyze");
+ group.addModule("awh");
group.addModule("dyndom");
group.addModule("filter");
group.addModule("lie");
print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
step, t,
- eprNORMAL, mdebin, fcd, groups, &(ir->opts));
+ eprNORMAL, mdebin, fcd, groups, &(ir->opts), ir->awh);
if (ir->bPull)
{
if (ir->nstcalcenergy > 0 && !bRerunMD)
{
print_ebin(mdoutf_get_fp_ene(outf), FALSE, FALSE, FALSE, fplog, step, t,
- eprAVER, mdebin, fcd, groups, &(ir->opts));
+ eprAVER, mdebin, fcd, groups, &(ir->opts), ir->awh);
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff