All structs have been converted to classes with RAII constructors.
Parameter objects can only be constructed from input data or
serializers.
Tests are adapted to use the new constructors, this means some code got
a bit ugly there.
This change is only refactoring, plus adding tests for the
serialization.
#include <cstring>
#include <algorithm>
+#include <vector>
#include "gromacs/fileio/enxio.h"
#include "gromacs/gmxlib/network.h"
*/
static bool anyDimUsesPull(const AwhBiasParams& awhBiasParams)
{
- for (int d = 0; d < awhBiasParams.ndim; d++)
+ for (const auto& awhDimParam : awhBiasParams.dimParams())
{
- const AwhDimParams& awhDimParams = awhBiasParams.dimParams[d];
- if (awhDimParams.eCoordProvider == AwhCoordinateProviderType::Pull)
+ if (awhDimParam.coordinateProvider() == AwhCoordinateProviderType::Pull)
{
return true;
}
*/
static bool anyDimUsesPull(const AwhParams& awhParams)
{
- for (int k = 0; k < awhParams.numBias; k++)
+ for (const auto& awhBiasParam : awhParams.awhBiasParams())
{
- const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[k];
- if (anyDimUsesPull(awhBiasParams))
+ if (anyDimUsesPull(awhBiasParam))
{
return true;
}
pull_t* pull_work,
int numFepLambdaStates,
int fepLambdaState) :
- seed_(awhParams.seed),
- nstout_(awhParams.nstOut),
+ seed_(awhParams.seed()),
+ nstout_(awhParams.nstout()),
commRecord_(commRecord),
multiSimRecord_(multiSimRecord),
pull_(pull_work),
}
int numSharingSimulations = 1;
- if (awhParams.shareBiasMultisim && isMultiSim(multiSimRecord_))
+ if (awhParams.shareBiasMultisim() && isMultiSim(multiSimRecord_))
{
numSharingSimulations = multiSimRecord_->numSimulations_;
}
/* Initialize all the biases */
- const double beta = 1 / (gmx::c_boltz * inputRecord.opts.ref_t[0]);
- for (int k = 0; k < awhParams.numBias; k++)
+ const double beta = 1 / (gmx::c_boltz * inputRecord.opts.ref_t[0]);
+ const auto& awhBiasParams = awhParams.awhBiasParams();
+ for (int k = 0; k < gmx::ssize(awhBiasParams); k++)
{
- const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[k];
-
std::vector<int> pullCoordIndex;
std::vector<DimParams> dimParams;
- for (int d = 0; d < awhBiasParams.ndim; d++)
+ for (const auto& awhDimParam : awhBiasParams[k].dimParams())
{
- const AwhDimParams& awhDimParams = awhBiasParams.dimParams[d];
- if (awhDimParams.eCoordProvider != AwhCoordinateProviderType::Pull
- && awhDimParams.eCoordProvider != AwhCoordinateProviderType::FreeEnergyLambda)
+ if (awhDimParam.coordinateProvider() != AwhCoordinateProviderType::Pull
+ && awhDimParam.coordinateProvider() != AwhCoordinateProviderType::FreeEnergyLambda)
{
GMX_THROW(
InvalidInputError("Currently only the pull code and lambda are supported "
"as coordinate providers"));
}
- if (awhDimParams.eCoordProvider == AwhCoordinateProviderType::Pull)
+ if (awhDimParam.coordinateProvider() == AwhCoordinateProviderType::Pull)
{
- const t_pull_coord& pullCoord = inputRecord.pull->coord[awhDimParams.coordIndex];
+ const t_pull_coord& pullCoord = inputRecord.pull->coord[awhDimParam.coordinateIndex()];
if (pullCoord.eGeom == PullGroupGeometry::DirectionPBC)
{
GMX_THROW(InvalidInputError(
"Pull geometry 'direction-periodic' is not supported by AWH"));
}
double conversionFactor = pull_coordinate_is_angletype(&pullCoord) ? gmx::c_deg2Rad : 1;
- pullCoordIndex.push_back(awhDimParams.coordIndex);
- dimParams.push_back(DimParams::pullDimParams(
- conversionFactor, awhDimParams.forceConstant, beta));
+ pullCoordIndex.push_back(awhDimParam.coordinateIndex());
+ dimParams.emplace_back(DimParams::pullDimParams(
+ conversionFactor, awhDimParam.forceConstant(), beta));
}
else
{
(MASTER(commRecord_) ? Bias::ThisRankWillDoIO::Yes : Bias::ThisRankWillDoIO::No);
biasCoupledToSystem_.emplace_back(Bias(k,
awhParams,
- awhParams.awhBiasParams[k],
+ awhBiasParams[k],
dimParams,
beta,
inputRecord.delta_t,
{
GMX_RELEASE_ASSERT(!anyDimUsesPull(awhParams) || pull_work, "Need a valid pull object");
- for (int k = 0; k < awhParams.numBias; k++)
+ for (const auto& biasParam : awhParams.awhBiasParams())
{
- const AwhBiasParams& biasParams = awhParams.awhBiasParams[k];
-
- for (int d = 0; d < biasParams.ndim; d++)
+ for (const auto& dimParam : biasParam.dimParams())
{
- if (biasParams.dimParams[d].eCoordProvider == AwhCoordinateProviderType::Pull)
+ if (dimParam.coordinateProvider() == AwhCoordinateProviderType::Pull)
{
register_external_pull_potential(
- pull_work, biasParams.dimParams[d].coordIndex, Awh::externalPotentialString());
+ pull_work, dimParam.coordinateIndex(), Awh::externalPotentialString());
}
}
}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2015,2016,2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
template<typename>
class ArrayRef;
struct AwhHistory;
-struct AwhParams;
+class AwhParams;
class Bias;
struct BiasCoupledToSystem;
class ForceWithVirial;
ThisRankWillDoIO thisRankWillDoIO,
BiasParams::DisableUpdateSkips disableUpdateSkips) :
dimParams_(dimParamsInit.begin(), dimParamsInit.end()),
- grid_(dimParamsInit, awhBiasParams.dimParams),
+ grid_(dimParamsInit, awhBiasParams.dimParams()),
params_(awhParams,
awhBiasParams,
dimParams_,
/* For a global update updateList covers all points, so reserve that */
updateList_.reserve(grid_.numPoints());
- state_.initGridPointState(awhBiasParams, dimParams_, grid_, params_, biasInitFilename, awhParams.numBias);
+ state_.initGridPointState(
+ awhBiasParams, dimParams_, grid_, params_, biasInitFilename, awhParams.numBias());
if (thisRankDoesIO_)
{
ndim(),
blockLength,
CorrelationGrid::BlockLengthMeasure::Time,
- awhParams.nstSampleCoord * mdTimeStep);
+ awhParams.nstSampleCoord() * mdTimeStep);
writer_ = std::make_unique<BiasWriter>(*this);
}
{
struct AwhBiasHistory;
-struct AwhBiasParams;
+class AwhBiasParams;
struct AwhHistory;
-struct AwhParams;
+class AwhParams;
struct AwhPointStateHistory;
class CorrelationGrid;
class BiasGrid;
{
/*! \brief
- * Modify x so that it is periodic in [-period/2, +period/2).
+ * Return x so that it is periodic in [-period/2, +period/2).
*
* x is modified by shifting its value by a +/- a period if
* needed. Thus, it is assumed that x is at most one period
* away from this interval. For period = 0, x is not modified.
*
- * \param[in,out] x Pointer to the value to modify.
- * \param[in] period The period, or 0 if not periodic.
+ * \param[in] x Pointer to the value to modify.
+ * \param[in] period The period, or 0 if not periodic.
+ * \returns Value that is within the period.
*/
-void centerPeriodicValueAroundZero(double* x, double period)
+double centerPeriodicValueAroundZero(const double x, double period)
{
GMX_ASSERT(period >= 0, "Periodic should not be negative");
const double halfPeriod = period * 0.5;
- if (*x >= halfPeriod)
+ double valueInPeriod = x;
+
+ if (valueInPeriod >= halfPeriod)
{
- *x -= period;
+ valueInPeriod -= period;
}
- else if (*x < -halfPeriod)
+ else if (valueInPeriod < -halfPeriod)
{
- *x += period;
+ valueInPeriod += period;
}
+ return valueInPeriod;
}
/*! \brief
if (period > 0)
{
- centerPeriodicValueAroundZero(&dev, period);
+ dev = centerPeriodicValueAroundZero(dev, period);
}
return dev;
if (axis_[d].period() > 0)
{
/* Do we always want the values to be centered around 0 ? */
- centerPeriodicValueAroundZero(&point.coordValue[d], axis_[d].period());
+ point.coordValue[d] =
+ centerPeriodicValueAroundZero(point.coordValue[d], axis_[d].period());
}
point.index[d] = indexWork[d];
}
}
-BiasGrid::BiasGrid(ArrayRef<const DimParams> dimParams, const AwhDimParams* awhDimParams)
+BiasGrid::BiasGrid(ArrayRef<const DimParams> dimParams, ArrayRef<const AwhDimParams> awhDimParams)
{
+ GMX_RELEASE_ASSERT(dimParams.size() == awhDimParams.size(), "Dimensions needs to be equal");
/* Define the discretization along each dimension */
awh_dvec period;
int numPoints = 1;
- for (size_t d = 0; d < dimParams.size(); d++)
+ for (int d = 0; d < gmx::ssize(awhDimParams); d++)
{
- double origin = dimParams[d].scaleUserInputToInternal(awhDimParams[d].origin);
- double end = dimParams[d].scaleUserInputToInternal(awhDimParams[d].end);
- if (awhDimParams[d].eCoordProvider == AwhCoordinateProviderType::Pull)
+ double origin = dimParams[d].scaleUserInputToInternal(awhDimParams[d].origin());
+ double end = dimParams[d].scaleUserInputToInternal(awhDimParams[d].end());
+ if (awhDimParams[d].coordinateProvider() == AwhCoordinateProviderType::Pull)
{
- period[d] = dimParams[d].scaleUserInputToInternal(awhDimParams[d].period);
+ period[d] = dimParams[d].scaleUserInputToInternal(awhDimParams[d].period());
static_assert(
c_numPointsPerSigma >= 1.0,
"The number of points per sigma should be at least 1.0 to get a uniformly "
namespace gmx
{
-struct AwhDimParams;
+class AwhDimParams;
/*! \internal
* \brief An axis, i.e. dimension, of the grid.
* coordinate living on the grid (determines the grid spacing).
* \param[in] awhDimParams Dimension params from inputrec.
*/
- BiasGrid(ArrayRef<const DimParams> dimParams, const AwhDimParams* awhDimParams);
+ BiasGrid(ArrayRef<const DimParams> dimParams, ArrayRef<const AwhDimParams> awhDimParams);
/*! \brief Returns the number of points in the grid.
*
* (this could be made a user-option but there is most likely no big need
* for tweaking this for most users).
*/
- switch (awhBiasParams.eTarget)
+ switch (awhBiasParams.targetDistribution())
{
case AwhTargetType::Constant: numStepsUpdateTarget = 0; break;
case AwhTargetType::Cutoff:
/* Updating the target generally requires updating the whole grid so to keep the cost
down we generally update the target less often than the free energy (unless the free
energy update step is set to > 100 samples). */
- numStepsUpdateTarget = std::max(100 % awhParams.numSamplesUpdateFreeEnergy,
- awhParams.numSamplesUpdateFreeEnergy)
- * awhParams.nstSampleCoord;
+ numStepsUpdateTarget = std::max(100 % awhParams.numSamplesUpdateFreeEnergy(),
+ awhParams.numSamplesUpdateFreeEnergy())
+ * awhParams.nstSampleCoord();
break;
case AwhTargetType::LocalBoltzmann:
/* The target distribution is set equal to the reference histogram which is updated every free energy update.
target distribution. One could avoid the global update by making a local target update function (and
postponing target updates for non-local points as for the free energy update). We avoid such additions
for now and accept that this target type always does global updates. */
- numStepsUpdateTarget = awhParams.numSamplesUpdateFreeEnergy * awhParams.nstSampleCoord;
+ numStepsUpdateTarget = awhParams.numSamplesUpdateFreeEnergy() * awhParams.nstSampleCoord();
break;
default: GMX_RELEASE_ASSERT(false, "Unknown AWH target type"); break;
}
/* Convert to number of steps using the sampling frequency. The
check interval should be a multiple of the update step
interval. */
- int numStepsUpdate = awhParams.numSamplesUpdateFreeEnergy * awhParams.nstSampleCoord;
- GMX_RELEASE_ASSERT(awhParams.numSamplesUpdateFreeEnergy > 0,
+ int numStepsUpdate = awhParams.numSamplesUpdateFreeEnergy() * awhParams.nstSampleCoord();
+ GMX_RELEASE_ASSERT(awhParams.numSamplesUpdateFreeEnergy() > 0,
"When checking for AWH coverings, the number of samples per AWH update need "
"to be > 0.");
- int numUpdatesCheck = std::max(1, minNumSamplesCover / awhParams.numSamplesUpdateFreeEnergy);
+ int numUpdatesCheck = std::max(1, minNumSamplesCover / awhParams.numSamplesUpdateFreeEnergy());
int numStepsCheck = numUpdatesCheck * numStepsUpdate;
GMX_RELEASE_ASSERT(numStepsCheck % numStepsUpdate == 0,
/* Get diffusion factor */
double maxCrossingTime = 0.;
std::vector<double> x;
+ const auto awhDimParams = awhBiasParams.dimParams();
for (size_t d = 0; d < gridAxis.size(); d++)
{
- GMX_RELEASE_ASSERT(awhBiasParams.dimParams[d].diffusion > 0, "We need positive diffusion");
+ GMX_RELEASE_ASSERT(awhDimParams[d].diffusion() > 0, "We need positive diffusion");
// With diffusion it takes on average T = L^2/2D time to cross length L
double axisLength = gridAxis[d].isFepLambdaAxis() ? 1.0 : gridAxis[d].length();
- double crossingTime = (axisLength * axisLength) / (2 * awhBiasParams.dimParams[d].diffusion);
+ double crossingTime = (axisLength * axisLength) / (2 * awhDimParams[d].diffusion());
maxCrossingTime = std::max(maxCrossingTime, crossingTime);
}
GMX_RELEASE_ASSERT(maxCrossingTime > 0, "We need at least one dimension with non-zero length");
- double errorInitialInKT = beta * awhBiasParams.errorInitial;
+ double errorInitialInKT = beta * awhBiasParams.initialErrorEstimate();
double histogramSize = maxCrossingTime / (gmx::square(errorInitialInKT) * samplingTimestep);
return histogramSize;
int numShared = 1;
- if (awhBiasParams.shareGroup > 0)
+ if (awhBiasParams.shareGroup() > 0)
{
/* We do not yet support sharing within a simulation */
int numSharedWithinThisSimulation = 1;
ArrayRef<const GridAxis> gridAxis,
int biasIndex) :
invBeta(beta > 0 ? 1 / beta : 0),
- numStepsSampleCoord_(awhParams.nstSampleCoord),
- numSamplesUpdateFreeEnergy_(awhParams.numSamplesUpdateFreeEnergy),
+ numStepsSampleCoord_(awhParams.nstSampleCoord()),
+ numSamplesUpdateFreeEnergy_(awhParams.numSamplesUpdateFreeEnergy()),
numStepsUpdateTarget_(calcTargetUpdateInterval(awhParams, awhBiasParams)),
numStepsCheckCovering_(calcCheckCoveringInterval(awhParams, dimParams, gridAxis)),
- eTarget(awhBiasParams.eTarget),
- freeEnergyCutoffInKT(beta * awhBiasParams.targetCutoff),
- temperatureScaleFactor(awhBiasParams.targetBetaScaling),
+ eTarget(awhBiasParams.targetDistribution()),
+ freeEnergyCutoffInKT(beta * awhBiasParams.targetCutoff()),
+ temperatureScaleFactor(awhBiasParams.targetBetaScaling()),
idealWeighthistUpdate(eTarget != AwhTargetType::LocalBoltzmann),
numSharedUpdate(getNumSharedUpdate(awhBiasParams, numSharingSimulations)),
updateWeight(numSamplesUpdateFreeEnergy_ * numSharedUpdate),
localWeightScaling(eTarget == AwhTargetType::LocalBoltzmann ? temperatureScaleFactor : 1),
- initialErrorInKT(beta * awhBiasParams.errorInitial),
+ initialErrorInKT(beta * awhBiasParams.initialErrorEstimate()),
initialHistogramSize(
getInitialHistogramSizeEstimate(awhBiasParams, gridAxis, beta, numStepsSampleCoord_ * mdTimeStep)),
- convolveForce(awhParams.ePotential == AwhPotentialType::Convolved),
+ convolveForce(awhParams.potential() == AwhPotentialType::Convolved),
biasIndex(biasIndex),
disableUpdateSkips_(disableUpdateSkips == DisableUpdateSkips::yes)
{
GMX_THROW(InvalidInputError("To use AWH, the beta=1/(k_B T) should be > 0"));
}
- for (int d = 0; d < awhBiasParams.ndim; d++)
+ const auto& awhDimParams = awhBiasParams.dimParams();
+ for (int d = 0; d < gmx::ssize(awhDimParams); d++)
{
/* The spacing in FEP dimensions is 1. The calculated coverRadius will be in lambda states
* (cf points in other dimensions). */
- double coverRadiusInNm = 0.5 * awhBiasParams.dimParams[d].coverDiameter;
+ double coverRadiusInNm = 0.5 * awhDimParams[d].coverDiameter();
double spacing = gridAxis[d].spacing();
coverRadius_[d] = spacing > 0 ? static_cast<int>(std::round(coverRadiusInNm / spacing)) : 0;
}
template<typename>
class ArrayRef;
-struct AwhBiasParams;
-struct AwhParams;
+class AwhBiasParams;
+class AwhParams;
struct DimParams;
class GridAxis;
enum class AwhTargetType : int;
#include "gromacs/mdrunutility/multisim.h"
#include "gromacs/mdtypes/awh_params.h"
#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/stringutil.h"
{
bool haveSharing = false;
- for (int k = 0; k < awhParams.numBias; k++)
+ const auto awhBiasParams = awhParams.awhBiasParams();
+ for (auto awhBiasParamIt = awhBiasParams.begin(); awhBiasParamIt != awhBiasParams.end(); ++awhBiasParamIt)
{
- int shareGroup = awhParams.awhBiasParams[k].shareGroup;
+ int shareGroup = awhBiasParamIt->shareGroup();
if (shareGroup > 0)
{
- for (int i = k + 1; i < awhParams.numBias; i++)
+ for (auto awhBiasParamIt2 = awhBiasParamIt + 1; awhBiasParamIt2 != awhBiasParams.end();
+ ++awhBiasParamIt2)
{
- if (awhParams.awhBiasParams[i].shareGroup == shareGroup)
+ if (awhBiasParamIt2->shareGroup() == shareGroup)
{
haveSharing = true;
}
* biases in order over the ranks and it does not restrict possibilities.
*/
int numShare = 0;
- for (int b = 0; b < awhParams.numBias; b++)
+ for (const auto& awhBiasParam : awhParams.awhBiasParams())
{
- int group = awhParams.awhBiasParams[b].shareGroup;
+ int group = awhBiasParam.shareGroup();
if (group > 0)
{
numShare++;
}
std::vector<int> intervals(numSim * 2);
- intervals[numSim * 0 + multiSimComm->simulationIndex_] = awhParams.nstSampleCoord;
- intervals[numSim * 1 + multiSimComm->simulationIndex_] = awhParams.numSamplesUpdateFreeEnergy;
+ intervals[numSim * 0 + multiSimComm->simulationIndex_] = awhParams.nstSampleCoord();
+ intervals[numSim * 1 + multiSimComm->simulationIndex_] = awhParams.numSamplesUpdateFreeEnergy();
gmx_sumi_sim(intervals.size(), intervals.data(), multiSimComm);
for (int sim = 1; sim < numSim; sim++)
{
/* Check the point sizes. This is a sufficient condition for running
* as shared multi-sim run. No physics checks are performed here.
*/
- for (int b = 0; b < awhParams.numBias; b++)
+ const auto& awhBiasParams = awhParams.awhBiasParams();
+ for (int b = 0; b < gmx::ssize(awhBiasParams); b++)
{
- if (awhParams.awhBiasParams[b].shareGroup > 0)
+ if (awhBiasParams[b].shareGroup() > 0)
{
std::vector<int64_t> pointSizes(numSim);
pointSizes[multiSimComm->simulationIndex_] = pointSize[b];
template<typename>
class ArrayRef;
-struct AwhParams;
+class AwhParams;
/*! \brief Returns if any bias is sharing within a simulation.
*
/* Modify PMF, free energy and the constant target distribution factor
* to user input values if there is data given.
*/
- if (awhBiasParams.bUserData)
+ if (awhBiasParams.userPMFEstimate())
{
readUserPmfAndTargetDistribution(dimParams, grid, filename, numBias, params.biasIndex, &points_);
setFreeEnergyToConvolvedPmf(dimParams, grid);
template<typename>
class ArrayRef;
struct AwhBiasHistory;
-struct AwhBiasParams;
class BiasParams;
class BiasGrid;
class GridAxis;
ArrayRef<const DimParams> dimParams,
const BiasGrid& grid)
{
- for (size_t d = 0; d < dimParams.size(); d++)
+ GMX_RELEASE_ASSERT(awhBiasParams.ndim() == dimParams.ssize(),
+ "Need to have identical size for dimensions");
+ const auto& awhDimParams = awhBiasParams.dimParams();
+ for (int d = 0; d < gmx::ssize(awhDimParams); d++)
{
- coordValue_[d] = dimParams[d].scaleUserInputToInternal(awhBiasParams.dimParams[d].coordValueInit);
+ coordValue_[d] = dimParams[d].scaleUserInputToInternal(awhDimParams[d].initialCoordinate());
}
/* The grid-point index is always the nearest point to the coordinate.
template<typename>
class ArrayRef;
-struct AwhBiasParams;
+class AwhBiasParams;
struct AwhBiasStateHistory;
class BiasParams;
class BiasGrid;
HistogramSize::HistogramSize(const AwhBiasParams& awhBiasParams, double histogramSizeInitial) :
numUpdates_(0),
histogramSize_(histogramSizeInitial),
- inInitialStage_(awhBiasParams.eGrowth == AwhHistogramGrowthType::ExponentialLinear),
- equilibrateHistogram_(awhBiasParams.equilibrateHistogram),
+ inInitialStage_(awhBiasParams.growthType() == AwhHistogramGrowthType::ExponentialLinear),
+ equilibrateHistogram_(awhBiasParams.equilibrateHistogram()),
logScaledSampleWeight_(0),
maxLogScaledSampleWeight_(0),
havePrintedAboutCovering_(false)
template<typename>
class ArrayRef;
struct AwhBiasStateHistory;
-struct AwhBiasParams;
+class AwhBiasParams;
class BiasParams;
class PointState;
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/enumerationhelpers.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/iserializer.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/stringutil.h"
bool usesPull = false;
bool usesFep = false;
- for (int b = 0; b < inputrec.awhParams->numBias; b++)
+ for (const auto& awhBiasParam : inputrec.awhParams->awhBiasParams())
{
- const auto& biasParams = inputrec.awhParams->awhBiasParams[b];
- for (int d = 0; d < biasParams.ndim; d++)
+ for (const auto& dimParam : awhBiasParam.dimParams())
{
- switch (biasParams.dimParams[d].eCoordProvider)
+ switch (dimParam.coordinateProvider())
{
case AwhCoordinateProviderType::Pull: usesPull = true; break;
case AwhCoordinateProviderType::FreeEnergyLambda: usesFep = true; break;
"computed at the slow MTS level");
}
- if (inputrec.awhParams->nstSampleCoord % inputrec.mtsLevels[awhMtsLevel].stepFactor != 0)
+ if (inputrec.awhParams->nstSampleCoord() % inputrec.mtsLevels[awhMtsLevel].stepFactor != 0)
{
warning_error(wi,
"With MTS applied to AWH, awh-nstsample should be a multiple of mts-factor");
* \param[in,out] wi Struct for bookeeping warnings.
*/
void checkPullDimParams(const std::string& prefix,
- AwhDimParams* dimParams,
+ const AwhDimParams& dimParams,
const pull_params_t& pull_params,
warninp_t wi)
{
- if (dimParams->coordIndex < 0)
+ if (dimParams.coordinateIndex() < 0)
{
gmx_fatal(FARGS,
"Failed to read a valid coordinate index for %s-coord-index. "
"Note that the pull coordinate indexing starts at 1.",
prefix.c_str());
}
- if (dimParams->coordIndex >= pull_params.ncoord)
+ if (dimParams.coordinateIndex() >= pull_params.ncoord)
{
gmx_fatal(FARGS,
"The given AWH coordinate index (%d) is larger than the number of pull "
"coordinates (%d)",
- dimParams->coordIndex + 1,
+ dimParams.coordinateIndex() + 1,
pull_params.ncoord);
}
- if (pull_params.coord[dimParams->coordIndex].rate != 0)
+ if (pull_params.coord[dimParams.coordinateIndex()].rate != 0)
{
auto message = formatString(
"Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate",
- dimParams->coordIndex + 1,
- pull_params.coord[dimParams->coordIndex].rate);
+ dimParams.coordinateIndex() + 1,
+ pull_params.coord[dimParams.coordinateIndex()].rate);
warning_error(wi, message);
}
- if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
+ if (gmx_within_tol(dimParams.end() - dimParams.origin(), 0, GMX_REAL_EPS))
{
auto message = formatString(
"The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
"This will result in only one point along this axis in the coordinate value grid.",
prefix.c_str(),
- dimParams->origin,
+ dimParams.origin(),
prefix.c_str(),
- dimParams->end);
+ dimParams.end());
warning(wi, message);
}
- if (dimParams->forceConstant <= 0)
+ if (dimParams.forceConstant() <= 0)
{
warning_error(wi, "The force AWH bias force constant should be > 0");
}
/* Grid params for each axis */
- PullGroupGeometry eGeom = pull_params.coord[dimParams->coordIndex].eGeom;
+ PullGroupGeometry eGeom = pull_params.coord[dimParams.coordinateIndex()].eGeom;
/* Check that the requested interval is in allowed range */
if (eGeom == PullGroupGeometry::Distance)
{
- if (dimParams->origin < 0 || dimParams->end < 0)
+ if (dimParams.origin() < 0 || dimParams.end() < 0)
{
gmx_fatal(FARGS,
"%s-start (%g) or %s-end (%g) set to a negative value. With pull "
"geometry distance coordinate values are non-negative. "
"Perhaps you want to use geometry %s instead?",
prefix.c_str(),
- dimParams->origin,
+ dimParams.origin(),
prefix.c_str(),
- dimParams->end,
+ dimParams.end(),
enumValueToString(PullGroupGeometry::Direction));
}
}
else if (eGeom == PullGroupGeometry::Angle || eGeom == PullGroupGeometry::AngleAxis)
{
- if (dimParams->origin < 0 || dimParams->end > 180)
+ if (dimParams.origin() < 0 || dimParams.end() > 180)
{
gmx_fatal(FARGS,
"%s-start (%g) and %s-end (%g) are outside of the allowed range "
"0 to 180 deg for pull geometries %s and %s ",
prefix.c_str(),
- dimParams->origin,
+ dimParams.origin(),
prefix.c_str(),
- dimParams->end,
+ dimParams.end(),
enumValueToString(PullGroupGeometry::Angle),
enumValueToString(PullGroupGeometry::AngleAxis));
}
}
else if (eGeom == PullGroupGeometry::Dihedral)
{
- if (dimParams->origin < -180 || dimParams->end > 180)
+ if (dimParams.origin() < -180 || dimParams.end() > 180)
{
gmx_fatal(FARGS,
"%s-start (%g) and %s-end (%g) are outside of the allowed range "
"-180 to 180 deg for pull geometry %s. ",
prefix.c_str(),
- dimParams->origin,
+ dimParams.origin(),
prefix.c_str(),
- dimParams->end,
+ dimParams.end(),
enumValueToString(PullGroupGeometry::Dihedral));
}
}
* \param[in,out] wi Struct for bookeeping warnings.
*/
void checkFepLambdaDimParams(const std::string& prefix,
- const AwhDimParams* dimParams,
+ const AwhDimParams& dimParams,
const t_lambda* lambdaParams,
const FreeEnergyPerturbationType efep,
warninp_t wi)
"AWH is not treated like other expanded ensemble methods. Do not use expanded.");
}
- if (dimParams->origin < 0)
+ if (dimParams.origin() < 0)
{
opt = prefix + "-start";
gmx_fatal(FARGS,
"When running AWH coupled to the free energy lambda state the lower lambda state "
"for AWH, %s (%.0f), must be >= 0.",
opt.c_str(),
- dimParams->origin);
+ dimParams.origin());
}
- if (dimParams->end >= lambdaParams->n_lambda)
+ if (dimParams.end() >= lambdaParams->n_lambda)
{
opt = prefix + "-end";
gmx_fatal(FARGS,
"When running AWH coupled to the free energy lambda state the upper lambda state "
"for AWH, %s (%.0f), must be < n_lambda (%d).",
opt.c_str(),
- dimParams->origin,
+ dimParams.origin(),
lambdaParams->n_lambda);
}
- if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
+ if (gmx_within_tol(dimParams.end() - dimParams.origin(), 0, GMX_REAL_EPS))
{
auto message = formatString(
"The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
"This will result in only one lambda point along this free energy lambda state "
"axis in the coordinate value grid.",
prefix.c_str(),
- dimParams->origin,
+ dimParams.origin(),
prefix.c_str(),
- dimParams->end);
+ dimParams.end());
warning(wi, message);
}
- if (dimParams->forceConstant != 0)
+ if (dimParams.forceConstant() != 0)
{
warning_error(
wi,
}
}
-/*! \brief
- * Read parameters of an AWH bias dimension.
- *
- * \param[in,out] inp Input file entries.
- * \param[in] prefix Prefix for dimension parameters.
- * \param[in,out] dimParams AWH dimensional parameters.
- * \param[in,out] wi Struct for bookeeping warnings.
- * \param[in] bComment True if comments should be printed.
- */
-void readDimParams(std::vector<t_inpfile>* inp,
- const std::string& prefix,
- AwhDimParams* dimParams,
- warninp_t wi,
- bool bComment)
-{
- std::string opt;
- if (bComment)
- {
- printStringNoNewline(
- inp,
- "The provider of the reaction coordinate, "
- "currently only 'pull' and 'fep-lambda' (free energy lambda state) is supported");
- }
- opt = prefix + "-coord-provider";
- dimParams->eCoordProvider = getEnum<AwhCoordinateProviderType>(inp, opt.c_str(), wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "The coordinate index for this dimension");
- }
- opt = prefix + "-coord-index";
- int coordIndexInput;
- coordIndexInput = get_eint(inp, opt, 1, wi);
-
- /* The pull coordinate indices start at 1 in the input file, at 0 internally */
- dimParams->coordIndex = coordIndexInput - 1;
-
- if (bComment)
- {
- printStringNoNewline(inp, "Start and end values for each coordinate dimension");
- }
- opt = prefix + "-start";
- dimParams->origin = get_ereal(inp, opt, 0., wi);
- opt = prefix + "-end";
- dimParams->end = get_ereal(inp, opt, 0., wi);
-
- if (bComment)
- {
- printStringNoNewline(
- inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
- }
- opt = prefix + "-force-constant";
- dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps, rad^2/ps or ps^-1)");
- }
- opt = prefix + "-diffusion";
- dimParams->diffusion = get_ereal(inp, opt, 0, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp,
- "Diameter that needs to be sampled around a point before it is "
- "considered covered. In FEP dimensions the cover diameter is "
- "specified in lambda states.");
- }
- opt = prefix + "-cover-diameter";
- dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
-}
-
/*! \brief
* Check the parameters of an AWH bias dimension.
*
* \param[in] ir Input parameter struct.
* \param[in,out] wi Struct for bookeeping warnings.
*/
-void checkDimParams(const std::string& prefix, AwhDimParams* dimParams, const t_inputrec* ir, warninp_t wi)
+void checkDimParams(const std::string& prefix, const AwhDimParams& dimParams, const t_inputrec& ir, warninp_t wi)
{
- if (dimParams->eCoordProvider == AwhCoordinateProviderType::Pull)
+ if (dimParams.coordinateProvider() == AwhCoordinateProviderType::Pull)
{
- if (!ir->bPull)
+ if (!ir.bPull)
{
gmx_fatal(FARGS,
"AWH biasing along a pull dimension is only compatible with COM pulling "
"turned on");
}
- checkPullDimParams(prefix, dimParams, *ir->pull, wi);
+ checkPullDimParams(prefix, dimParams, *ir.pull, wi);
}
- else if (dimParams->eCoordProvider == AwhCoordinateProviderType::FreeEnergyLambda)
+ else if (dimParams.coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda)
{
- if (ir->efep == FreeEnergyPerturbationType::No)
+ if (ir.efep == FreeEnergyPerturbationType::No)
{
gmx_fatal(FARGS,
"AWH biasing along a free energy lambda state dimension is only compatible "
"with free energy turned on");
}
- checkFepLambdaDimParams(prefix, dimParams, ir->fepvals.get(), ir->efep, wi);
+ checkFepLambdaDimParams(prefix, dimParams, ir.fepvals.get(), ir.efep, wi);
}
else
{
}
}
-/*! \brief
- * Check consistency of input at the AWH bias level.
- *
- * \param[in] awhBiasParams AWH bias parameters.
- * \param[in,out] wi Struct for bookkeeping warnings.
- */
-void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
-{
- /* Covering diameter and sharing warning. */
- for (int d = 0; d < awhBiasParams.ndim; d++)
- {
- double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
- if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
- {
- warning(wi,
- "The covering diameter is only relevant to set for bias sharing simulations.");
- }
- }
-}
-
-/*! \brief
- * Read parameters of an AWH bias.
- *
- * \param[in,out] inp Input file entries.
- * \param[in,out] awhBiasParams AWH dimensional parameters.
- * \param[in] prefix Prefix for bias parameters.
- * \param[in,out] wi Struct for bookeeping warnings.
- * \param[in] bComment True if comments should be printed.
- */
-void readBiasParams(std::vector<t_inpfile>* inp,
- AwhBiasParams* awhBiasParams,
- const std::string& prefix,
- warninp_t wi,
- bool bComment)
-{
- if (bComment)
- {
- printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
- }
-
- std::string opt = prefix + "-error-init";
- awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp,
- "Growth rate of the reference histogram determining the bias update "
- "size: exp-linear or linear");
- }
- opt = prefix + "-growth";
- awhBiasParams->eGrowth = getEnum<AwhHistogramGrowthType>(inp, opt.c_str(), wi);
-
- if (bComment)
- {
- printStringNoNewline(inp,
- "Start the simulation by equilibrating histogram towards the target "
- "distribution: no or yes");
- }
- opt = prefix + "-equilibrate-histogram";
- awhBiasParams->equilibrateHistogram = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
-
- if (bComment)
- {
- printStringNoNewline(
- inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
- }
- opt = prefix + "-target";
- awhBiasParams->eTarget = getEnum<AwhTargetType>(inp, opt.c_str(), wi);
-
- if (bComment)
- {
- printStringNoNewline(inp,
- "Boltzmann beta scaling factor for target distribution types "
- "'boltzmann' and 'boltzmann-local'");
- }
- opt = prefix + "-target-beta-scaling";
- awhBiasParams->targetBetaScaling = get_ereal(inp, opt, 0, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
- }
- opt = prefix + "-target-cutoff";
- awhBiasParams->targetCutoff = get_ereal(inp, opt, 0, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
- }
- opt = prefix + "-user-data";
- awhBiasParams->bUserData = getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No;
-
- if (bComment)
- {
- printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
- }
- opt = prefix + "-share-group";
- awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "Dimensionality of the coordinate");
- }
- opt = prefix + "-ndim";
- awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
-
- /* Check this before starting to read the AWH dimension parameters. */
- if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
- {
- gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), awhBiasParams->ndim, c_biasMaxNumDim);
- }
- snew(awhBiasParams->dimParams, awhBiasParams->ndim);
- for (int d = 0; d < awhBiasParams->ndim; d++)
- {
- bComment = bComment && d == 0;
- std::string prefixdim = prefix + formatString("-dim%d", d + 1);
- readDimParams(inp, prefixdim, &awhBiasParams->dimParams[d], wi, bComment);
- }
-}
-
/*! \brief
* Check the parameters of an AWH bias.
*
* \param[in] ir Input parameter struct.
* \param[in,out] wi Struct for bookeeping warnings.
*/
-void checkBiasParams(const AwhBiasParams* awhBiasParams, const std::string& prefix, const t_inputrec* ir, warninp_t wi)
+void checkBiasParams(const AwhBiasParams& awhBiasParams, const std::string& prefix, const t_inputrec& ir, warninp_t wi)
{
std::string opt = prefix + "-error-init";
- if (awhBiasParams->errorInitial <= 0)
+ if (awhBiasParams.initialErrorEstimate() <= 0)
{
gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
}
opt = prefix + "-equilibrate-histogram";
- if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != AwhHistogramGrowthType::ExponentialLinear)
+ if (awhBiasParams.equilibrateHistogram()
+ && awhBiasParams.growthType() != AwhHistogramGrowthType::ExponentialLinear)
{
auto message =
formatString("Option %s will only have an effect for histogram growth type '%s'.",
warning(wi, message);
}
- if ((awhBiasParams->eTarget == AwhTargetType::LocalBoltzmann)
- && (awhBiasParams->eGrowth == AwhHistogramGrowthType::ExponentialLinear))
+ if ((awhBiasParams.targetDistribution() == AwhTargetType::LocalBoltzmann)
+ && (awhBiasParams.growthType() == AwhHistogramGrowthType::ExponentialLinear))
{
auto message = formatString(
"Target type '%s' combined with histogram growth type '%s' is not "
}
opt = prefix + "-target-beta-scaling";
- switch (awhBiasParams->eTarget)
+ switch (awhBiasParams.targetDistribution())
{
case AwhTargetType::Boltzmann:
case AwhTargetType::LocalBoltzmann:
- if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
+ if (awhBiasParams.targetBetaScaling() < 0 || awhBiasParams.targetBetaScaling() > 1)
{
gmx_fatal(FARGS,
"%s = %g is not useful for target type %s.",
opt.c_str(),
- awhBiasParams->targetBetaScaling,
- enumValueToString(awhBiasParams->eTarget));
+ awhBiasParams.targetBetaScaling(),
+ enumValueToString(awhBiasParams.targetDistribution()));
}
break;
default:
- if (awhBiasParams->targetBetaScaling != 0)
+ if (awhBiasParams.targetBetaScaling() != 0)
{
gmx_fatal(
FARGS,
"Value for %s (%g) set explicitly but will not be used for target type %s.",
opt.c_str(),
- awhBiasParams->targetBetaScaling,
- enumValueToString(awhBiasParams->eTarget));
+ awhBiasParams.targetBetaScaling(),
+ enumValueToString(awhBiasParams.targetDistribution()));
}
break;
}
opt = prefix + "-target-cutoff";
- switch (awhBiasParams->eTarget)
+ switch (awhBiasParams.targetDistribution())
{
case AwhTargetType::Cutoff:
- if (awhBiasParams->targetCutoff <= 0)
+ if (awhBiasParams.targetCutoff() <= 0)
{
gmx_fatal(FARGS,
"%s = %g is not useful for target type %s.",
opt.c_str(),
- awhBiasParams->targetCutoff,
- enumValueToString(awhBiasParams->eTarget));
+ awhBiasParams.targetCutoff(),
+ enumValueToString(awhBiasParams.targetDistribution()));
}
break;
default:
- if (awhBiasParams->targetCutoff != 0)
+ if (awhBiasParams.targetCutoff() != 0)
{
gmx_fatal(
FARGS,
"Value for %s (%g) set explicitly but will not be used for target type %s.",
opt.c_str(),
- awhBiasParams->targetCutoff,
- enumValueToString(awhBiasParams->eTarget));
+ awhBiasParams.targetCutoff(),
+ enumValueToString(awhBiasParams.targetDistribution()));
}
break;
}
opt = prefix + "-share-group";
- if (awhBiasParams->shareGroup < 0)
+ if (awhBiasParams.shareGroup() < 0)
{
warning_error(wi, "AWH bias share-group should be >= 0");
}
opt = prefix + "-ndim";
- if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
+ if (awhBiasParams.ndim() <= 0 || awhBiasParams.ndim() > c_biasMaxNumDim)
{
- gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), awhBiasParams->ndim, c_biasMaxNumDim);
+ gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), awhBiasParams.ndim(), c_biasMaxNumDim);
}
- if (awhBiasParams->ndim > 2)
+ if (awhBiasParams.ndim() > 2)
{
warning_note(wi,
"For awh-dim > 2 the estimate based on the diffusion and the initial error is "
"currently only a rough guideline."
" You should verify its usefulness for your system before production runs!");
}
- for (int d = 0; d < awhBiasParams->ndim; d++)
+ for (int d = 0; d < awhBiasParams.ndim(); d++)
{
std::string prefixdim = prefix + formatString("-dim%d", d + 1);
- checkDimParams(prefixdim, &awhBiasParams->dimParams[d], ir, wi);
+ checkDimParams(prefixdim, awhBiasParams.dimParams()[d], ir, wi);
}
+}
- /* Check consistencies here that cannot be checked at read time at a lower level. */
- checkInputConsistencyAwhBias(*awhBiasParams, wi);
+/*! \brief
+ * Check consistency of input at the AWH bias level.
+ *
+ * \param[in] awhBiasParams AWH bias parameters.
+ * \param[in,out] wi Struct for bookkeeping warnings.
+ */
+void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
+{
+ /* Covering diameter and sharing warning. */
+ auto awhBiasDimensionParams = awhBiasParams.dimParams();
+ for (const auto& dimensionParam : awhBiasDimensionParams)
+ {
+ double coverDiameter = dimensionParam.coverDiameter();
+ if (awhBiasParams.shareGroup() <= 0 && coverDiameter > 0)
+ {
+ warning(wi,
+ "The covering diameter is only relevant to set for bias sharing simulations.");
+ }
+ }
}
/*! \brief
* Check that we have a shared bias when requesting multisim sharing.
*/
bool haveSharedBias = false;
- for (int k1 = 0; k1 < awhParams.numBias; k1++)
+ auto awhBiasParams = awhParams.awhBiasParams();
+ for (int k1 = 0; k1 < awhParams.numBias(); k1++)
{
- const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1];
+ const AwhBiasParams& awhBiasParams1 = awhBiasParams[k1];
- if (awhBiasParams1.shareGroup > 0)
+ if (awhBiasParams1.shareGroup() > 0)
{
haveSharedBias = true;
}
/* k1 is the reference AWH, k2 is the AWH we compare with (can be equal to k1) */
- for (int k2 = k1; k2 < awhParams.numBias; k2++)
+ for (int k2 = k1; k2 < awhParams.numBias(); k2++)
{
- for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
+ const AwhBiasParams& awhBiasParams2 = awhBiasParams[k2];
+ const auto& dimParams1 = awhBiasParams1.dimParams();
+ const auto& dimParams2 = awhBiasParams2.dimParams();
+ for (int d1 = 0; d1 < gmx::ssize(dimParams1); d1++)
{
- if (awhBiasParams1.dimParams[d1].eCoordProvider == AwhCoordinateProviderType::FreeEnergyLambda)
+ if (dimParams1[d1].coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda)
{
continue;
}
- const AwhBiasParams& awhBiasParams2 = awhParams.awhBiasParams[k2];
-
/* d1 is the reference dimension of the reference AWH. d2 is the dim index of the AWH to compare with. */
- for (int d2 = 0; d2 < awhBiasParams2.ndim; d2++)
+ for (int d2 = 0; d2 < gmx::ssize(dimParams2); d2++)
{
- if (awhBiasParams2.dimParams[d2].eCoordProvider == AwhCoordinateProviderType::FreeEnergyLambda)
+ if (dimParams2[d2].coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda)
{
continue;
}
/* Give an error if (d1, k1) is different from (d2, k2) but the pull coordinate is the same */
if ((d1 != d2 || k1 != k2)
- && (awhBiasParams1.dimParams[d1].coordIndex == awhBiasParams2.dimParams[d2].coordIndex))
+ && (dimParams1[d1].coordinateIndex() == dimParams2[d2].coordinateIndex()))
{
char errormsg[STRLEN];
sprintf(errormsg,
"dimensions (awh%d-dim%d and awh%d-dim%d). "
"If this is really what you want to do you will have to duplicate "
"this pull coordinate.",
- awhBiasParams1.dimParams[d1].coordIndex + 1,
+ dimParams1[d1].coordinateIndex() + 1,
k1 + 1,
d1 + 1,
k2 + 1,
}
}
- if (awhParams.shareBiasMultisim && !haveSharedBias)
+ if (awhParams.shareBiasMultisim() && !haveSharedBias)
{
warning(wi,
"Sharing of biases over multiple simulations is requested, but no bias is marked "
"this (yet)");
}
}
+
} // namespace
-AwhParams* readAwhParams(std::vector<t_inpfile>* inp, warninp_t wi)
+AwhDimParams::AwhDimParams(std::vector<t_inpfile>* inp,
+ const std::string& prefix,
+ const t_inputrec& ir,
+ warninp_t wi,
+ bool bComment)
+{
+ std::string opt;
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp,
+ "The provider of the reaction coordinate, "
+ "currently only 'pull' and 'fep-lambda' (free energy lambda state) is supported");
+ }
+
+ opt = prefix + "-coord-provider";
+ eCoordProvider_ = getEnum<AwhCoordinateProviderType>(inp, opt.c_str(), wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "The coordinate index for this dimension");
+ }
+ opt = prefix + "-coord-index";
+ int coordIndexInput;
+ coordIndexInput = get_eint(inp, opt, 1, wi);
+ if (coordIndexInput < 1)
+ {
+ gmx_fatal(FARGS,
+ "Failed to read a valid coordinate index for %s. "
+ "Note that the pull coordinate indexing starts at 1.",
+ opt.c_str());
+ }
+
+ /* The pull coordinate indices start at 1 in the input file, at 0 internally */
+ coordIndex_ = coordIndexInput - 1;
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Start and end values for each coordinate dimension");
+ }
+
+ opt = prefix + "-start";
+ origin_ = get_ereal(inp, opt, 0., wi);
+
+ opt = prefix + "-end";
+ end_ = get_ereal(inp, opt, 0., wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
+ }
+ opt = prefix + "-force-constant";
+ forceConstant_ = get_ereal(inp, opt, 0, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps or ps^-1)");
+ }
+ opt = prefix + "-diffusion";
+ double diffusionValue = get_ereal(inp, opt, 0, wi);
+ if (diffusionValue <= 0)
+ {
+ const double diffusion_default = 1e-5;
+ auto message = formatString(
+ "%s not explicitly set by user. You can choose to use a default "
+ "value (%g nm^2/ps or rad^2/ps) but this may very well be "
+ "non-optimal for your system!",
+ opt.c_str(),
+ diffusion_default);
+ warning(wi, message);
+ diffusionValue = diffusion_default;
+ }
+ diffusion_ = diffusionValue;
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Diameter that needs to be sampled around a point before it is "
+ "considered covered. In FEP dimensions the cover diameter is "
+ "specified in lambda states.");
+ }
+ opt = prefix + "-cover-diameter";
+ coverDiameter_ = get_ereal(inp, opt, 0, wi);
+
+ checkDimParams(prefix, *this, ir, wi);
+}
+
+AwhDimParams::AwhDimParams(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(serializer->reading(),
+ "Can not use writing serializer for creating datastructure");
+ serializer->doEnumAsInt(&eCoordProvider_);
+ serializer->doInt(&coordIndex_);
+ serializer->doDouble(&origin_);
+ serializer->doDouble(&end_);
+ serializer->doDouble(&period_);
+ serializer->doDouble(&forceConstant_);
+ serializer->doDouble(&diffusion_);
+ serializer->doDouble(&coordValueInit_);
+ serializer->doDouble(&coverDiameter_);
+}
+
+void AwhDimParams::serialize(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(!serializer->reading(),
+ "Can not use reading serializer for writing datastructure");
+ serializer->doEnumAsInt(&eCoordProvider_);
+ serializer->doInt(&coordIndex_);
+ serializer->doDouble(&origin_);
+ serializer->doDouble(&end_);
+ serializer->doDouble(&period_);
+ serializer->doDouble(&forceConstant_);
+ serializer->doDouble(&diffusion_);
+ serializer->doDouble(&coordValueInit_);
+ serializer->doDouble(&coverDiameter_);
+}
+
+AwhBiasParams::AwhBiasParams(std::vector<t_inpfile>* inp,
+ const std::string& prefix,
+ const t_inputrec& ir,
+ warninp_t wi,
+ bool bComment)
+{
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
+ }
+
+ std::string opt = prefix + "-error-init";
+ errorInitial_ = get_ereal(inp, opt, 10, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Growth rate of the reference histogram determining the bias update "
+ "size: exp-linear or linear");
+ }
+ opt = prefix + "-growth";
+ eGrowth_ = getEnum<AwhHistogramGrowthType>(inp, opt.c_str(), wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Start the simulation by equilibrating histogram towards the target "
+ "distribution: no or yes");
+ }
+ opt = prefix + "-equilibrate-histogram";
+ equilibrateHistogram_ = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
+
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
+ }
+ opt = prefix + "-target";
+ eTarget_ = getEnum<AwhTargetType>(inp, opt.c_str(), wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Boltzmann beta scaling factor for target distribution types "
+ "'boltzmann' and 'boltzmann-local'");
+ }
+ opt = prefix + "-target-beta-scaling";
+ targetBetaScaling_ = get_ereal(inp, opt, 0, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
+ }
+ opt = prefix + "-target-cutoff";
+ targetCutoff_ = get_ereal(inp, opt, 0, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
+ }
+ opt = prefix + "-user-data";
+ bUserData_ = getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No;
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
+ }
+ opt = prefix + "-share-group";
+ shareGroup_ = get_eint(inp, opt, 0, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Dimensionality of the coordinate");
+ }
+ opt = prefix + "-ndim";
+ int ndim = get_eint(inp, opt, 0, wi);
+
+ /* Check this before starting to read the AWH dimension parameters. */
+ if (ndim <= 0 || ndim > c_biasMaxNumDim)
+ {
+ gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), ndim, c_biasMaxNumDim);
+ }
+ for (int d = 0; d < ndim; d++)
+ {
+ bComment = bComment && d == 0;
+ std::string prefixdim = prefix + formatString("-dim%d", d + 1);
+ dimParams_.emplace_back(inp, prefixdim, ir, wi, bComment);
+ }
+ checkBiasParams(*this, prefix, ir, wi);
+}
+
+AwhBiasParams::AwhBiasParams(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(serializer->reading(),
+ "Can not use writing serializer to create datastructure");
+ serializer->doEnumAsInt(&eTarget_);
+ serializer->doDouble(&targetBetaScaling_);
+ serializer->doDouble(&targetCutoff_);
+ serializer->doEnumAsInt(&eGrowth_);
+ int temp = 0;
+ serializer->doInt(&temp);
+ bUserData_ = static_cast<bool>(temp);
+ serializer->doDouble(&errorInitial_);
+ int numDimensions = dimParams_.size();
+ serializer->doInt(&numDimensions);
+ serializer->doInt(&shareGroup_);
+ serializer->doBool(&equilibrateHistogram_);
+
+ for (int k = 0; k < numDimensions; k++)
+ {
+ dimParams_.emplace_back(serializer);
+ }
+ /* Check consistencies here that cannot be checked at read time at a lower level. */
+ checkInputConsistencyAwhBias(*this, nullptr);
+}
+
+void AwhBiasParams::serialize(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(!serializer->reading(),
+ "Can not use reading serializer to write datastructure");
+ serializer->doEnumAsInt(&eTarget_);
+ serializer->doDouble(&targetBetaScaling_);
+ serializer->doDouble(&targetCutoff_);
+ serializer->doEnumAsInt(&eGrowth_);
+ int temp = static_cast<int>(bUserData_);
+ serializer->doInt(&temp);
+ serializer->doDouble(&errorInitial_);
+ int numDimensions = ndim();
+ serializer->doInt(&numDimensions);
+ serializer->doInt(&shareGroup_);
+ serializer->doBool(&equilibrateHistogram_);
+
+ for (int k = 0; k < numDimensions; k++)
+ {
+ dimParams_[k].serialize(serializer);
+ }
+}
+
+AwhParams::AwhParams(std::vector<t_inpfile>* inp, const t_inputrec& ir, warninp_t wi)
{
- AwhParams* awhParams;
- snew(awhParams, 1);
std::string opt;
/* Parameters common for all biases */
printStringNoNewline(inp, "The way to apply the biasing potential: convolved or umbrella");
- opt = "awh-potential";
-
awhParams->ePotential = getEnum<AwhPotentialType>(inp, opt.c_str(), wi);
+ opt = "awh-potential";
+
potentialEnum_ = getEnum<AwhPotentialType>(inp, opt.c_str(), wi);
printStringNoNewline(inp,
"The random seed used for sampling the umbrella center in the case of "
"umbrella type potential");
- opt = "awh-seed";
- awhParams->seed = get_eint(inp, opt, -1, wi);
- if (awhParams->seed == -1)
+ opt = "awh-seed";
+ seed_ = get_eint(inp, opt, -1, wi);
+ if (seed_ == -1)
{
- awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
- fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
+ seed_ = static_cast<int>(gmx::makeRandomSeed());
+ fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", seed_);
}
printStringNoNewline(inp, "Data output interval in number of steps");
- opt = "awh-nstout";
- awhParams->nstOut = get_eint(inp, opt, 100000, wi);
+ opt = "awh-nstout";
+ nstOut_ = get_eint(inp, opt, 100000, wi);
printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
- opt = "awh-nstsample";
- awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
+ opt = "awh-nstsample";
+ nstSampleCoord_ = get_eint(inp, opt, 10, wi);
printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
- opt = "awh-nsamples-update";
- awhParams->numSamplesUpdateFreeEnergy = get_eint(inp, opt, 10, wi);
+ opt = "awh-nsamples-update";
+ numSamplesUpdateFreeEnergy_ = get_eint(inp, opt, 10, wi);
printStringNoNewline(
inp, "When true, biases with share-group>0 are shared between multiple simulations");
- opt = "awh-share-multisim";
- awhParams->shareBiasMultisim = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
+ opt = "awh-share-multisim";
+ shareBiasMultisim_ = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
printStringNoNewline(inp, "The number of independent AWH biases");
- opt = "awh-nbias";
- awhParams->numBias = get_eint(inp, opt, 1, wi);
+ opt = "awh-nbias";
+ int numBias = get_eint(inp, opt, 1, wi);
/* Check this before starting to read the AWH biases. */
- if (awhParams->numBias <= 0)
+ if (numBias <= 0)
{
gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
}
/* Read the parameters specific to each AWH bias */
- snew(awhParams->awhBiasParams, awhParams->numBias);
-
- for (int k = 0; k < awhParams->numBias; k++)
+ for (int k = 0; k < numBias; k++)
{
bool bComment = (k == 0);
std::string prefixawh = formatString("awh%d", k + 1);
- readBiasParams(inp, &awhParams->awhBiasParams[k], prefixawh, wi, bComment);
+ awhBiasParams_.emplace_back(inp, prefixawh, ir, wi, bComment);
}
-
- return awhParams;
+ checkAwhParams(*this, ir, wi);
+ checkInputConsistencyAwh(*this, wi);
}
-void checkAwhParams(const AwhParams* awhParams, const t_inputrec* ir, warninp_t wi)
+AwhParams::AwhParams(ISerializer* serializer)
{
- std::string opt;
-
- checkMtsConsistency(*ir, wi);
-
- opt = "awh-nstout";
- if (awhParams->nstOut <= 0)
- {
- auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
- opt.c_str(),
- awhParams->nstOut);
- warning_error(wi, message);
- }
- /* This restriction can be removed by changing a flag of print_ebin() */
- if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
- {
- auto message = formatString(
- "%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(), awhParams->nstOut, ir->nstenergy);
- warning_error(wi, message);
- }
-
- opt = "awh-nsamples-update";
- if (awhParams->numSamplesUpdateFreeEnergy <= 0)
- {
- warning_error(wi, opt + " needs to be an integer > 0");
- }
-
- for (int k = 0; k < awhParams->numBias; k++)
- {
- std::string prefixawh = formatString("awh%d", k + 1);
- checkBiasParams(&awhParams->awhBiasParams[k], prefixawh, ir, wi);
+ GMX_RELEASE_ASSERT(serializer->reading(),
+ "Can not use writing serializer to read AWH parameters");
+ int numberOfBiases = awhBiasParams_.size();
+ serializer->doInt(&numberOfBiases);
+ serializer->doInt(&nstOut_);
+ serializer->doInt64(&seed_);
+ serializer->doInt(&nstSampleCoord_);
+ serializer->doInt(&numSamplesUpdateFreeEnergy_);
+ serializer->doEnumAsInt(&potentialEnum_);
+ serializer->doBool(&shareBiasMultisim_);
+
+ if (numberOfBiases > 0)
+ {
+ for (int k = 0; k < numberOfBiases; k++)
+ {
+ awhBiasParams_.emplace_back(serializer);
+ }
}
+ checkInputConsistencyAwh(*this, nullptr);
+}
- /* Do a final consistency check before returning */
- checkInputConsistencyAwh(*awhParams, wi);
-
- if (ir->init_step != 0)
- {
- warning_error(wi, "With AWH init-step should be 0");
+void AwhParams::serialize(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(!serializer->reading(),
+ "Can not use reading serializer to write AWH parameters");
+ int numberOfBiases = numBias();
+ serializer->doInt(&numberOfBiases);
+ serializer->doInt(&nstOut_);
+ serializer->doInt64(&seed_);
+ serializer->doInt(&nstSampleCoord_);
+ serializer->doInt(&numSamplesUpdateFreeEnergy_);
+ serializer->doEnumAsInt(&potentialEnum_);
+ serializer->doBool(&shareBiasMultisim_);
+
+ if (numberOfBiases > 0)
+ {
+ for (int k = 0; k < numberOfBiases; k++)
+ {
+ awhBiasParams_[k].serialize(serializer);
+ }
}
}
* \param[in] awhParams AWH parameters.
* \param[in,out] wi Struct for bookeeping warnings.
*/
-static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t wi)
+static void checkInputConsistencyInterval(const AwhParams& awhParams, warninp_t wi)
{
- for (int k = 0; k < awhParams->numBias; k++)
+ const auto& awhBiasParams = awhParams.awhBiasParams();
+ for (int k = 0; k < gmx::ssize(awhBiasParams); k++)
{
- AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
- for (int d = 0; d < awhBiasParams->ndim; d++)
+ const auto& dimParams = awhBiasParams[k].dimParams();
+ for (int d = 0; d < gmx::ssize(dimParams); d++)
{
- AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
- int coordIndex = dimParams->coordIndex;
- double origin = dimParams->origin, end = dimParams->end, period = dimParams->period;
- double coordValueInit = dimParams->coordValueInit;
+ int coordIndex = dimParams[d].coordinateIndex();
+ double origin = dimParams[d].origin(), end = dimParams[d].end(),
+ period = dimParams[d].period();
+ double coordValueInit = dimParams[d].initialCoordinate();
if ((period == 0) && (origin > end))
{
static void setStateDependentAwhPullDimParams(AwhDimParams* dimParams,
const int biasIndex,
const int dimIndex,
- const pull_params_t* pull_params,
+ const pull_params_t& pull_params,
pull_t* pull_work,
const t_pbc& pbc,
const tensor& compressibility,
warninp_t wi)
{
- const t_pull_coord& pullCoordParams = pull_params->coord[dimParams->coordIndex];
+ const t_pull_coord& pullCoordParams = pull_params.coord[dimParams->coordinateIndex()];
if (pullCoordParams.eGeom == PullGroupGeometry::DirectionPBC)
{
enumValueToString(PullGroupGeometry::Direction));
}
- dimParams->period = get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin);
+ dimParams->setPeriod(
+ get_pull_coord_period(pullCoordParams, pbc, dimParams->end() - dimParams->origin()));
// We would like to check for scaling, but we don't have the full inputrec available here
- if (dimParams->period > 0
+ if (dimParams->period() > 0
&& !(pullCoordParams.eGeom == PullGroupGeometry::Angle
|| pullCoordParams.eGeom == PullGroupGeometry::Dihedral))
{
{
std::string mesg = gmx::formatString(
"AWH dimension %d of bias %d is periodic with pull geometry '%s', "
- "while you should are applying pressure scaling to the "
+ "while you should be applying pressure scaling to the "
"corresponding box vector, this is not supported.",
dimIndex + 1,
biasIndex + 1,
}
/* The initial coordinate value, converted to external user units. */
- dimParams->coordValueInit = get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
-
- dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(&pullCoordParams);
+ double initialCoordinate = get_pull_coord_value(pull_work, dimParams->coordinateIndex(), &pbc);
+ initialCoordinate *= pull_conversion_factor_internal2userinput(&pullCoordParams);
+ dimParams->setInitialCoordinate(initialCoordinate);
}
void setStateDependentAwhParams(AwhParams* awhParams,
t_pbc pbc;
set_pbc(&pbc, pbcType, box);
- for (int k = 0; k < awhParams->numBias; k++)
+ auto awhBiasParams = awhParams->awhBiasParams();
+ for (int k = 0; k < awhParams->numBias(); k++)
{
- AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
- for (int d = 0; d < awhBiasParams->ndim; d++)
+ auto awhBiasDimensionParams = awhBiasParams[k].dimParams();
+ for (int d = 0; d < awhBiasParams[k].ndim(); d++)
{
- AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
- if (dimParams->eCoordProvider == AwhCoordinateProviderType::Pull)
+ AwhDimParams* dimParams = &awhBiasDimensionParams[d];
+ if (dimParams->coordinateProvider() == AwhCoordinateProviderType::Pull)
{
setStateDependentAwhPullDimParams(
- dimParams, k, d, &pull_params, pull_work, pbc, compressibility, wi);
+ dimParams, k, d, pull_params, pull_work, pbc, compressibility, wi);
}
else
{
- dimParams->coordValueInit = initLambda;
+ dimParams->setInitialCoordinate(initLambda);
checkFepLambdaDimDecouplingConsistency(mtop, wi);
}
}
}
- checkInputConsistencyInterval(awhParams, wi);
+ checkInputConsistencyInterval(*awhParams, wi);
/* Register AWH as external potential with pull (for AWH dimensions that use pull as
* reaction coordinate provider) to check consistency. */
Awh::registerAwhWithPull(*awhParams, pull_work);
}
+void checkAwhParams(const AwhParams& awhParams, const t_inputrec& ir, warninp_t wi)
+{
+ std::string opt;
+ checkMtsConsistency(ir, wi);
+
+ opt = "awh-nstout";
+ if (awhParams.nstout() <= 0)
+ {
+ auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
+ opt.c_str(),
+ awhParams.nstout());
+ warning_error(wi, message);
+ }
+ /* This restriction can be removed by changing a flag of print_ebin() */
+ if (ir.nstenergy == 0 || awhParams.nstout() % ir.nstenergy != 0)
+ {
+ auto message = formatString(
+ "%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(), awhParams.nstout(), ir.nstenergy);
+ warning_error(wi, message);
+ }
+
+ opt = "awh-nsamples-update";
+ if (awhParams.numSamplesUpdateFreeEnergy() <= 0)
+ {
+ warning_error(wi, opt + " needs to be an integer > 0");
+ }
+
+ for (int k = 0; k < awhParams.numBias(); k++)
+ {
+ std::string prefixawh = formatString("awh%d", k + 1);
+ checkBiasParams(awhParams.awhBiasParams()[k], prefixawh, ir, wi);
+ }
+
+ if (ir.init_step != 0)
+ {
+ warning_error(wi, "With AWH init-step should be 0");
+ }
+}
+
} // namespace gmx
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
namespace gmx
{
-struct AwhParams;
-/*! \brief Allocate and initialize the AWH parameters with values from the input file.
- *
- * \param[in,out] inp Input file entries.
- * \param[in,out] wi Struct for bookeeping warnings.
- * \returns AWH parameters.
- */
-AwhParams* readAwhParams(std::vector<t_inpfile>* inp, warninp_t wi);
+class AwhParams;
+class ISerializer;
/*! \brief Check the AWH parameters.
*
- * \param[in,out] awhParams The AWH parameters.
+ * \param[in] awhParams The AWH parameters.
* \param[in] inputrec Input parameter struct.
* \param[in,out] wi Struct for bookeeping warnings.
*/
-void checkAwhParams(const AwhParams* awhParams, const t_inputrec* inputrec, warninp_t wi);
+void checkAwhParams(const AwhParams& awhParams, const t_inputrec& inputrec, warninp_t wi);
/*! \brief
#
# This file is part of the GROMACS molecular simulation package.
#
-# Copyright (c) 2017,2020, by the GROMACS development team, led by
+# Copyright (c) 2017,2020,2021, by the GROMACS development team, led by
# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
# and including many others, as listed in the AUTHORS file in the
# top-level source directory and at http://www.gromacs.org.
gmx_add_unit_test(AwhTest awh-test
CPP_SOURCE_FILES
+ awh_setup.cpp
bias.cpp
biasgrid.cpp
biasstate.cpp
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * 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 "awh_setup.h"
+
+#include <cmath>
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+#include <gmock/gmock.h>
+#include <gmock/gmock-matchers.h>
+#include <gtest/gtest.h>
+
+#include "gromacs/applied_forces/awh/bias.h"
+#include "gromacs/applied_forces/awh/correlationgrid.h"
+#include "gromacs/applied_forces/awh/pointstate.h"
+#include "gromacs/mdtypes/awh_params.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/inmemoryserializer.h"
+#include "gromacs/utility/stringutil.h"
+
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+
+namespace gmx
+{
+
+namespace test
+{
+
+using ::testing::Eq;
+using ::testing::Pointwise;
+
+std::vector<char> awhDimParamSerialized(AwhCoordinateProviderType inputCoordinateProvider,
+ int inputCoordIndex,
+ double inputOrigin,
+ double inputEnd,
+ double inputPeriod,
+ double inputDiffusion)
+{
+ AwhCoordinateProviderType eCoordProvider = inputCoordinateProvider;
+ int coordIndex = inputCoordIndex;
+ double forceConstant = 10;
+ double period = inputPeriod;
+ double diffusion = inputDiffusion;
+ double origin = inputOrigin;
+ double end = inputEnd;
+ double coordValueInit = inputOrigin;
+ double coverDiameter = 0;
+
+ gmx::InMemorySerializer serializer;
+ serializer.doEnumAsInt(&eCoordProvider);
+ serializer.doInt(&coordIndex);
+ serializer.doDouble(&origin);
+ serializer.doDouble(&end);
+ serializer.doDouble(&period);
+ serializer.doDouble(&forceConstant);
+ serializer.doDouble(&diffusion);
+ serializer.doDouble(&coordValueInit);
+ serializer.doDouble(&coverDiameter);
+ return serializer.finishAndGetBuffer();
+}
+
+/*! \internal \brief
+ * Prepare a memory buffer with serialized AwhBiasParams.
+ *
+ * \param[in] eawhgrowth Way to grow potential.
+ * \param[in] beta Value for 1/(kB*T).
+ * \param[in] inputErrorScaling Factor for initial error scaling.
+ * \param[in] dimensionParameterBuffers Buffers containing the dimension parameters.
+ * \param[in] inputUserData If there is a user provided PMF estimate.
+ */
+static std::vector<char> awhBiasParamSerialized(AwhHistogramGrowthType eawhgrowth,
+ double beta,
+ double inputErrorScaling,
+ ArrayRef<const std::vector<char>> dimensionParameterBuffers,
+ bool inputUserData)
+{
+ int ndim = dimensionParameterBuffers.size();
+ AwhTargetType eTarget = AwhTargetType::Constant;
+ double targetBetaScaling = 0;
+ double targetCutoff = 0;
+ AwhHistogramGrowthType eGrowth = eawhgrowth;
+ bool bUserData = inputUserData;
+ double errorInitial = inputErrorScaling / beta;
+ int shareGroup = 0;
+ bool equilibrateHistogram = false;
+
+ gmx::InMemorySerializer serializer;
+ serializer.doEnumAsInt(&eTarget);
+ serializer.doDouble(&targetBetaScaling);
+ serializer.doDouble(&targetCutoff);
+ serializer.doEnumAsInt(&eGrowth);
+ int temp = static_cast<int>(bUserData);
+ serializer.doInt(&temp);
+ serializer.doDouble(&errorInitial);
+ serializer.doInt(&ndim);
+ serializer.doInt(&shareGroup);
+ serializer.doBool(&equilibrateHistogram);
+
+ auto awhDimBuffer = awhDimParamSerialized();
+ auto awhBiasBuffer = serializer.finishAndGetBuffer();
+ for (const auto& dimParamBuffer : dimensionParameterBuffers)
+ {
+ awhBiasBuffer.insert(awhBiasBuffer.end(), dimParamBuffer.begin(), dimParamBuffer.end());
+ }
+ return awhBiasBuffer;
+}
+
+/*! \internal \brief
+ * Prepare a memory buffer with serialized AwhParams.
+ *
+ * \param[in] eawhgrowth Way to grow potential.
+ * \param[in] eawhpotential Which potential to use.
+ * \param[in] beta Value for 1/(kB*T).
+ * \param[in] inputErrorScaling Factor for initial error scaling.
+ * \param[in] inputSeed Seed value to use.
+ * \param[in] dimensionParameterBuffers Buffers containing the dimension parameters.
+ * \param[in] inputUserData If there is a user provided PMF estimate.
+ */
+static std::vector<char> awhParamSerialized(AwhHistogramGrowthType eawhgrowth,
+ AwhPotentialType eawhpotential,
+ double beta,
+ double inputErrorScaling,
+ int64_t inputSeed,
+ ArrayRef<const std::vector<char>> dimensionParameterBuffers,
+ bool inputUserData)
+{
+ int numBias = 1;
+ int64_t seed = inputSeed;
+ int nstOut = 0;
+ int nstSampleCoord = 1;
+ int numSamplesUpdateFreeEnergy = 10;
+ AwhPotentialType ePotential = eawhpotential;
+ bool shareBiasMultisim = false;
+
+ gmx::InMemorySerializer serializer;
+ serializer.doInt(&numBias);
+ serializer.doInt(&nstOut);
+ serializer.doInt64(&seed);
+ serializer.doInt(&nstSampleCoord);
+ serializer.doInt(&numSamplesUpdateFreeEnergy);
+ serializer.doEnumAsInt(&ePotential);
+ serializer.doBool(&shareBiasMultisim);
+
+ auto awhParamBuffer = serializer.finishAndGetBuffer();
+ auto awhBiasBuffer = awhBiasParamSerialized(
+ eawhgrowth, beta, inputErrorScaling, dimensionParameterBuffers, inputUserData);
+
+ awhParamBuffer.insert(awhParamBuffer.end(), awhBiasBuffer.begin(), awhBiasBuffer.end());
+
+ return awhParamBuffer;
+}
+
+AwhTestParameters::AwhTestParameters(ISerializer* serializer) : awhParams(serializer) {}
+/*! \brief
+ * Helper function to set up the C-style AWH parameters for the test.
+ *
+ * Builds the test input data from serialized data.
+ */
+AwhTestParameters getAwhTestParameters(AwhHistogramGrowthType eawhgrowth,
+ AwhPotentialType eawhpotential,
+ ArrayRef<const std::vector<char>> dimensionParameterBuffers,
+ bool inputUserData,
+ double beta,
+ bool useAwhFep,
+ double inputErrorScaling,
+ int numFepLambdaStates)
+{
+ double convFactor = 1;
+ double k = 1000;
+ int64_t seed = 93471803;
+
+ auto awhParamBuffer = awhParamSerialized(
+ eawhgrowth, eawhpotential, beta, inputErrorScaling, seed, dimensionParameterBuffers, inputUserData);
+ gmx::InMemoryDeserializer deserializer(awhParamBuffer, false);
+ AwhTestParameters params(&deserializer);
+
+ params.beta = beta;
+
+ if (useAwhFep)
+ {
+ params.dimParams.emplace_back(DimParams::fepLambdaDimParams(numFepLambdaStates, params.beta));
+ }
+ else
+ {
+ params.dimParams.emplace_back(DimParams::pullDimParams(convFactor, k, params.beta));
+ }
+ return params;
+}
+
+TEST(SerializationTest, CanSerializeDimParams)
+{
+ auto awhDimBuffer = awhDimParamSerialized();
+ gmx::InMemoryDeserializer deserializer(awhDimBuffer, false);
+ AwhDimParams awhDimParams(&deserializer);
+ EXPECT_EQ(awhDimParams.coordinateProvider(), AwhCoordinateProviderType::Pull);
+ EXPECT_EQ(awhDimParams.coordinateIndex(), 0);
+ EXPECT_FLOAT_EQ(awhDimParams.forceConstant(), 10);
+ EXPECT_FLOAT_EQ(awhDimParams.period(), 0);
+ EXPECT_FLOAT_EQ(awhDimParams.diffusion(), 0.34690997);
+ EXPECT_FLOAT_EQ(awhDimParams.origin(), 0.5);
+ EXPECT_FLOAT_EQ(awhDimParams.end(), 1.5);
+ EXPECT_FLOAT_EQ(awhDimParams.initialCoordinate(), awhDimParams.origin());
+ EXPECT_FLOAT_EQ(awhDimParams.coverDiameter(), 0);
+
+ gmx::InMemorySerializer serializer;
+ awhDimParams.serialize(&serializer);
+ EXPECT_THAT(awhDimBuffer, Pointwise(Eq(), serializer.finishAndGetBuffer()));
+}
+
+TEST(SerializationTest, CanSerializeBiasParams)
+{
+ auto awhDimBuffer = awhDimParamSerialized();
+ auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
+ auto awhBiasBuffer = awhBiasParamSerialized(
+ AwhHistogramGrowthType::ExponentialLinear, 0.4, 0.5, awhDimArrayRef, false);
+ gmx::InMemoryDeserializer deserializer(awhBiasBuffer, false);
+ AwhBiasParams awhBiasParams(&deserializer);
+ EXPECT_EQ(awhBiasParams.ndim(), 1);
+ EXPECT_EQ(awhBiasParams.targetDistribution(), AwhTargetType::Constant);
+ EXPECT_FLOAT_EQ(awhBiasParams.targetBetaScaling(), 0);
+ EXPECT_FLOAT_EQ(awhBiasParams.targetCutoff(), 0);
+ EXPECT_EQ(awhBiasParams.growthType(), AwhHistogramGrowthType::ExponentialLinear);
+ EXPECT_EQ(awhBiasParams.userPMFEstimate(), 0);
+ EXPECT_FLOAT_EQ(awhBiasParams.initialErrorEstimate(), 0.5 / 0.4);
+ EXPECT_EQ(awhBiasParams.shareGroup(), 0);
+ EXPECT_EQ(awhBiasParams.equilibrateHistogram(), false);
+ const auto& awhDimParams = awhBiasParams.dimParams()[0];
+ EXPECT_EQ(awhDimParams.coordinateProvider(), AwhCoordinateProviderType::Pull);
+ EXPECT_EQ(awhDimParams.coordinateIndex(), 0);
+ EXPECT_FLOAT_EQ(awhDimParams.forceConstant(), 10);
+ EXPECT_FLOAT_EQ(awhDimParams.period(), 0);
+ EXPECT_FLOAT_EQ(awhDimParams.diffusion(), 0.34690997);
+ EXPECT_FLOAT_EQ(awhDimParams.origin(), 0.5);
+ EXPECT_FLOAT_EQ(awhDimParams.end(), 1.5);
+ EXPECT_FLOAT_EQ(awhDimParams.initialCoordinate(), awhDimParams.origin());
+ EXPECT_FLOAT_EQ(awhDimParams.coverDiameter(), 0);
+
+ gmx::InMemorySerializer serializer;
+ awhBiasParams.serialize(&serializer);
+ EXPECT_THAT(awhBiasBuffer, Pointwise(Eq(), serializer.finishAndGetBuffer()));
+}
+
+TEST(SerializationTest, CanSerializeAwhParams)
+{
+ auto awhDimBuffer = awhDimParamSerialized();
+ auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
+ auto awhParamBuffer = awhParamSerialized(
+ AwhHistogramGrowthType::ExponentialLinear, AwhPotentialType::Convolved, 0.4, 0.5, 1337, awhDimArrayRef, false);
+ gmx::InMemoryDeserializer deserializer(awhParamBuffer, false);
+ AwhParams awhParams(&deserializer);
+ EXPECT_EQ(awhParams.numBias(), 1);
+ EXPECT_EQ(awhParams.seed(), 1337);
+ EXPECT_EQ(awhParams.nstout(), 0);
+ EXPECT_EQ(awhParams.nstSampleCoord(), 1);
+ EXPECT_EQ(awhParams.numSamplesUpdateFreeEnergy(), 10);
+ EXPECT_EQ(awhParams.potential(), AwhPotentialType::Convolved);
+ EXPECT_EQ(awhParams.shareBiasMultisim(), false);
+ const auto& awhBiasParams = awhParams.awhBiasParams()[0];
+ EXPECT_EQ(awhBiasParams.ndim(), 1);
+ EXPECT_EQ(awhBiasParams.targetDistribution(), AwhTargetType::Constant);
+ EXPECT_FLOAT_EQ(awhBiasParams.targetBetaScaling(), 0);
+ EXPECT_FLOAT_EQ(awhBiasParams.targetCutoff(), 0);
+ EXPECT_EQ(awhBiasParams.growthType(), AwhHistogramGrowthType::ExponentialLinear);
+ EXPECT_EQ(awhBiasParams.userPMFEstimate(), 0);
+ EXPECT_FLOAT_EQ(awhBiasParams.initialErrorEstimate(), 0.5 / 0.4);
+ EXPECT_EQ(awhBiasParams.shareGroup(), 0);
+ EXPECT_EQ(awhBiasParams.equilibrateHistogram(), false);
+ const auto& awhDimParams = awhBiasParams.dimParams()[0];
+ EXPECT_EQ(awhDimParams.coordinateProvider(), AwhCoordinateProviderType::Pull);
+ EXPECT_EQ(awhDimParams.coordinateIndex(), 0);
+ EXPECT_FLOAT_EQ(awhDimParams.forceConstant(), 10);
+ EXPECT_FLOAT_EQ(awhDimParams.period(), 0);
+ EXPECT_FLOAT_EQ(awhDimParams.diffusion(), 0.34690997);
+ EXPECT_FLOAT_EQ(awhDimParams.origin(), 0.5);
+ EXPECT_FLOAT_EQ(awhDimParams.end(), 1.5);
+ EXPECT_FLOAT_EQ(awhDimParams.initialCoordinate(), awhDimParams.origin());
+ EXPECT_FLOAT_EQ(awhDimParams.coverDiameter(), 0);
+
+ gmx::InMemorySerializer serializer;
+ awhParams.serialize(&serializer);
+ EXPECT_THAT(awhParamBuffer, Pointwise(Eq(), serializer.finishAndGetBuffer()));
+}
+
+} // namespace test
+} // namespace gmx
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * 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.
+ */
+#ifndef GMX_AWH_TEST_SETUP_H
+#define GMX_AWH_TEST_SETUP_H
+
+#include "gmxpre.h"
+
+#include <memory>
+#include <vector>
+
+#include "gromacs/applied_forces/awh/bias.h"
+#include "gromacs/mdtypes/awh_params.h"
+
+namespace gmx
+{
+
+template<typename>
+class ArrayRef;
+
+namespace test
+{
+
+/*! \internal \brief
+ * Prepare a memory buffer with serialized AwhDimParams.
+ */
+std::vector<char> awhDimParamSerialized(
+ AwhCoordinateProviderType inputCoordinateProvider = AwhCoordinateProviderType::Pull,
+ int inputCoordIndex = 0,
+ double inputOrigin = 0.5,
+ double inputEnd = 1.5,
+ double inputPeriod = 0,
+ // Correction for removal of GaussianGeometryFactor/2 in histogram size
+ double inputDiffusion = 0.1 / (0.144129616073222 * 2));
+
+/*! \internal \brief
+ * Struct that gathers all input for setting up and using a Bias
+ */
+struct AwhTestParameters
+{
+ explicit AwhTestParameters(ISerializer* serializer);
+ //! Move constructor
+ AwhTestParameters(AwhTestParameters&& o) noexcept :
+ beta(o.beta),
+ awhParams(std::move(o.awhParams)),
+ dimParams(std::move(o.dimParams))
+ {
+ }
+ //! 1/(kB*T).
+ double beta;
+
+ //! AWH parameters, this is the struct to actually use.
+ AwhParams awhParams;
+ //! Dimension parameters for setting up Bias.
+ std::vector<DimParams> dimParams;
+};
+
+/*! \brief
+ * Helper function to set up the C-style AWH parameters for the test.
+ *
+ * Builds the test input data from serialized data.
+ */
+AwhTestParameters getAwhTestParameters(AwhHistogramGrowthType eawhgrowth,
+ AwhPotentialType eawhpotential,
+ ArrayRef<const std::vector<char>> dimensionParameterBuffers,
+ bool inputUserData,
+ double beta,
+ bool useAwhFep,
+ double inputErrorScaling,
+ int numFepLambdaStates);
+
+} // namespace test
+} // namespace gmx
+
+#endif
#include <vector>
#include <gmock/gmock.h>
+#include <gmock/gmock-matchers.h>
#include <gtest/gtest.h>
#include "gromacs/applied_forces/awh/correlationgrid.h"
#include "gromacs/applied_forces/awh/pointstate.h"
#include "gromacs/mdtypes/awh_params.h"
+#include "gromacs/utility/inmemoryserializer.h"
#include "gromacs/utility/stringutil.h"
+#include "gromacs/applied_forces/awh/tests/awh_setup.h"
#include "testutils/refdata.h"
#include "testutils/testasserts.h"
namespace test
{
-/*! \internal \brief
- * Struct that gathers all input for setting up and using a Bias
- */
-struct AwhTestParameters
-{
- AwhTestParameters() = default;
- //! Move constructor
- AwhTestParameters(AwhTestParameters&& o) noexcept :
- beta(o.beta),
- awhDimParams(o.awhDimParams),
- awhBiasParams(o.awhBiasParams),
- awhParams(o.awhParams),
- dimParams(std::move(o.dimParams))
- {
- awhBiasParams.dimParams = &awhDimParams;
- awhParams.awhBiasParams = &awhBiasParams;
- }
- double beta; //!< 1/(kB*T)
-
- AwhDimParams awhDimParams; //!< Dimension parameters pointed to by \p awhBiasParams
- AwhBiasParams awhBiasParams; //!< Bias parameters pointed to by \[ awhParams
- AwhParams awhParams; //!< AWH parameters, this is the struct to actually use
-
- std::vector<DimParams> dimParams; //!< Dimension parameters for setting up Bias
-};
-
-//! Helper function to set up the C-style AWH parameters for the test
-static AwhTestParameters getAwhTestParameters(AwhHistogramGrowthType eawhgrowth, AwhPotentialType eawhpotential)
-{
- AwhTestParameters params;
-
- params.beta = 0.4;
-
- AwhDimParams& awhDimParams = params.awhDimParams;
-
- awhDimParams.period = 0;
- // Correction for removal of GaussianGeometryFactor/2 in histogram size
- awhDimParams.diffusion = 0.1 / (0.144129616073222 * 2);
- awhDimParams.origin = 0.5;
- awhDimParams.end = 1.5;
- awhDimParams.coordValueInit = awhDimParams.origin;
- awhDimParams.coverDiameter = 0;
- awhDimParams.eCoordProvider = AwhCoordinateProviderType::Pull;
-
- AwhBiasParams& awhBiasParams = params.awhBiasParams;
-
- awhBiasParams.ndim = 1;
- awhBiasParams.dimParams = &awhDimParams;
- awhBiasParams.eTarget = AwhTargetType::Constant;
- awhBiasParams.targetBetaScaling = 0;
- awhBiasParams.targetCutoff = 0;
- awhBiasParams.eGrowth = eawhgrowth;
- awhBiasParams.bUserData = FALSE;
- awhBiasParams.errorInitial = 0.5 / params.beta;
- awhBiasParams.shareGroup = 0;
- awhBiasParams.equilibrateHistogram = FALSE;
-
- double convFactor = 1;
- double k = 1000;
- int64_t seed = 93471803;
-
- params.dimParams.push_back(DimParams::pullDimParams(convFactor, k, params.beta));
-
- AwhParams& awhParams = params.awhParams;
-
- awhParams.numBias = 1;
- awhParams.awhBiasParams = &awhBiasParams;
- awhParams.seed = seed;
- awhParams.nstOut = 0;
- awhParams.nstSampleCoord = 1;
- awhParams.numSamplesUpdateFreeEnergy = 10;
- awhParams.ePotential = eawhpotential;
- awhParams.shareBiasMultisim = FALSE;
-
- return params;
-}
-
//! Database of 21 test coordinates that represent a trajectory */
-const double g_coords[] = { 0.62, 0.70, 0.68, 0.80, 0.93, 0.87, 1.16, 1.14, 0.95, 0.89, 0.91,
- 0.86, 0.88, 0.79, 0.75, 0.82, 0.74, 0.70, 0.68, 0.71, 0.73 };
+constexpr double g_coords[] = { 0.62, 0.70, 0.68, 0.80, 0.93, 0.87, 1.16, 1.14, 0.95, 0.89, 0.91,
+ 0.86, 0.88, 0.79, 0.75, 0.82, 0.74, 0.70, 0.68, 0.71, 0.73 };
//! Convenience typedef: growth type enum, potential type enum, disable update skips
typedef std::tuple<AwhHistogramGrowthType, AwhPotentialType, BiasParams::DisableUpdateSkips> BiasTestParameters;
*/
class BiasTest : public ::testing::TestWithParam<BiasTestParameters>
{
+private:
+ //! Storage for test parameters.
+ std::unique_ptr<AwhTestParameters> params_;
+
public:
//! Random seed for AWH MC sampling
int64_t seed_;
-
//! Coordinates representing a trajectory in time
std::vector<double> coordinates_;
//! The awh Bias
* The idea is to, among other things, have part of the interval
* not covered by samples.
*/
- const AwhTestParameters params = getAwhTestParameters(eawhgrowth, eawhpotential);
+ auto awhDimBuffer = awhDimParamSerialized();
+ auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
+ params_ = std::make_unique<AwhTestParameters>(getAwhTestParameters(
+ eawhgrowth, eawhpotential, awhDimArrayRef, 0, 0.4, false, 0.5, 0));
- seed_ = params.awhParams.seed;
+ seed_ = params_->awhParams.seed();
- double mdTimeStep = 0.1;
+ constexpr double mdTimeStep = 0.1;
int numSamples = coordinates_.size() - 1; // No sample taken at step 0
- GMX_RELEASE_ASSERT(numSamples % params.awhParams.numSamplesUpdateFreeEnergy == 0,
+ GMX_RELEASE_ASSERT(numSamples % params_->awhParams.numSamplesUpdateFreeEnergy() == 0,
"This test is intended to reproduce the situation when the might need "
"to write output during a normal AWH run, therefore the number of "
"samples should be a multiple of the free-energy update interval (but "
"the test should also runs fine without this condition).");
bias_ = std::make_unique<Bias>(-1,
- params.awhParams,
- params.awhBiasParams,
- params.dimParams,
- params.beta,
+ params_->awhParams,
+ params_->awhParams.awhBiasParams()[0],
+ params_->dimParams,
+ params_->beta,
mdTimeStep,
1,
"",
// Test that we detect coverings and exit the initial stage at the correct step
TEST(BiasTest, DetectsCovering)
{
+ const std::vector<char> serializedAwhParametersPerDim = awhDimParamSerialized();
+ auto awhDimArrayRef = gmx::arrayRefFromArray(&serializedAwhParametersPerDim, 1);
const AwhTestParameters params = getAwhTestParameters(AwhHistogramGrowthType::ExponentialLinear,
- AwhPotentialType::Convolved);
- const AwhDimParams& awhDimParams = params.awhParams.awhBiasParams[0].dimParams[0];
+ AwhPotentialType::Convolved,
+ awhDimArrayRef,
+ 0,
+ 0.4,
+ false,
+ 0.5,
+ 0);
+ const AwhDimParams& awhDimParams = params.awhParams.awhBiasParams()[0].dimParams()[0];
- const double mdTimeStep = 0.1;
+ constexpr double mdTimeStep = 0.1;
Bias bias(-1,
params.awhParams,
- params.awhBiasParams,
+ params.awhParams.awhBiasParams()[0],
params.dimParams,
params.beta,
mdTimeStep,
* coordinate range in a semi-realistic way. The period is 4*pi=12.57.
* We get out of the initial stage after 4 coverings at step 300.
*/
- const int64_t exitStepRef = 300;
- const double midPoint = 0.5 * (awhDimParams.end + awhDimParams.origin);
- const double halfWidth = 0.5 * (awhDimParams.end - awhDimParams.origin);
+ constexpr int64_t exitStepRef = 300;
+ const double midPoint = 0.5 * (awhDimParams.end() + awhDimParams.origin());
+ const double halfWidth = 0.5 * (awhDimParams.end() - awhDimParams.origin());
bool inInitialStage = bias.state().inInitialStage();
/* Normally this loop exits at exitStepRef, but we extend with failure */
nullptr,
step,
step,
- params.awhParams.seed,
+ params.awhParams.seed(),
nullptr);
inInitialStage = bias.state().inInitialStage();
#include "gromacs/mdtypes/awh_params.h"
#include "gromacs/utility/stringutil.h"
+#include "gromacs/applied_forces/awh/tests/awh_setup.h"
#include "testutils/refdata.h"
#include "testutils/testasserts.h"
{
//! The number of lambda states to use in the tests.
-const int numLambdaStates = 16;
+static constexpr int c_numLambdaStates = 16;
-/*! \internal \brief
- * Struct that gathers all input for setting up and using a Bias
- */
-struct AwhFepLambdaStateTestParameters
-{
- AwhFepLambdaStateTestParameters() = default;
- //! Move constructor
- AwhFepLambdaStateTestParameters(AwhFepLambdaStateTestParameters&& o) noexcept :
- beta(o.beta),
- awhDimParams(o.awhDimParams),
- awhBiasParams(o.awhBiasParams),
- awhParams(o.awhParams),
- dimParams(std::move(o.dimParams))
- {
- awhBiasParams.dimParams = &awhDimParams;
- awhParams.awhBiasParams = &awhBiasParams;
- }
- double beta; //!< 1/(kB*T)
-
- AwhDimParams awhDimParams; //!< Dimension parameters pointed to by \p awhBiasParams
- AwhBiasParams awhBiasParams; //!< Bias parameters pointed to by \[ awhParams
- AwhParams awhParams; //!< AWH parameters, this is the struct to actually use
-
- std::vector<DimParams> dimParams; //!< Dimension parameters for setting up Bias
-};
-
-//! Helper function to set up the C-style AWH parameters for the test
-static AwhFepLambdaStateTestParameters getAwhFepLambdaTestParameters(AwhHistogramGrowthType eawhgrowth,
- AwhPotentialType eawhpotential)
-{
- AwhFepLambdaStateTestParameters params;
-
- params.beta = 0.4;
-
- AwhDimParams& awhDimParams = params.awhDimParams;
-
- awhDimParams.period = 0;
- // Correction for removal of GaussianGeometryFactor/2 in histogram size
- awhDimParams.diffusion = 1e-4 / (0.12927243028700 * 2);
- awhDimParams.origin = 0;
- awhDimParams.end = numLambdaStates - 1;
- awhDimParams.coordValueInit = awhDimParams.origin;
- awhDimParams.coverDiameter = 0;
- awhDimParams.eCoordProvider = AwhCoordinateProviderType::FreeEnergyLambda;
-
- AwhBiasParams& awhBiasParams = params.awhBiasParams;
-
- awhBiasParams.ndim = 1;
- awhBiasParams.dimParams = &awhDimParams;
- awhBiasParams.eTarget = AwhTargetType::Constant;
- awhBiasParams.targetBetaScaling = 0;
- awhBiasParams.targetCutoff = 0;
- awhBiasParams.eGrowth = eawhgrowth;
- awhBiasParams.bUserData = FALSE;
- awhBiasParams.errorInitial = 1.0 / params.beta;
- awhBiasParams.shareGroup = 0;
- awhBiasParams.equilibrateHistogram = FALSE;
-
- int64_t seed = 93471803;
-
- params.dimParams.push_back(DimParams::fepLambdaDimParams(numLambdaStates, params.beta));
-
- AwhParams& awhParams = params.awhParams;
-
- awhParams.numBias = 1;
- awhParams.awhBiasParams = &awhBiasParams;
- awhParams.seed = seed;
- awhParams.nstOut = 0;
- awhParams.nstSampleCoord = 1;
- awhParams.numSamplesUpdateFreeEnergy = 10;
- awhParams.ePotential = eawhpotential;
- awhParams.shareBiasMultisim = FALSE;
-
- return params;
-}
//! Convenience typedef: growth type enum, potential type enum, disable update skips
typedef std::tuple<AwhHistogramGrowthType, AwhPotentialType, BiasParams::DisableUpdateSkips> BiasTestParameters;
*/
class BiasFepLambdaStateTest : public ::testing::TestWithParam<BiasTestParameters>
{
+private:
+ //! Storage for test parameters.
+ std::unique_ptr<AwhTestParameters> params_;
+
public:
//! Random seed for AWH MC sampling
int64_t seed_;
/* Set up a basic AWH setup with a single, 1D bias with parameters
* such that we can measure the effects of different parameters.
*/
- const AwhFepLambdaStateTestParameters params =
- getAwhFepLambdaTestParameters(eawhgrowth, eawhpotential);
-
- seed_ = params.awhParams.seed;
+ constexpr AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::FreeEnergyLambda;
+ constexpr int coordIndex = 0;
+ constexpr double origin = 0;
+ constexpr double end = c_numLambdaStates - 1;
+ constexpr double period = 0;
+ // Correction for removal of GaussianGeometryFactor/2 in histogram size
+ constexpr double diffusion = 1e-4 / (0.12927243028700 * 2);
+ const auto awhDimBuffer =
+ awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
+ const auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
+ params_ = std::make_unique<AwhTestParameters>(getAwhTestParameters(
+ eawhgrowth, eawhpotential, awhDimArrayRef, false, 0.4, true, 1.0, c_numLambdaStates));
+
+ seed_ = params_->awhParams.seed();
double mdTimeStep = 0.1;
bias_ = std::make_unique<Bias>(-1,
- params.awhParams,
- params.awhBiasParams,
- params.dimParams,
- params.beta,
+ params_->awhParams,
+ params_->awhParams.awhBiasParams()[0],
+ params_->dimParams,
+ params_->beta,
mdTimeStep,
1,
"",
int nSteps = 501;
/* Some energies to use as base values (to which some noise is added later on). */
- std::vector<double> neighborLambdaEnergies(numLambdaStates);
- std::vector<double> neighborLambdaDhdl(numLambdaStates);
+ std::vector<double> neighborLambdaEnergies(c_numLambdaStates);
+ std::vector<double> neighborLambdaDhdl(c_numLambdaStates);
const double magnitude = 12.0;
- for (int i = 0; i < numLambdaStates; i++)
+ for (int i = 0; i < c_numLambdaStates; i++)
{
neighborLambdaEnergies[i] = magnitude * std::sin(i * 0.1);
neighborLambdaDhdl[i] = magnitude * std::cos(i * 0.1);
// Test that we detect coverings and exit the initial stage at the correct step
TEST(BiasFepLambdaStateTest, DetectsCovering)
{
- const AwhFepLambdaStateTestParameters params = getAwhFepLambdaTestParameters(
- AwhHistogramGrowthType::ExponentialLinear, AwhPotentialType::Convolved);
+ constexpr AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::FreeEnergyLambda;
+ constexpr int coordIndex = 0;
+ constexpr double origin = 0;
+ constexpr double end = c_numLambdaStates - 1;
+ constexpr double period = 0;
+ constexpr double diffusion = 1e-4 / (0.12927243028700 * 2);
+ auto awhDimBuffer =
+ awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
+ auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
+ const AwhTestParameters params(getAwhTestParameters(AwhHistogramGrowthType::ExponentialLinear,
+ AwhPotentialType::Convolved,
+ awhDimArrayRef,
+ false,
+ 0.4,
+ true,
+ 1.0,
+ c_numLambdaStates));
const double mdTimeStep = 0.1;
Bias bias(-1,
params.awhParams,
- params.awhBiasParams,
+ params.awhParams.awhBiasParams()[0],
params.dimParams,
params.beta,
mdTimeStep,
bool inInitialStage = bias.state().inInitialStage();
/* Some energies to use as base values (to which some noise is added later on). */
- std::vector<double> neighborLambdaEnergies(numLambdaStates);
- std::vector<double> neighborLambdaDhdl(numLambdaStates);
+ std::vector<double> neighborLambdaEnergies(c_numLambdaStates);
+ std::vector<double> neighborLambdaDhdl(c_numLambdaStates);
const double magnitude = 12.0;
- for (int i = 0; i < numLambdaStates; i++)
+ for (int i = 0; i < c_numLambdaStates; i++)
{
neighborLambdaEnergies[i] = magnitude * std::sin(i * 0.1);
neighborLambdaDhdl[i] = magnitude * std::cos(i * 0.1);
nullptr,
step,
step,
- params.awhParams.seed,
+ params.awhParams.seed(),
nullptr);
inInitialStage = bias.state().inInitialStage();
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2017,2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/mdtypes/awh_params.h"
#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/inmemoryserializer.h"
+#include "gromacs/applied_forces/awh/tests/awh_setup.h"
#include "testutils/testasserts.h"
namespace gmx
*/
const int numDim = 2;
- std::vector<AwhDimParams> awhDimParams(numDim);
-
- awhDimParams[0].origin = -5;
- awhDimParams[0].end = 5;
- awhDimParams[0].period = 10;
-
- awhDimParams[1].origin = 0.5;
- awhDimParams[1].end = 2.0;
- awhDimParams[1].period = 0;
+ std::vector<AwhDimParams> awhDimParams;
+ AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::Pull;
+ double diffusion = 0.1;
+ {
+ int coordIndex = 0;
+ double origin = -5;
+ double end = 5;
+ double period = 10;
+ auto awhDimBuffer =
+ awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
+ gmx::InMemoryDeserializer serializer(awhDimBuffer, false);
+ awhDimParams.emplace_back(AwhDimParams(&serializer));
+ }
+ {
+ int coordIndex = 1;
+ double origin = 0.5;
+ double end = 2;
+ double period = 0;
+ auto awhDimBuffer =
+ awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
+ gmx::InMemoryDeserializer serializer(awhDimBuffer, false);
+ awhDimParams.emplace_back(AwhDimParams(&serializer));
+ }
const real conversionFactor = 1;
const real beta = 3.0;
dimParams.push_back(DimParams::pullDimParams(conversionFactor, 1 / (beta * 0.7 * 0.7), beta));
dimParams.push_back(DimParams::pullDimParams(conversionFactor, 1 / (beta * 0.1 * 0.1), beta));
- BiasGrid grid(dimParams, awhDimParams.data());
+ BiasGrid grid(dimParams, awhDimParams);
const int numPoints = grid.numPoints();
#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/smalloc.h"
+#include "gromacs/applied_forces/awh/tests/awh_setup.h"
#include "testutils/testasserts.h"
#include "testutils/testfilemanager.h"
namespace test
{
-/*! \internal \brief
- * Struct that gathers all input for setting up and using a Bias
- */
-struct AwhTestParameters
-{
- double beta; //!< 1/(kB*T)
-
- AwhDimParams awhDimParams[2]; //!< Dimension parameters pointed to by \p awhBiasParams
- AwhBiasParams awhBiasParams; //!< Bias parameters pointed to by \[ awhParams
- AwhParams awhParams; //!< AWH parameters, this is the struct to actually use
-};
-
-//! Helper function to set up the C-style AWH parameters for the test
-static AwhTestParameters getAwhTestParameters()
-{
- AwhTestParameters params;
-
- params.beta = 1.0;
-
- AwhParams& awhParams = params.awhParams;
- snew(params.awhParams.awhBiasParams, 1);
- AwhBiasParams& awhBiasParams = params.awhParams.awhBiasParams[0];
- snew(awhBiasParams.dimParams, 2);
-
- AwhDimParams& awhDimParams0 = awhBiasParams.dimParams[0];
-
- awhDimParams0.period = 0;
- awhDimParams0.diffusion = 0.1;
- awhDimParams0.origin = 0.5;
- awhDimParams0.end = 1.5;
- awhDimParams0.coordValueInit = awhDimParams0.origin;
- awhDimParams0.coverDiameter = 0;
- awhDimParams0.eCoordProvider = AwhCoordinateProviderType::Pull;
-
- AwhDimParams& awhDimParams1 = awhBiasParams.dimParams[1];
-
- awhDimParams1.period = 0;
- awhDimParams1.diffusion = 0.1;
- awhDimParams1.origin = 0.8;
- awhDimParams1.end = 1.3;
- awhDimParams1.coordValueInit = awhDimParams1.origin;
- awhDimParams1.coverDiameter = 0;
- awhDimParams1.eCoordProvider = AwhCoordinateProviderType::Pull;
-
- awhBiasParams.ndim = 2;
- awhBiasParams.eTarget = AwhTargetType::Constant;
- awhBiasParams.targetBetaScaling = 0;
- awhBiasParams.targetCutoff = 0;
- awhBiasParams.eGrowth = AwhHistogramGrowthType::Linear;
- awhBiasParams.bUserData = TRUE;
- awhBiasParams.errorInitial = 0.5;
- awhBiasParams.shareGroup = 0;
- awhBiasParams.equilibrateHistogram = FALSE;
-
- awhParams.numBias = 1;
- awhParams.seed = 93471803;
- awhParams.nstOut = 0;
- awhParams.nstSampleCoord = 1;
- awhParams.numSamplesUpdateFreeEnergy = 10;
- awhParams.ePotential = AwhPotentialType::Convolved;
- awhParams.shareBiasMultisim = FALSE;
-
- return params;
-}
-
/*! \brief Test fixture for testing Bias updates
*/
class BiasStateTest : public ::testing::TestWithParam<const char*>
{
+private:
+ std::unique_ptr<AwhTestParameters> params_;
+
public:
std::unique_ptr<BiasState> biasState_; //!< The bias state
BiasStateTest()
{
- AwhTestParameters params = getAwhTestParameters();
- const AwhParams& awhParams = params.awhParams;
- const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[0];
+ std::vector<std::vector<char>> awhDimParameters;
+ AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::Pull;
+ double diffusion = 0.1;
+ {
+ int coordIndex = 0;
+ double origin = 0.5;
+ double end = 1.5;
+ double period = 0;
+ awhDimParameters.emplace_back(awhDimParamSerialized(
+ coordinateProvider, coordIndex, origin, end, period, diffusion));
+ }
+ {
+ int coordIndex = 1;
+ double origin = 0.8;
+ double end = 1.3;
+ double period = 0;
+ awhDimParameters.emplace_back(awhDimParamSerialized(
+ coordinateProvider, coordIndex, origin, end, period, diffusion));
+ }
+ params_ = std::make_unique<AwhTestParameters>(getAwhTestParameters(
+ AwhHistogramGrowthType::Linear, AwhPotentialType::Convolved, awhDimParameters, 1, 1.0, false, 0.5, 0));
+ const AwhParams& awhParams = params_->awhParams;
+ const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams()[0];
std::vector<DimParams> dimParams;
- dimParams.push_back(DimParams::pullDimParams(1.0, 15.0, params.beta));
- dimParams.push_back(DimParams::pullDimParams(1.0, 15.0, params.beta));
- BiasGrid grid(dimParams, awhBiasParams.dimParams);
+ dimParams.push_back(DimParams::pullDimParams(1.0, 15.0, params_->beta));
+ dimParams.push_back(DimParams::pullDimParams(1.0, 15.0, params_->beta));
+ BiasGrid grid(dimParams, awhBiasParams.dimParams());
BiasParams biasParams(
awhParams, awhBiasParams, dimParams, 1.0, 1.0, BiasParams::DisableUpdateSkips::no, 1, grid.axis(), 0);
biasState_ = std::make_unique<BiasState>(awhBiasParams, 1.0, dimParams, grid);
// Here we initialize the grid point state using the input file
std::string filename = gmx::test::TestFileManager::getInputFilePath(GetParam());
biasState_->initGridPointState(
- awhBiasParams, dimParams, grid, biasParams, filename, params.awhParams.numBias);
-
- sfree(params.awhParams.awhBiasParams[0].dimParams);
- sfree(params.awhParams.awhBiasParams);
+ awhBiasParams, dimParams, grid, biasParams, filename, params_->awhParams.numBias());
}
};
#include <memory>
#include <vector>
+#include "gromacs/applied_forces/awh/read_params.h"
#include "gromacs/fileio/filetypes.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/gmxfio_xdr.h"
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/inmemoryserializer.h"
+#include "gromacs/utility/iserializer.h"
#include "gromacs/utility/keyvaluetreebuilder.h"
#include "gromacs/utility/keyvaluetreeserializer.h"
#include "gromacs/utility/smalloc.h"
}
}
-static void do_awhBias(gmx::ISerializer* serializer, gmx::AwhBiasParams* awhBiasParams, int gmx_unused file_version)
-{
- serializer->doEnumAsInt(&awhBiasParams->eTarget);
- serializer->doDouble(&awhBiasParams->targetBetaScaling);
- serializer->doDouble(&awhBiasParams->targetCutoff);
- serializer->doEnumAsInt(&awhBiasParams->eGrowth);
- if (serializer->reading())
- {
- int temp = 0;
- serializer->doInt(&temp);
- awhBiasParams->bUserData = static_cast<bool>(temp);
- }
- else
- {
- int temp = static_cast<int>(awhBiasParams->bUserData);
- serializer->doInt(&temp);
- }
- serializer->doDouble(&awhBiasParams->errorInitial);
- serializer->doInt(&awhBiasParams->ndim);
- serializer->doInt(&awhBiasParams->shareGroup);
- serializer->doBool(&awhBiasParams->equilibrateHistogram);
-
- if (serializer->reading())
- {
- snew(awhBiasParams->dimParams, awhBiasParams->ndim);
- }
-
- for (int d = 0; d < awhBiasParams->ndim; d++)
- {
- gmx::AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
-
- serializer->doEnumAsInt(&dimParams->eCoordProvider);
- serializer->doInt(&dimParams->coordIndex);
- serializer->doDouble(&dimParams->origin);
- serializer->doDouble(&dimParams->end);
- serializer->doDouble(&dimParams->period);
- serializer->doDouble(&dimParams->forceConstant);
- serializer->doDouble(&dimParams->diffusion);
- serializer->doDouble(&dimParams->coordValueInit);
- serializer->doDouble(&dimParams->coverDiameter);
- }
-}
-
-static void do_awh(gmx::ISerializer* serializer, gmx::AwhParams* awhParams, int gmx_unused file_version)
-{
- serializer->doInt(&awhParams->numBias);
- serializer->doInt(&awhParams->nstOut);
- serializer->doInt64(&awhParams->seed);
- serializer->doInt(&awhParams->nstSampleCoord);
- serializer->doInt(&awhParams->numSamplesUpdateFreeEnergy);
- serializer->doEnumAsInt(&awhParams->ePotential);
- serializer->doBool(&awhParams->shareBiasMultisim);
-
- if (awhParams->numBias > 0)
- {
- if (serializer->reading())
- {
- snew(awhParams->awhBiasParams, awhParams->numBias);
- }
-
- for (int k = 0; k < awhParams->numBias; k++)
- {
- do_awhBias(serializer, &awhParams->awhBiasParams[k], file_version);
- }
- }
-}
-
static void do_pull(gmx::ISerializer* serializer, pull_params_t* pull, int file_version, PullingAlgorithm ePullOld)
{
PullGroupGeometry eGeomOld = PullGroupGeometry::Count;
{
if (serializer->reading())
{
- snew(ir->awhParams, 1);
+ ir->awhParams = std::make_unique<gmx::AwhParams>(serializer);
+ }
+ else
+ {
+ ir->awhParams->serialize(serializer);
}
- do_awh(serializer, ir->awhParams, file_version);
}
}
else
#include <memory>
#include <string>
+#include "gromacs/applied_forces/awh/read_params.h"
#include "gromacs/commandline/pargs.h"
#include "gromacs/fileio/enxio.h"
#include "gromacs/fileio/gmxfio.h"
*/
void initializeAwhOutputFile(int subBlockStart,
int numSubBlocks,
- const AwhBiasParams* awhBiasParams,
+ const AwhBiasParams& awhBiasParams,
AwhGraphSelection graphSelection,
EnergyUnit energyUnit,
real kTValue);
*/
void initializeFrictionOutputFile(int subBlockStart,
int numSubBlocks,
- const AwhBiasParams* awhBiasParams,
+ const AwhBiasParams& awhBiasParams,
EnergyUnit energyUnit,
real kTValue);
{
public:
//! Constructor
- AwhReader(const AwhParams* awhParams,
+ AwhReader(const AwhParams& awhParams,
int numFileOptions,
const t_filenm* filenames,
AwhGraphSelection awhGraphSelection,
constexpr int maxAwhGraphs = 6;
/*! \brief Constructs a legend for a standard awh output file */
-std::vector<std::string> makeLegend(const AwhBiasParams* awhBiasParams,
+std::vector<std::string> makeLegend(const AwhBiasParams& awhBiasParams,
OutputFileType outputFileType,
size_t numLegend)
{
std::vector<std::string> legend;
/* Give legends to dimensions higher than the first */
- for (int d = 1; d < awhBiasParams->ndim; d++)
+ for (int d = 1; d < awhBiasParams.ndim(); d++)
{
legend.push_back(gmx::formatString("dim%d", d));
}
}
break;
case OutputFileType::Friction:
- for (int i0 = 0; i0 < awhBiasParams->ndim; i0++)
+ for (int i0 = 0; i0 < awhBiasParams.ndim(); i0++)
{
for (int i1 = 0; i1 <= i0; i1++)
{
void OutputFile::initializeAwhOutputFile(int subblockStart,
int numSubBlocks,
- const AwhBiasParams* awhBiasParams,
+ const AwhBiasParams& awhBiasParams,
AwhGraphSelection graphSelection,
EnergyUnit energyUnit,
real kTValue)
{
/* The first subblock with actual graph y-values is index 1 + ndim */
- numDim_ = awhBiasParams->ndim;
+ numDim_ = awhBiasParams.ndim();
firstGraphSubBlock_ = subblockStart + 1 + numDim_;
if (graphSelection == AwhGraphSelection::All)
{
if (graphSelection == AwhGraphSelection::All)
{
yLabel_ += gmx::formatString(
- ", (nm\\S-%d\\N or rad\\S-%d\\N), (-)", awhBiasParams->ndim, awhBiasParams->ndim);
+ ", (nm\\S-%d\\N or rad\\S-%d\\N), (-)", awhBiasParams.ndim(), awhBiasParams.ndim());
}
}
/*! \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,
+ const AwhBiasParams& awhBiasParams,
EnergyUnit energyUnit,
real kTValue)
{
/* The first subblock with actual graph y-values is index 1 + ndim */
- numDim_ = awhBiasParams->ndim;
+ numDim_ = awhBiasParams.ndim();
int numTensorElements = (numDim_ * (numDim_ + 1)) / 2;
/* The friction tensor elements are always the last subblocks */
}
}
-AwhReader::AwhReader(const AwhParams* awhParams,
+AwhReader::AwhReader(const AwhParams& awhParams,
int numFileOptions,
const t_filenm* filenames,
AwhGraphSelection awhGraphSelection,
/* Keep track of the first subblock of this AWH */
int subblockStart = 0;
- for (int k = 0; k < awhParams->numBias; k++)
+ for (int k = 0; k < awhParams.numBias(); k++)
{
- AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
+ const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams()[k];
int numSubBlocks = static_cast<int>(block->sub[subblockStart].fval[0]);
std::unique_ptr<OutputFile> awhOutputFile(new OutputFile(
- opt2fn("-o", numFileOptions, filenames), "AWH", awhParams->numBias, k));
+ opt2fn("-o", numFileOptions, filenames), "AWH", awhParams.numBias(), k));
awhOutputFile->initializeAwhOutputFile(
subblockStart, numSubBlocks, awhBiasParams, awhGraphSelection, energyUnit, kT);
if (outputFriction)
{
frictionOutputFile = std::make_unique<OutputFile>(
- opt2fn("-fric", numFileOptions, filenames), "Friction tensor", awhParams->numBias, k);
+ opt2fn("-fric", numFileOptions, filenames), "Friction tensor", awhParams.numBias(), k);
frictionOutputFile->initializeFrictionOutputFile(
subblockStart, numSubBlocks, awhBiasParams, energyUnit, kT);
AwhGraphSelection awhGraphSelection =
(moreGraphs ? AwhGraphSelection::All : AwhGraphSelection::Pmf);
EnergyUnit energyUnit = (kTUnit ? EnergyUnit::KT : EnergyUnit::KJPerMol);
- awhReader = std::make_unique<AwhReader>(ir.awhParams,
+ awhReader = std::make_unique<AwhReader>(*ir.awhParams,
nfile,
fnm,
awhGraphSelection,
copy_mat(ir->compress, compressibility);
}
setStateDependentAwhParams(
- ir->awhParams,
+ ir->awhParams.get(),
*ir->pull,
pull,
state.box,
*/
#include "gmxpre.h"
-#include "gromacs/utility/enumerationhelpers.h"
#include "readir.h"
#include <cctype>
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/calc_verletbuf.h"
#include "gromacs/mdrun/mdmodules.h"
+#include "gromacs/mdtypes/awh_params.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/multipletimestepping.h"
ir->bDoAwh = (getEnum<Boolean>(&inp, "awh", wi) != Boolean::No);
if (ir->bDoAwh)
{
- ir->awhParams = gmx::readAwhParams(&inp, wi);
+ ir->awhParams = std::make_unique<gmx::AwhParams>(&inp, *ir, wi);
}
/* Enforced rotation */
if (ir->bDoAwh)
{
- gmx::checkAwhParams(ir->awhParams, ir, wi);
+ gmx::checkAwhParams(*ir->awhParams, *ir, wi);
}
sfree(dumstr[0]);
#include <numeric>
#include "gromacs/applied_forces/awh/awh.h"
+#include "gromacs/applied_forces/awh/read_params.h"
#include "gromacs/commandline/filenm.h"
#include "gromacs/domdec/collect.h"
#include "gromacs/domdec/dlbtiming.h"
// a user can't just do multi-sim with single-sim orientation restraints.
bool usingEnsembleRestraints =
(fcdata.disres->nsystems > 1) || ((ms != nullptr) && (fcdata.orires->nr != 0));
- bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && (ms != nullptr));
+ bool awhUsesMultiSim = (ir->bDoAwh && ir->awhParams->shareBiasMultisim() && (ms != nullptr));
// Replica exchange, ensemble restraints and AWH need all
// simulations to remain synchronized, so they need
}
if (ir->bDoAwh)
{
- nstfep = std::gcd(ir->awhParams->nstSampleCoord, nstfep);
+ nstfep = std::gcd(ir->awhParams->nstSampleCoord(), nstfep);
}
}
#ifndef GMX_MDTYPES_AWH_PARAMS_H
#define GMX_MDTYPES_AWH_PARAMS_H
+#include <vector>
+
#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/classhelpers.h"
+
+struct t_inpfile;
+struct t_inputrec;
+struct pull_params_t;
+using warninp_t = struct warninp*;
namespace gmx
{
+class ISerializer;
//! Target distribution enum.
enum class AwhTargetType : int
{
//! String for AWH bias reaction coordinate provider.
const char* enumValueToString(AwhCoordinateProviderType enumValue);
-/*! \cond INTERNAL */
-
-//! Parameters for an AWH coordinate dimension.
-struct AwhDimParams
+class AwhDimParams
{
- AwhCoordinateProviderType eCoordProvider; /**< The module providing the reaction coordinate. */
- int coordIndex; /**< Index of reaction coordinate in the provider. */
- double origin; /**< Start value of the interval. */
- double end; /**< End value of the interval. */
- double period; /**< The period of this dimension (= 0 if not periodic). */
- double forceConstant; /**< The force constant in kJ/mol/nm^2, kJ/mol/rad^2 */
- double diffusion; /**< Estimated diffusion constant in units of nm^2/ps, rad^2/ps or ps^-1. */
- double coordValueInit; /**< The initial coordinate value. */
- double coverDiameter; /**< The diameter that needs to be sampled around a point before it is considered covered. */
-};
+public:
+ //! Constructor from input file.
+ AwhDimParams(std::vector<t_inpfile>* inp,
+ const std::string& prefix,
+ const t_inputrec& ir,
+ warninp_t wi,
+ bool bComment);
+ //! Constructor to generate from file reading.
+ explicit AwhDimParams(ISerializer* serializer);
-//! Parameters for an AWH bias.
-struct AwhBiasParams
-{
- // TODO: Turn dimParams into a std::vector when moved into AWH module
- int ndim; /**< Dimension of the coordinate space. */
- AwhDimParams* dimParams; /**< AWH parameters per dimension. */
- AwhTargetType eTarget; /**< Type of target distribution. */
- double targetBetaScaling; /**< Beta scaling value for Boltzmann type target distributions. */
- double targetCutoff; /**< Free energy cutoff value for cutoff type target distribution in kJ/mol.*/
- AwhHistogramGrowthType eGrowth; /**< How the biasing histogram grows. */
- bool bUserData; /**< Is there a user-defined initial PMF estimate and target estimate? */
- double errorInitial; /**< Estimated initial free energy error in kJ/mol. */
- int shareGroup; /**< When >0, the bias is shared with biases of the same group and across multiple simulations when shareBiasMultisim=true */
- gmx_bool equilibrateHistogram; /**< True if the simulation starts out by equilibrating the histogram. */
+ //! Move constructor.
+ AwhDimParams(AwhDimParams&&) = default;
+ //! Move assignment operator.
+ AwhDimParams& operator=(AwhDimParams&&) = default;
+ //! Delete copy constructor.
+ AwhDimParams(const AwhDimParams&) = delete;
+ //! Delete copy assignment.
+ AwhDimParams& operator=(const AwhDimParams&) = delete;
+
+ //! Which module is providing the reaction coordinate.
+ AwhCoordinateProviderType coordinateProvider() const { return eCoordProvider_; }
+ //! Index for reaction coordinate in provider.
+ int coordinateIndex() const { return coordIndex_; }
+ //! Start value for interval.
+ double origin() const { return origin_; }
+ //! End value for interval.
+ double end() const { return end_; }
+ //! Period for the dimension.
+ double period() const { return period_; }
+ //! Set period value dependent on state.
+ void setPeriod(double period) { period_ = period; }
+ //! Force constant for this dimension.
+ double forceConstant() const { return forceConstant_; }
+ //! Estimated diffusion constant.
+ double diffusion() const { return diffusion_; }
+ //! Initial value for coordinate.
+ double initialCoordinate() const { return coordValueInit_; }
+ //! Set initial coordinate value dependent on state.
+ void setInitialCoordinate(double initialCoordinate) { coordValueInit_ = initialCoordinate; }
+ //! Diameter needed to be sampled.
+ double coverDiameter() const { return coverDiameter_; }
+ //! Write datastructure.
+ void serialize(ISerializer* serializer);
+
+private:
+ //! The module providing the reaction coordinate.
+ AwhCoordinateProviderType eCoordProvider_;
+ //! Index of reaction coordinate in the provider.
+ int coordIndex_;
+ //! Start value of the interval.
+ double origin_;
+ //! End value of the interval.
+ double end_;
+ //! The period of this dimension (= 0 if not periodic).
+ double period_;
+ //! The force constant in kJ/mol/nm^2, kJ/mol/rad^2
+ double forceConstant_;
+ //! Estimated diffusion constant in units of nm^2/ps or rad^2/ps or ps^-1.
+ double diffusion_;
+ //! The initial coordinate value.
+ double coordValueInit_;
+ //! The diameter that needs to be sampled around a point before it is considered covered.
+ double coverDiameter_;
};
-//! Parameters for AWH.
-struct AwhParams
+class AwhBiasParams
{
- // TODO: Turn awhBiasParams into a std::vector when moved into AWH module
- int numBias; /**< The number of AWH biases.*/
- AwhBiasParams* awhBiasParams; /**< AWH bias parameters.*/
- int64_t seed; /**< Random seed.*/
- int nstOut; /**< Output step interval.*/
- int nstSampleCoord; /**< Number of samples per coordinate sample (also used for PMF) */
- int numSamplesUpdateFreeEnergy; /**< Number of samples per free energy update. */
- AwhPotentialType ePotential; /**< Type of potential. */
- gmx_bool shareBiasMultisim; /**< When true, share biases with shareGroup>0 between multi-simulations */
+public:
+ //! Constructor from input file.
+ AwhBiasParams(std::vector<t_inpfile>* inp,
+ const std::string& prefix,
+ const t_inputrec& ir,
+ warninp_t wi,
+ bool bComment);
+ //! Constructor to generate from file reading.
+ explicit AwhBiasParams(ISerializer* serializer);
+
+ //! Move constructor.
+ AwhBiasParams(AwhBiasParams&&) = default;
+ //! Move assignment operator.
+ AwhBiasParams& operator=(AwhBiasParams&&) = default;
+ //! Delete copy constructor.
+ AwhBiasParams(const AwhBiasParams&) = delete;
+ //! Delete copy assignment.
+ AwhBiasParams& operator=(const AwhBiasParams&) = delete;
+
+ //! Which target distribution is searched.
+ AwhTargetType targetDistribution() const { return eTarget_; }
+ //! Beta scaling to reach target distribution.
+ double targetBetaScaling() const { return targetBetaScaling_; }
+ //! Cutoff for target.
+ double targetCutoff() const { return targetCutoff_; }
+ //! Which kind of growth to use.
+ AwhHistogramGrowthType growthType() const { return eGrowth_; }
+ //! User provided PMF estimate.
+ bool userPMFEstimate() const { return bUserData_; }
+ //! Estimated initial free energy error in kJ/mol.
+ double initialErrorEstimate() const { return errorInitial_; }
+ //! Dimensions of coordinate space.
+ int ndim() const { return dimParams_.size(); }
+ //! Number of groups to share this bias with.
+ int shareGroup() const { return shareGroup_; }
+ //! If the simulation starts with equilibrating histogram.
+ bool equilibrateHistogram() const { return equilibrateHistogram_; }
+ //! Access to dimension parameters.
+ ArrayRef<AwhDimParams> dimParams() { return dimParams_; }
+ //! Const access to dimension parameters.
+ ArrayRef<const AwhDimParams> dimParams() const { return dimParams_; }
+ //! Write datastructure.
+ void serialize(ISerializer* serializer);
+
+private:
+ //! AWH parameters per dimension.
+ std::vector<AwhDimParams> dimParams_;
+ //! Type of target distribution.
+ AwhTargetType eTarget_;
+ //! Beta scaling value for Boltzmann type target distributions.
+ double targetBetaScaling_;
+ //! Free energy cutoff value for cutoff type target distribution in kJ/mol.
+ double targetCutoff_;
+ //! How the biasing histogram grows.
+ AwhHistogramGrowthType eGrowth_;
+ //! Is there a user-defined initial PMF estimate and target estimate?
+ bool bUserData_;
+ //! Estimated initial free energy error in kJ/mol.
+ double errorInitial_;
+ //! When >0, the bias is shared with biases of the same group and across multiple simulations when shareBiasMultisim=true
+ int shareGroup_;
+ //! True if the simulation starts out by equilibrating the histogram.
+ bool equilibrateHistogram_;
};
+/*! \brief
+ * Structure holding parameter information for AWH.
+ */
+class AwhParams
+{
+public:
+ //! Constructor from input file.
+ AwhParams(std::vector<t_inpfile>* inp, const t_inputrec& ir, warninp_t wi);
+ //! Constructor used to generate awh parameter from file reading.
+ explicit AwhParams(ISerializer* serializer);
-/*! \endcond */
+ //! Move constructor.
+ AwhParams(AwhParams&&) = default;
+ //! Move assignment operator.
+ AwhParams& operator=(AwhParams&&) = default;
+ //! Delete copy constructor.
+ AwhParams(const AwhParams&) = delete;
+ //! Delete copy assignment.
+ AwhParams& operator=(const AwhParams&) = delete;
+
+ //! Get number of biases.
+ int numBias() const { return awhBiasParams_.size(); }
+ //! Get access to bias parameters.
+ ArrayRef<AwhBiasParams> awhBiasParams() { return awhBiasParams_; }
+ //! Const access to bias parameters.
+ ArrayRef<const AwhBiasParams> awhBiasParams() const { return awhBiasParams_; }
+ //! What king of potential is being used. \todo should use actual enum class.
+ AwhPotentialType potential() const { return potentialEnum_; }
+ //! Seed used for starting AWH.
+ int64_t seed() const { return seed_; }
+ //! Output step interval.
+ int nstout() const { return nstOut_; }
+ //! Number of samples per coordinate sample.
+ int nstSampleCoord() const { return nstSampleCoord_; }
+ //! Number of samples per free energy update.
+ int numSamplesUpdateFreeEnergy() const { return numSamplesUpdateFreeEnergy_; }
+ //! If biases are shared in multisim.
+ bool shareBiasMultisim() const { return shareBiasMultisim_; }
+ //! Serialize awh parameters.
+ void serialize(ISerializer* serializer);
+
+private:
+ //! AWH bias parameters.
+ std::vector<AwhBiasParams> awhBiasParams_;
+ //! Random seed.
+ int64_t seed_;
+ //! Output step interval.
+ int nstOut_;
+ //! Number of samples per coordinate sample (also used for PMF)
+ int nstSampleCoord_;
+ //! Number of samples per free energy update.
+ int numSamplesUpdateFreeEnergy_;
+ //! Type of potential.
+ AwhPotentialType potentialEnum_;
+ //! Whether to share biases with shareGroup>0 between multi-simulations.
+ bool shareBiasMultisim_;
+};
} // namespace gmx
#include <memory>
#include <numeric>
+#include "gromacs/applied_forces/awh/read_params.h"
#include "gromacs/math/veccompare.h"
#include "gromacs/math/vecdump.h"
#include "gromacs/mdtypes/awh_params.h"
}
}
-static void pr_awh_bias_dim(FILE* fp, int indent, gmx::AwhDimParams* awhDimParams, const char* prefix)
+static void pr_awh_bias_dim(FILE* fp, int indent, const gmx::AwhDimParams& awhDimParams, const char* prefix)
{
pr_indent(fp, indent);
indent++;
fprintf(fp, "%s:\n", prefix);
- PS("coord-provider", enumValueToString(awhDimParams->eCoordProvider));
- PI("coord-index", awhDimParams->coordIndex + 1);
- PR("start", awhDimParams->origin);
- PR("end", awhDimParams->end);
- PR("period", awhDimParams->period);
- PR("force-constant", awhDimParams->forceConstant);
- PR("diffusion", awhDimParams->diffusion);
- PR("cover-diameter", awhDimParams->coverDiameter);
+ PS("coord-provider", enumValueToString(awhDimParams.coordinateProvider()));
+ PI("coord-index", awhDimParams.coordinateIndex() + 1);
+ PR("start", awhDimParams.origin());
+ PR("end", awhDimParams.end());
+ PR("period", awhDimParams.period());
+ PR("force-constant", awhDimParams.forceConstant());
+ PR("diffusion", awhDimParams.diffusion());
+ PR("cover-diameter", awhDimParams.coverDiameter());
}
-static void pr_awh_bias(FILE* fp, int indent, gmx::AwhBiasParams* awhBiasParams, const char* prefix)
+static void pr_awh_bias(FILE* fp, int indent, const gmx::AwhBiasParams& awhBiasParams, const char* prefix)
{
char opt[STRLEN];
sprintf(opt, "%s-error-init", prefix);
- PR(opt, awhBiasParams->errorInitial);
+ PR(opt, awhBiasParams.initialErrorEstimate());
sprintf(opt, "%s-growth", prefix);
- PS(opt, enumValueToString(awhBiasParams->eGrowth));
+ PS(opt, enumValueToString(awhBiasParams.growthType()));
sprintf(opt, "%s-target", prefix);
- PS(opt, enumValueToString(awhBiasParams->eTarget));
+ PS(opt, enumValueToString(awhBiasParams.targetDistribution()));
sprintf(opt, "%s-target-beta-scalng", prefix);
- PR(opt, awhBiasParams->targetBetaScaling);
+ PR(opt, awhBiasParams.targetBetaScaling());
sprintf(opt, "%s-target-cutoff", prefix);
- PR(opt, awhBiasParams->targetCutoff);
+ PR(opt, awhBiasParams.targetCutoff());
sprintf(opt, "%s-user-data", prefix);
- PS(opt, EBOOL(awhBiasParams->bUserData));
+ PS(opt, EBOOL(awhBiasParams.userPMFEstimate()));
sprintf(opt, "%s-share-group", prefix);
- PI(opt, awhBiasParams->shareGroup);
+ PI(opt, awhBiasParams.shareGroup());
sprintf(opt, "%s-equilibrate-histogram", prefix);
- PS(opt, EBOOL(awhBiasParams->equilibrateHistogram));
+ PS(opt, EBOOL(awhBiasParams.equilibrateHistogram()));
sprintf(opt, "%s-ndim", prefix);
- PI(opt, awhBiasParams->ndim);
+ PI(opt, awhBiasParams.ndim());
- for (int d = 0; d < awhBiasParams->ndim; d++)
+ int d = 0;
+ for (const auto& dimParam : awhBiasParams.dimParams())
{
char prefixdim[STRLEN];
sprintf(prefixdim, "%s-dim%d", prefix, d + 1);
- pr_awh_bias_dim(fp, indent, &awhBiasParams->dimParams[d], prefixdim);
+ pr_awh_bias_dim(fp, indent, dimParam, prefixdim);
+ d++;
}
}
static void pr_awh(FILE* fp, int indent, gmx::AwhParams* awhParams)
{
- PS("awh-potential", enumValueToString(awhParams->ePotential));
- PI("awh-seed", awhParams->seed);
- PI("awh-nstout", awhParams->nstOut);
- PI("awh-nstsample", awhParams->nstSampleCoord);
- PI("awh-nsamples-update", awhParams->numSamplesUpdateFreeEnergy);
- PS("awh-share-bias-multisim", EBOOL(awhParams->shareBiasMultisim));
- PI("awh-nbias", awhParams->numBias);
+ PS("awh-potential", enumValueToString(awhParams->potential()));
+ PI("awh-seed", awhParams->seed());
+ PI("awh-nstout", awhParams->nstout());
+ PI("awh-nstsample", awhParams->nstSampleCoord());
+ PI("awh-nsamples-update", awhParams->numSamplesUpdateFreeEnergy());
+ PS("awh-share-bias-multisim", EBOOL(awhParams->shareBiasMultisim()));
+ PI("awh-nbias", awhParams->numBias());
- for (int k = 0; k < awhParams->numBias; k++)
+ int k = 0;
+ for (const auto& awhBiasParam : awhParams->awhBiasParams())
{
auto prefix = gmx::formatString("awh%d", k + 1);
- pr_awh_bias(fp, indent, &awhParams->awhBiasParams[k], prefix.c_str());
+ pr_awh_bias(fp, indent, awhBiasParam, prefix.c_str());
+ k++;
}
}
PS("awh", EBOOL(ir->bDoAwh));
if (ir->bDoAwh)
{
- pr_awh(fp, indent, ir->awhParams);
+ pr_awh(fp, indent, ir->awhParams.get());
}
/* ENFORCED ROTATION */
}
static void cmp_awhDimParams(FILE* fp,
- const gmx::AwhDimParams* dimp1,
- const gmx::AwhDimParams* dimp2,
+ const gmx::AwhDimParams& dimp1,
+ const gmx::AwhDimParams& dimp2,
int dimIndex,
real ftol,
real abstol)
/* Note that we have double index here, but the compare functions only
* support one index, so here we only print the dim index and not the bias.
*/
- cmp_int(fp, "inputrec.awhParams->bias?->dim->coord_index", dimIndex, dimp1->coordIndex, dimp2->coordIndex);
- cmp_double(fp, "inputrec->awhParams->bias?->dim->period", dimIndex, dimp1->period, dimp2->period, ftol, abstol);
- cmp_double(fp, "inputrec->awhParams->bias?->dim->diffusion", dimIndex, dimp1->diffusion, dimp2->diffusion, ftol, abstol);
- cmp_double(fp, "inputrec->awhParams->bias?->dim->origin", dimIndex, dimp1->origin, dimp2->origin, ftol, abstol);
- cmp_double(fp, "inputrec->awhParams->bias?->dim->end", dimIndex, dimp1->end, dimp2->end, ftol, abstol);
+ cmp_int(fp,
+ "inputrec.awhParams->bias?->dim->coord_index",
+ dimIndex,
+ dimp1.coordinateIndex(),
+ dimp2.coordinateIndex());
+ cmp_double(fp, "inputrec->awhParams->bias?->dim->period", dimIndex, dimp1.period(), dimp2.period(), ftol, abstol);
+ cmp_double(fp,
+ "inputrec->awhParams->bias?->dim->diffusion",
+ dimIndex,
+ dimp1.diffusion(),
+ dimp2.diffusion(),
+ ftol,
+ abstol);
+ cmp_double(fp, "inputrec->awhParams->bias?->dim->origin", dimIndex, dimp1.origin(), dimp2.origin(), ftol, abstol);
+ cmp_double(fp, "inputrec->awhParams->bias?->dim->end", dimIndex, dimp1.end(), dimp2.end(), ftol, abstol);
cmp_double(fp,
"inputrec->awhParams->bias?->dim->coord_value_init",
dimIndex,
- dimp1->coordValueInit,
- dimp2->coordValueInit,
+ dimp1.initialCoordinate(),
+ dimp2.initialCoordinate(),
ftol,
abstol);
cmp_double(fp,
"inputrec->awhParams->bias?->dim->coverDiameter",
dimIndex,
- dimp1->coverDiameter,
- dimp2->coverDiameter,
+ dimp1.coverDiameter(),
+ dimp2.coverDiameter(),
ftol,
abstol);
}
static void cmp_awhBiasParams(FILE* fp,
- const gmx::AwhBiasParams* bias1,
- const gmx::AwhBiasParams* bias2,
+ const gmx::AwhBiasParams& bias1,
+ const gmx::AwhBiasParams& bias2,
int biasIndex,
real ftol,
real abstol)
{
- cmp_int(fp, "inputrec->awhParams->ndim", biasIndex, bias1->ndim, bias2->ndim);
- cmpEnum<gmx::AwhTargetType>(fp, "inputrec->awhParams->biaseTarget", bias1->eTarget, bias2->eTarget);
+ cmp_int(fp, "inputrec->awhParams->ndim", biasIndex, bias1.ndim(), bias2.ndim());
+ cmpEnum<gmx::AwhTargetType>(
+ fp, "inputrec->awhParams->biaseTarget", bias1.targetDistribution(), bias2.targetDistribution());
cmp_double(fp,
"inputrec->awhParams->biastargetBetaScaling",
biasIndex,
- bias1->targetBetaScaling,
- bias2->targetBetaScaling,
+ bias1.targetBetaScaling(),
+ bias2.targetBetaScaling(),
ftol,
abstol);
cmp_double(fp,
"inputrec->awhParams->biastargetCutoff",
biasIndex,
- bias1->targetCutoff,
- bias2->targetCutoff,
+ bias1.targetCutoff(),
+ bias2.targetCutoff(),
ftol,
abstol);
cmpEnum<gmx::AwhHistogramGrowthType>(
- fp, "inputrec->awhParams->biaseGrowth", bias1->eGrowth, bias2->eGrowth);
- cmp_bool(fp, "inputrec->awhParams->biasbUserData", biasIndex, bias1->bUserData, bias2->bUserData);
+ fp, "inputrec->awhParams->biaseGrowth", bias1.growthType(), bias2.growthType());
+ cmp_bool(fp, "inputrec->awhParams->biasbUserData", biasIndex, bias1.userPMFEstimate(), bias2.userPMFEstimate());
cmp_double(fp,
"inputrec->awhParams->biaserror_initial",
biasIndex,
- bias1->errorInitial,
- bias2->errorInitial,
+ bias1.initialErrorEstimate(),
+ bias2.initialErrorEstimate(),
ftol,
abstol);
- cmp_int(fp, "inputrec->awhParams->biasShareGroup", biasIndex, bias1->shareGroup, bias2->shareGroup);
+ cmp_int(fp, "inputrec->awhParams->biasShareGroup", biasIndex, bias1.shareGroup(), bias2.shareGroup());
- for (int dim = 0; dim < std::min(bias1->ndim, bias2->ndim); dim++)
+ const auto dimParams1 = bias1.dimParams();
+ const auto dimParams2 = bias2.dimParams();
+ for (int dim = 0; dim < std::min(bias1.ndim(), bias2.ndim()); dim++)
{
- cmp_awhDimParams(fp, &bias1->dimParams[dim], &bias2->dimParams[dim], dim, ftol, abstol);
+ cmp_awhDimParams(fp, dimParams1[dim], dimParams2[dim], dim, ftol, abstol);
}
}
-static void cmp_awhParams(FILE* fp, const gmx::AwhParams* awh1, const gmx::AwhParams* awh2, real ftol, real abstol)
+static void cmp_awhParams(FILE* fp, const gmx::AwhParams& awh1, const gmx::AwhParams& awh2, real ftol, real abstol)
{
- cmp_int(fp, "inputrec->awhParams->nbias", -1, awh1->numBias, awh2->numBias);
- cmp_int64(fp, "inputrec->awhParams->seed", awh1->seed, awh2->seed);
- cmp_int(fp, "inputrec->awhParams->nstout", -1, awh1->nstOut, awh2->nstOut);
- cmp_int(fp, "inputrec->awhParams->nstsample_coord", -1, awh1->nstSampleCoord, awh2->nstSampleCoord);
+ cmp_int(fp, "inputrec->awhParams->nbias", -1, awh1.numBias(), awh2.numBias());
+ cmp_int64(fp, "inputrec->awhParams->seed", awh1.seed(), awh2.seed());
+ cmp_int(fp, "inputrec->awhParams->nstout", -1, awh1.nstout(), awh2.nstout());
+ cmp_int(fp, "inputrec->awhParams->nstsample_coord", -1, awh1.nstSampleCoord(), awh2.nstSampleCoord());
cmp_int(fp,
"inputrec->awhParams->nsamples_update_free_energy",
-1,
- awh1->numSamplesUpdateFreeEnergy,
- awh2->numSamplesUpdateFreeEnergy);
- cmpEnum<gmx::AwhPotentialType>(fp, "inputrec->awhParams->ePotential", awh1->ePotential, awh2->ePotential);
- cmp_bool(fp, "inputrec->awhParams->shareBiasMultisim", -1, awh1->shareBiasMultisim, awh2->shareBiasMultisim);
+ awh1.numSamplesUpdateFreeEnergy(),
+ awh2.numSamplesUpdateFreeEnergy());
+ cmpEnum<gmx::AwhPotentialType>(
+ fp, "inputrec->awhParams->ePotential", awh1.potential(), awh2.potential());
+ cmp_bool(fp, "inputrec->awhParams->shareBiasMultisim", -1, awh1.shareBiasMultisim(), awh2.shareBiasMultisim());
- if (awh1->numBias == awh2->numBias)
+ if (awh1.numBias() == awh2.numBias())
{
- for (int bias = 0; bias < awh1->numBias; bias++)
+ const auto awhBiasParams1 = awh1.awhBiasParams();
+ const auto awhBiasParams2 = awh2.awhBiasParams();
+ for (int bias = 0; bias < awh1.numBias(); bias++)
{
- cmp_awhBiasParams(fp, &awh1->awhBiasParams[bias], &awh2->awhBiasParams[bias], bias, ftol, abstol);
+ cmp_awhBiasParams(fp, awhBiasParams1[bias], awhBiasParams2[bias], bias, ftol, abstol);
}
}
}
cmp_bool(fp, "inputrec->bDoAwh", -1, ir1->bDoAwh, ir2->bDoAwh);
if (ir1->bDoAwh && ir2->bDoAwh)
{
- cmp_awhParams(fp, ir1->awhParams, ir2->awhParams, ftol, abstol);
+ cmp_awhParams(fp, *ir1->awhParams, *ir2->awhParams, ftol, abstol);
}
cmpEnum(fp, "inputrec->eDisre", ir1->eDisre, ir2->eDisre);
namespace gmx
{
class Awh;
-struct AwhParams;
+class AwhParams;
class KeyValueTreeObject;
struct MtsLevel;
} // namespace gmx
//! Whether to use AWH biasing for PMF calculations
gmx_bool bDoAwh;
//! AWH biasing parameters
- gmx::AwhParams* awhParams;
+ std::unique_ptr<gmx::AwhParams> awhParams;
/* Enforced rotation data */
//! Whether to calculate enforced rotation potential(s)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff