From d4610e91e7815756f728e68785b58b884aead153 Mon Sep 17 00:00:00 2001 From: Paul Bauer Date: Wed, 10 Mar 2021 20:06:31 +0000 Subject: [PATCH] Make AWH parameters proper C++ 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. --- src/gromacs/applied_forces/awh/awh.cpp | 53 +- src/gromacs/applied_forces/awh/awh.h | 4 +- src/gromacs/applied_forces/awh/bias.cpp | 7 +- src/gromacs/applied_forces/awh/bias.h | 4 +- src/gromacs/applied_forces/awh/biasgrid.cpp | 38 +- src/gromacs/applied_forces/awh/biasgrid.h | 4 +- src/gromacs/applied_forces/awh/biasparams.cpp | 44 +- src/gromacs/applied_forces/awh/biasparams.h | 4 +- .../applied_forces/awh/biassharing.cpp | 24 +- src/gromacs/applied_forces/awh/biassharing.h | 2 +- src/gromacs/applied_forces/awh/biasstate.cpp | 2 +- src/gromacs/applied_forces/awh/biasstate.h | 1 - src/gromacs/applied_forces/awh/coordstate.cpp | 7 +- src/gromacs/applied_forces/awh/coordstate.h | 2 +- .../applied_forces/awh/histogramsize.cpp | 4 +- .../applied_forces/awh/histogramsize.h | 2 +- .../applied_forces/awh/read_params.cpp | 850 ++++++++++-------- src/gromacs/applied_forces/awh/read_params.h | 16 +- .../applied_forces/awh/tests/CMakeLists.txt | 3 +- .../applied_forces/awh/tests/awh_setup.cpp | 322 +++++++ .../applied_forces/awh/tests/awh_setup.h | 106 +++ src/gromacs/applied_forces/awh/tests/bias.cpp | 132 +-- .../awh/tests/bias_fep_lambda_state.cpp | 143 +-- .../applied_forces/awh/tests/biasgrid.cpp | 38 +- .../applied_forces/awh/tests/biasstate.cpp | 106 +-- src/gromacs/fileio/tpxio.cpp | 76 +- src/gromacs/gmxana/gmx_awh.cpp | 35 +- src/gromacs/gmxpreprocess/grompp.cpp | 2 +- src/gromacs/gmxpreprocess/readir.cpp | 6 +- src/gromacs/mdrun/md.cpp | 5 +- src/gromacs/mdtypes/awh_params.h | 233 ++++- src/gromacs/mdtypes/inputrec.cpp | 161 ++-- src/gromacs/mdtypes/inputrec.h | 4 +- 33 files changed, 1490 insertions(+), 950 deletions(-) create mode 100644 src/gromacs/applied_forces/awh/tests/awh_setup.cpp create mode 100644 src/gromacs/applied_forces/awh/tests/awh_setup.h diff --git a/src/gromacs/applied_forces/awh/awh.cpp b/src/gromacs/applied_forces/awh/awh.cpp index 176c781954..d987b09f2c 100644 --- a/src/gromacs/applied_forces/awh/awh.cpp +++ b/src/gromacs/applied_forces/awh/awh.cpp @@ -54,6 +54,7 @@ #include #include +#include #include "gromacs/fileio/enxio.h" #include "gromacs/gmxlib/network.h" @@ -115,10 +116,9 @@ struct BiasCoupledToSystem */ 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; } @@ -133,10 +133,9 @@ static bool anyDimUsesPull(const AwhBiasParams& awhBiasParams) */ 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; } @@ -180,8 +179,8 @@ Awh::Awh(FILE* fplog, 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), @@ -210,41 +209,39 @@ Awh::Awh(FILE* fplog, } 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 pullCoordIndex; std::vector 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 { @@ -257,7 +254,7 @@ Awh::Awh(FILE* fplog, (MASTER(commRecord_) ? Bias::ThisRankWillDoIO::Yes : Bias::ThisRankWillDoIO::No); biasCoupledToSystem_.emplace_back(Bias(k, awhParams, - awhParams.awhBiasParams[k], + awhBiasParams[k], dimParams, beta, inputRecord.delta_t, @@ -482,16 +479,14 @@ void Awh::registerAwhWithPull(const AwhParams& awhParams, pull_t* pull_work) { 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()); } } } diff --git a/src/gromacs/applied_forces/awh/awh.h b/src/gromacs/applied_forces/awh/awh.h index 30b0a25b1a..c1d3632b02 100644 --- a/src/gromacs/applied_forces/awh/awh.h +++ b/src/gromacs/applied_forces/awh/awh.h @@ -1,7 +1,7 @@ /* * 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. @@ -88,7 +88,7 @@ namespace gmx template class ArrayRef; struct AwhHistory; -struct AwhParams; +class AwhParams; class Bias; struct BiasCoupledToSystem; class ForceWithVirial; diff --git a/src/gromacs/applied_forces/awh/bias.cpp b/src/gromacs/applied_forces/awh/bias.cpp index 4d1ae3ce84..afb2ae1ecf 100644 --- a/src/gromacs/applied_forces/awh/bias.cpp +++ b/src/gromacs/applied_forces/awh/bias.cpp @@ -370,7 +370,7 @@ Bias::Bias(int biasIndexInCollection, ThisRankWillDoIO thisRankWillDoIO, BiasParams::DisableUpdateSkips disableUpdateSkips) : dimParams_(dimParamsInit.begin(), dimParamsInit.end()), - grid_(dimParamsInit, awhBiasParams.dimParams), + grid_(dimParamsInit, awhBiasParams.dimParams()), params_(awhParams, awhBiasParams, dimParams_, @@ -390,7 +390,8 @@ Bias::Bias(int biasIndexInCollection, /* 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_) { @@ -406,7 +407,7 @@ Bias::Bias(int biasIndexInCollection, ndim(), blockLength, CorrelationGrid::BlockLengthMeasure::Time, - awhParams.nstSampleCoord * mdTimeStep); + awhParams.nstSampleCoord() * mdTimeStep); writer_ = std::make_unique(*this); } diff --git a/src/gromacs/applied_forces/awh/bias.h b/src/gromacs/applied_forces/awh/bias.h index f3386facf3..39eeb991b4 100644 --- a/src/gromacs/applied_forces/awh/bias.h +++ b/src/gromacs/applied_forces/awh/bias.h @@ -73,9 +73,9 @@ namespace gmx { struct AwhBiasHistory; -struct AwhBiasParams; +class AwhBiasParams; struct AwhHistory; -struct AwhParams; +class AwhParams; struct AwhPointStateHistory; class CorrelationGrid; class BiasGrid; diff --git a/src/gromacs/applied_forces/awh/biasgrid.cpp b/src/gromacs/applied_forces/awh/biasgrid.cpp index 1f4ce66d5f..bd3fe10c58 100644 --- a/src/gromacs/applied_forces/awh/biasgrid.cpp +++ b/src/gromacs/applied_forces/awh/biasgrid.cpp @@ -69,29 +69,33 @@ namespace { /*! \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 @@ -180,7 +184,7 @@ double getDeviationPeriodic(double x, double x0, double period) if (period > 0) { - centerPeriodicValueAroundZero(&dev, period); + dev = centerPeriodicValueAroundZero(dev, period); } return dev; @@ -737,7 +741,8 @@ void BiasGrid::initPoints() 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]; @@ -811,18 +816,19 @@ GridAxis::GridAxis(double origin, double end, double period, int numPoints, bool } } -BiasGrid::BiasGrid(ArrayRef dimParams, const AwhDimParams* awhDimParams) +BiasGrid::BiasGrid(ArrayRef dimParams, ArrayRef 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 " diff --git a/src/gromacs/applied_forces/awh/biasgrid.h b/src/gromacs/applied_forces/awh/biasgrid.h index 0a4bbc3034..2c6777dcb2 100644 --- a/src/gromacs/applied_forces/awh/biasgrid.h +++ b/src/gromacs/applied_forces/awh/biasgrid.h @@ -65,7 +65,7 @@ namespace gmx { -struct AwhDimParams; +class AwhDimParams; /*! \internal * \brief An axis, i.e. dimension, of the grid. @@ -199,7 +199,7 @@ public: * coordinate living on the grid (determines the grid spacing). * \param[in] awhDimParams Dimension params from inputrec. */ - BiasGrid(ArrayRef dimParams, const AwhDimParams* awhDimParams); + BiasGrid(ArrayRef dimParams, ArrayRef awhDimParams); /*! \brief Returns the number of points in the grid. * diff --git a/src/gromacs/applied_forces/awh/biasparams.cpp b/src/gromacs/applied_forces/awh/biasparams.cpp index 9524c7047b..f6d206ee11 100644 --- a/src/gromacs/applied_forces/awh/biasparams.cpp +++ b/src/gromacs/applied_forces/awh/biasparams.cpp @@ -81,7 +81,7 @@ int64_t calcTargetUpdateInterval(const AwhParams& awhParams, const AwhBiasParams * (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: @@ -89,9 +89,9 @@ int64_t calcTargetUpdateInterval(const AwhParams& awhParams, const AwhBiasParams /* 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. @@ -100,7 +100,7 @@ int64_t calcTargetUpdateInterval(const AwhParams& awhParams, const AwhBiasParams 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; } @@ -152,11 +152,11 @@ int64_t calcCheckCoveringInterval(const AwhParams& awhParams, /* 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, @@ -182,16 +182,17 @@ double getInitialHistogramSizeEstimate(const AwhBiasParams& awhBiasParams, /* Get diffusion factor */ double maxCrossingTime = 0.; std::vector 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; @@ -209,7 +210,7 @@ int getNumSharedUpdate(const AwhBiasParams& awhBiasParams, int numSharingSimulat int numShared = 1; - if (awhBiasParams.shareGroup > 0) + if (awhBiasParams.shareGroup() > 0) { /* We do not yet support sharing within a simulation */ int numSharedWithinThisSimulation = 1; @@ -231,21 +232,21 @@ BiasParams::BiasParams(const AwhParams& awhParams, ArrayRef 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) { @@ -254,11 +255,12 @@ BiasParams::BiasParams(const AwhParams& awhParams, 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(std::round(coverRadiusInNm / spacing)) : 0; } diff --git a/src/gromacs/applied_forces/awh/biasparams.h b/src/gromacs/applied_forces/awh/biasparams.h index b416d98556..ec72ace256 100644 --- a/src/gromacs/applied_forces/awh/biasparams.h +++ b/src/gromacs/applied_forces/awh/biasparams.h @@ -62,8 +62,8 @@ namespace gmx template class ArrayRef; -struct AwhBiasParams; -struct AwhParams; +class AwhBiasParams; +class AwhParams; struct DimParams; class GridAxis; enum class AwhTargetType : int; diff --git a/src/gromacs/applied_forces/awh/biassharing.cpp b/src/gromacs/applied_forces/awh/biassharing.cpp index 73bc4a6e8d..50d2603778 100644 --- a/src/gromacs/applied_forces/awh/biassharing.cpp +++ b/src/gromacs/applied_forces/awh/biassharing.cpp @@ -51,6 +51,7 @@ #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" @@ -61,14 +62,16 @@ bool haveBiasSharingWithinSimulation(const AwhParams& awhParams) { 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; } @@ -90,9 +93,9 @@ void biasesAreCompatibleForSharingBetweenSimulations(const AwhParams& awhP * 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++; @@ -117,8 +120,8 @@ void biasesAreCompatibleForSharingBetweenSimulations(const AwhParams& awhP } std::vector 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++) { @@ -137,9 +140,10 @@ void biasesAreCompatibleForSharingBetweenSimulations(const AwhParams& awhP /* 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 pointSizes(numSim); pointSizes[multiSimComm->simulationIndex_] = pointSize[b]; diff --git a/src/gromacs/applied_forces/awh/biassharing.h b/src/gromacs/applied_forces/awh/biassharing.h index 0b79170937..1d06d0ec69 100644 --- a/src/gromacs/applied_forces/awh/biassharing.h +++ b/src/gromacs/applied_forces/awh/biassharing.h @@ -58,7 +58,7 @@ namespace gmx template class ArrayRef; -struct AwhParams; +class AwhParams; /*! \brief Returns if any bias is sharing within a simulation. * diff --git a/src/gromacs/applied_forces/awh/biasstate.cpp b/src/gromacs/applied_forces/awh/biasstate.cpp index d7c35b2351..1206622f1d 100644 --- a/src/gromacs/applied_forces/awh/biasstate.cpp +++ b/src/gromacs/applied_forces/awh/biasstate.cpp @@ -1831,7 +1831,7 @@ void BiasState::initGridPointState(const AwhBiasParams& awhBiasParams, /* 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); diff --git a/src/gromacs/applied_forces/awh/biasstate.h b/src/gromacs/applied_forces/awh/biasstate.h index 97a632ffbd..859f524af6 100644 --- a/src/gromacs/applied_forces/awh/biasstate.h +++ b/src/gromacs/applied_forces/awh/biasstate.h @@ -75,7 +75,6 @@ namespace gmx template class ArrayRef; struct AwhBiasHistory; -struct AwhBiasParams; class BiasParams; class BiasGrid; class GridAxis; diff --git a/src/gromacs/applied_forces/awh/coordstate.cpp b/src/gromacs/applied_forces/awh/coordstate.cpp index 544a940bb4..f29eee36b0 100644 --- a/src/gromacs/applied_forces/awh/coordstate.cpp +++ b/src/gromacs/applied_forces/awh/coordstate.cpp @@ -67,9 +67,12 @@ CoordState::CoordState(const AwhBiasParams& awhBiasParams, ArrayRef 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. diff --git a/src/gromacs/applied_forces/awh/coordstate.h b/src/gromacs/applied_forces/awh/coordstate.h index 05995f826c..d9dff8fa68 100644 --- a/src/gromacs/applied_forces/awh/coordstate.h +++ b/src/gromacs/applied_forces/awh/coordstate.h @@ -61,7 +61,7 @@ namespace gmx template class ArrayRef; -struct AwhBiasParams; +class AwhBiasParams; struct AwhBiasStateHistory; class BiasParams; class BiasGrid; diff --git a/src/gromacs/applied_forces/awh/histogramsize.cpp b/src/gromacs/applied_forces/awh/histogramsize.cpp index 059d8cc6de..15d1ace4ba 100644 --- a/src/gromacs/applied_forces/awh/histogramsize.cpp +++ b/src/gromacs/applied_forces/awh/histogramsize.cpp @@ -67,8 +67,8 @@ namespace gmx 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) diff --git a/src/gromacs/applied_forces/awh/histogramsize.h b/src/gromacs/applied_forces/awh/histogramsize.h index 9c60adedfe..f0c5725d33 100644 --- a/src/gromacs/applied_forces/awh/histogramsize.h +++ b/src/gromacs/applied_forces/awh/histogramsize.h @@ -64,7 +64,7 @@ namespace gmx template class ArrayRef; struct AwhBiasStateHistory; -struct AwhBiasParams; +class AwhBiasParams; class BiasParams; class PointState; diff --git a/src/gromacs/applied_forces/awh/read_params.cpp b/src/gromacs/applied_forces/awh/read_params.cpp index a889c0e715..c2b8e33d7d 100644 --- a/src/gromacs/applied_forces/awh/read_params.cpp +++ b/src/gromacs/applied_forces/awh/read_params.cpp @@ -57,6 +57,7 @@ #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" @@ -118,12 +119,11 @@ void checkMtsConsistency(const t_inputrec& inputrec, warninp_t wi) 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; @@ -145,7 +145,7 @@ void checkMtsConsistency(const t_inputrec& inputrec, warninp_t wi) "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"); @@ -161,96 +161,96 @@ void checkMtsConsistency(const t_inputrec& inputrec, warninp_t wi) * \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)); } } @@ -266,7 +266,7 @@ void checkPullDimParams(const std::string& prefix, * \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) @@ -303,39 +303,39 @@ void checkFepLambdaDimParams(const std::string& prefix, "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, @@ -368,78 +368,6 @@ void checkFepLambdaDimDecouplingConsistency(const gmx_mtop_t& mtop, warninp_t 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* 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(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. * @@ -448,27 +376,27 @@ void readDimParams(std::vector* inp, * \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 { @@ -478,126 +406,6 @@ void checkDimParams(const std::string& prefix, AwhDimParams* dimParams, const t_ } } -/*! \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* 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(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(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(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(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. * @@ -606,16 +414,17 @@ void readBiasParams(std::vector* inp, * \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'.", @@ -624,8 +433,8 @@ void checkBiasParams(const AwhBiasParams* awhBiasParams, const std::string& pref 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 " @@ -638,84 +447,102 @@ void checkBiasParams(const AwhBiasParams* awhBiasParams, const std::string& pref } 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 @@ -730,36 +557,38 @@ void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi) * 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, @@ -767,7 +596,7 @@ void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi) "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, @@ -779,7 +608,7 @@ void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi) } } - if (awhParams.shareBiasMultisim && !haveSharedBias) + if (awhParams.shareBiasMultisim() && !haveSharedBias) { warning(wi, "Sharing of biases over multiple simulations is requested, but no bias is marked " @@ -794,110 +623,365 @@ void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi) "this (yet)"); } } + } // namespace -AwhParams* readAwhParams(std::vector* inp, warninp_t wi) +AwhDimParams::AwhDimParams(std::vector* 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(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* 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(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(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(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(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(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(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* 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(inp, opt.c_str(), wi); + opt = "awh-potential"; + potentialEnum_ = getEnum(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(gmx::makeRandomSeed()); - fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed); + seed_ = static_cast(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(inp, opt.c_str(), wi) != Boolean::No); + opt = "awh-share-multisim"; + shareBiasMultisim_ = (getEnum(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); + } } } @@ -1005,17 +1089,18 @@ static bool valueIsInInterval(double origin, double end, double period, double v * \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)) { @@ -1099,13 +1184,13 @@ static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t 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) { @@ -1118,9 +1203,10 @@ static void setStateDependentAwhPullDimParams(AwhDimParams* dimParams, 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)) { @@ -1136,7 +1222,7 @@ static void setStateDependentAwhPullDimParams(AwhDimParams* dimParams, { 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, @@ -1146,9 +1232,9 @@ static void setStateDependentAwhPullDimParams(AwhDimParams* dimParams, } /* 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, @@ -1183,29 +1269,69 @@ 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 diff --git a/src/gromacs/applied_forces/awh/read_params.h b/src/gromacs/applied_forces/awh/read_params.h index bb1cd63fa9..22247ae480 100644 --- a/src/gromacs/applied_forces/awh/read_params.h +++ b/src/gromacs/applied_forces/awh/read_params.h @@ -2,7 +2,7 @@ * 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. @@ -59,23 +59,17 @@ enum class PbcType : int; 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* 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 diff --git a/src/gromacs/applied_forces/awh/tests/CMakeLists.txt b/src/gromacs/applied_forces/awh/tests/CMakeLists.txt index bbce544ace..864dcabab8 100644 --- a/src/gromacs/applied_forces/awh/tests/CMakeLists.txt +++ b/src/gromacs/applied_forces/awh/tests/CMakeLists.txt @@ -1,7 +1,7 @@ # # 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. @@ -34,6 +34,7 @@ gmx_add_unit_test(AwhTest awh-test CPP_SOURCE_FILES + awh_setup.cpp bias.cpp biasgrid.cpp biasstate.cpp diff --git a/src/gromacs/applied_forces/awh/tests/awh_setup.cpp b/src/gromacs/applied_forces/awh/tests/awh_setup.cpp new file mode 100644 index 0000000000..edb67b5e51 --- /dev/null +++ b/src/gromacs/applied_forces/awh/tests/awh_setup.cpp @@ -0,0 +1,322 @@ +/* + * 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 + +#include +#include +#include + +#include +#include +#include + +#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 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 awhBiasParamSerialized(AwhHistogramGrowthType eawhgrowth, + double beta, + double inputErrorScaling, + ArrayRef> 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(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 awhParamSerialized(AwhHistogramGrowthType eawhgrowth, + AwhPotentialType eawhpotential, + double beta, + double inputErrorScaling, + int64_t inputSeed, + ArrayRef> 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> 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 diff --git a/src/gromacs/applied_forces/awh/tests/awh_setup.h b/src/gromacs/applied_forces/awh/tests/awh_setup.h new file mode 100644 index 0000000000..18bd4be1a9 --- /dev/null +++ b/src/gromacs/applied_forces/awh/tests/awh_setup.h @@ -0,0 +1,106 @@ +/* + * 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 +#include + +#include "gromacs/applied_forces/awh/bias.h" +#include "gromacs/mdtypes/awh_params.h" + +namespace gmx +{ + +template +class ArrayRef; + +namespace test +{ + +/*! \internal \brief + * Prepare a memory buffer with serialized AwhDimParams. + */ +std::vector 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; +}; + +/*! \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> dimensionParameterBuffers, + bool inputUserData, + double beta, + bool useAwhFep, + double inputErrorScaling, + int numFepLambdaStates); + +} // namespace test +} // namespace gmx + +#endif diff --git a/src/gromacs/applied_forces/awh/tests/bias.cpp b/src/gromacs/applied_forces/awh/tests/bias.cpp index 5abe89de9d..85e033ef9f 100644 --- a/src/gromacs/applied_forces/awh/tests/bias.cpp +++ b/src/gromacs/applied_forces/awh/tests/bias.cpp @@ -43,13 +43,16 @@ #include #include +#include #include #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" @@ -59,86 +62,9 @@ namespace gmx 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; //!< 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 BiasTestParameters; @@ -147,10 +73,13 @@ typedef std::tuple { +private: + //! Storage for test parameters. + std::unique_ptr params_; + public: //! Random seed for AWH MC sampling int64_t seed_; - //! Coordinates representing a trajectory in time std::vector coordinates_; //! The awh Bias @@ -183,24 +112,27 @@ public: * 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(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(-1, - params.awhParams, - params.awhBiasParams, - params.dimParams, - params.beta, + params_->awhParams, + params_->awhParams.awhBiasParams()[0], + params_->dimParams, + params_->beta, mdTimeStep, 1, "", @@ -297,15 +229,23 @@ INSTANTIATE_TEST_CASE_P(WithParameters, // Test that we detect coverings and exit the initial stage at the correct step TEST(BiasTest, DetectsCovering) { + const std::vector 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, @@ -317,9 +257,9 @@ TEST(BiasTest, DetectsCovering) * 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 */ @@ -341,7 +281,7 @@ TEST(BiasTest, DetectsCovering) nullptr, step, step, - params.awhParams.seed, + params.awhParams.seed(), nullptr); inInitialStage = bias.state().inInitialStage(); diff --git a/src/gromacs/applied_forces/awh/tests/bias_fep_lambda_state.cpp b/src/gromacs/applied_forces/awh/tests/bias_fep_lambda_state.cpp index b53777cfd8..787b88f960 100644 --- a/src/gromacs/applied_forces/awh/tests/bias_fep_lambda_state.cpp +++ b/src/gromacs/applied_forces/awh/tests/bias_fep_lambda_state.cpp @@ -49,6 +49,7 @@ #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" @@ -59,83 +60,8 @@ namespace test { //! 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; //!< 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 BiasTestParameters; @@ -144,6 +70,10 @@ typedef std::tuple { +private: + //! Storage for test parameters. + std::unique_ptr params_; + public: //! Random seed for AWH MC sampling int64_t seed_; @@ -177,18 +107,28 @@ public: /* 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(getAwhTestParameters( + eawhgrowth, eawhpotential, awhDimArrayRef, false, 0.4, true, 1.0, c_numLambdaStates)); + + seed_ = params_->awhParams.seed(); double mdTimeStep = 0.1; bias_ = std::make_unique(-1, - params.awhParams, - params.awhBiasParams, - params.dimParams, - params.beta, + params_->awhParams, + params_->awhParams.awhBiasParams()[0], + params_->dimParams, + params_->beta, mdTimeStep, 1, "", @@ -221,10 +161,10 @@ TEST_P(BiasFepLambdaStateTest, ForcesBiasPmf) int nSteps = 501; /* Some energies to use as base values (to which some noise is added later on). */ - std::vector neighborLambdaEnergies(numLambdaStates); - std::vector neighborLambdaDhdl(numLambdaStates); + std::vector neighborLambdaEnergies(c_numLambdaStates); + std::vector 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); @@ -294,14 +234,29 @@ INSTANTIATE_TEST_CASE_P(WithParameters, // 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, @@ -314,10 +269,10 @@ TEST(BiasFepLambdaStateTest, DetectsCovering) bool inInitialStage = bias.state().inInitialStage(); /* Some energies to use as base values (to which some noise is added later on). */ - std::vector neighborLambdaEnergies(numLambdaStates); - std::vector neighborLambdaDhdl(numLambdaStates); + std::vector neighborLambdaEnergies(c_numLambdaStates); + std::vector 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); @@ -341,7 +296,7 @@ TEST(BiasFepLambdaStateTest, DetectsCovering) nullptr, step, step, - params.awhParams.seed, + params.awhParams.seed(), nullptr); inInitialStage = bias.state().inInitialStage(); diff --git a/src/gromacs/applied_forces/awh/tests/biasgrid.cpp b/src/gromacs/applied_forces/awh/tests/biasgrid.cpp index 9337d0632d..d008eac9ec 100644 --- a/src/gromacs/applied_forces/awh/tests/biasgrid.cpp +++ b/src/gromacs/applied_forces/awh/tests/biasgrid.cpp @@ -1,7 +1,7 @@ /* * 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. @@ -46,7 +46,9 @@ #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 @@ -71,15 +73,29 @@ TEST(biasGridTest, neighborhood) */ const int numDim = 2; - std::vector 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; + 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; @@ -89,7 +105,7 @@ TEST(biasGridTest, neighborhood) 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(); diff --git a/src/gromacs/applied_forces/awh/tests/biasstate.cpp b/src/gromacs/applied_forces/awh/tests/biasstate.cpp index c78d55061a..fdffe31cfe 100644 --- a/src/gromacs/applied_forces/awh/tests/biasstate.cpp +++ b/src/gromacs/applied_forces/awh/tests/biasstate.cpp @@ -51,6 +51,7 @@ #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" @@ -60,87 +61,45 @@ namespace gmx 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 { +private: + std::unique_ptr params_; + public: std::unique_ptr biasState_; //!< The bias state BiasStateTest() { - AwhTestParameters params = getAwhTestParameters(); - const AwhParams& awhParams = params.awhParams; - const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[0]; + std::vector> 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(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.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(awhBiasParams, 1.0, dimParams, grid); @@ -148,10 +107,7 @@ public: // 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()); } }; diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index 2336b231ff..bcaaca2b88 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -49,6 +49,7 @@ #include #include +#include "gromacs/applied_forces/awh/read_params.h" #include "gromacs/fileio/filetypes.h" #include "gromacs/fileio/gmxfio.h" #include "gromacs/fileio/gmxfio_xdr.h" @@ -75,6 +76,7 @@ #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" @@ -655,73 +657,6 @@ static void do_fepvals(gmx::ISerializer* serializer, t_lambda* fepvals, int file } } -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(temp); - } - else - { - int temp = static_cast(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; @@ -1568,9 +1503,12 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v { if (serializer->reading()) { - snew(ir->awhParams, 1); + ir->awhParams = std::make_unique(serializer); + } + else + { + ir->awhParams->serialize(serializer); } - do_awh(serializer, ir->awhParams, file_version); } } else diff --git a/src/gromacs/gmxana/gmx_awh.cpp b/src/gromacs/gmxana/gmx_awh.cpp index 690d932eea..f83e9ba3a6 100644 --- a/src/gromacs/gmxana/gmx_awh.cpp +++ b/src/gromacs/gmxana/gmx_awh.cpp @@ -52,6 +52,7 @@ #include #include +#include "gromacs/applied_forces/awh/read_params.h" #include "gromacs/commandline/pargs.h" #include "gromacs/fileio/enxio.h" #include "gromacs/fileio/gmxfio.h" @@ -122,7 +123,7 @@ public: */ void initializeAwhOutputFile(int subBlockStart, int numSubBlocks, - const AwhBiasParams* awhBiasParams, + const AwhBiasParams& awhBiasParams, AwhGraphSelection graphSelection, EnergyUnit energyUnit, real kTValue); @@ -137,7 +138,7 @@ public: */ void initializeFrictionOutputFile(int subBlockStart, int numSubBlocks, - const AwhBiasParams* awhBiasParams, + const AwhBiasParams& awhBiasParams, EnergyUnit energyUnit, real kTValue); @@ -210,7 +211,7 @@ class AwhReader { public: //! Constructor - AwhReader(const AwhParams* awhParams, + AwhReader(const AwhParams& awhParams, int numFileOptions, const t_filenm* filenames, AwhGraphSelection awhGraphSelection, @@ -246,7 +247,7 @@ enum class OutputFileType constexpr int maxAwhGraphs = 6; /*! \brief Constructs a legend for a standard awh output file */ -std::vector makeLegend(const AwhBiasParams* awhBiasParams, +std::vector makeLegend(const AwhBiasParams& awhBiasParams, OutputFileType outputFileType, size_t numLegend) { @@ -259,7 +260,7 @@ std::vector makeLegend(const AwhBiasParams* awhBiasParams, std::vector 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)); } @@ -278,7 +279,7 @@ std::vector makeLegend(const AwhBiasParams* awhBiasParams, } 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++) { @@ -315,13 +316,13 @@ OutputFile::OutputFile(const std::string& filename, const std::string& baseTitle 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) { @@ -347,19 +348,19 @@ void OutputFile::initializeAwhOutputFile(int subblockStart, 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 */ @@ -392,7 +393,7 @@ void OutputFile::initializeFrictionOutputFile(int subBlockStart } } -AwhReader::AwhReader(const AwhParams* awhParams, +AwhReader::AwhReader(const AwhParams& awhParams, int numFileOptions, const t_filenm* filenames, AwhGraphSelection awhGraphSelection, @@ -413,14 +414,14 @@ AwhReader::AwhReader(const AwhParams* awhParams, /* 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(block->sub[subblockStart].fval[0]); std::unique_ptr 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); @@ -429,7 +430,7 @@ AwhReader::AwhReader(const AwhParams* awhParams, if (outputFriction) { frictionOutputFile = std::make_unique( - 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); @@ -630,7 +631,7 @@ int gmx_awh(int argc, char* argv[]) AwhGraphSelection awhGraphSelection = (moreGraphs ? AwhGraphSelection::All : AwhGraphSelection::Pmf); EnergyUnit energyUnit = (kTUnit ? EnergyUnit::KT : EnergyUnit::KJPerMol); - awhReader = std::make_unique(ir.awhParams, + awhReader = std::make_unique(*ir.awhParams, nfile, fnm, awhGraphSelection, diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 67abdfc7ec..9bd9996437 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -2412,7 +2412,7 @@ int gmx_grompp(int argc, char* argv[]) copy_mat(ir->compress, compressibility); } setStateDependentAwhParams( - ir->awhParams, + ir->awhParams.get(), *ir->pull, pull, state.box, diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index eaf605d85b..fba69a3e46 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -37,7 +37,6 @@ */ #include "gmxpre.h" -#include "gromacs/utility/enumerationhelpers.h" #include "readir.h" #include @@ -61,6 +60,7 @@ #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" @@ -2223,7 +2223,7 @@ void get_ir(const char* mdparin, ir->bDoAwh = (getEnum(&inp, "awh", wi) != Boolean::No); if (ir->bDoAwh) { - ir->awhParams = gmx::readAwhParams(&inp, wi); + ir->awhParams = std::make_unique(&inp, *ir, wi); } /* Enforced rotation */ @@ -2784,7 +2784,7 @@ void get_ir(const char* mdparin, if (ir->bDoAwh) { - gmx::checkAwhParams(ir->awhParams, ir, wi); + gmx::checkAwhParams(*ir->awhParams, *ir, wi); } sfree(dumstr[0]); diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 3c8663e92b..cb086641f4 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -53,6 +53,7 @@ #include #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" @@ -285,7 +286,7 @@ void gmx::LegacySimulator::do_md() // 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 @@ -599,7 +600,7 @@ void gmx::LegacySimulator::do_md() } if (ir->bDoAwh) { - nstfep = std::gcd(ir->awhParams->nstSampleCoord, nstfep); + nstfep = std::gcd(ir->awhParams->nstSampleCoord(), nstfep); } } diff --git a/src/gromacs/mdtypes/awh_params.h b/src/gromacs/mdtypes/awh_params.h index debbb84672..729707e0e0 100644 --- a/src/gromacs/mdtypes/awh_params.h +++ b/src/gromacs/mdtypes/awh_params.h @@ -51,12 +51,22 @@ #ifndef GMX_MDTYPES_AWH_PARAMS_H #define GMX_MDTYPES_AWH_PARAMS_H +#include + #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 { @@ -103,53 +113,196 @@ enum class AwhCoordinateProviderType : 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* 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* 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 dimParams() { return dimParams_; } + //! Const access to dimension parameters. + ArrayRef dimParams() const { return dimParams_; } + //! Write datastructure. + void serialize(ISerializer* serializer); + +private: + //! AWH parameters per dimension. + std::vector 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* 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() { return awhBiasParams_; } + //! Const access to bias parameters. + ArrayRef 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_; + //! 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 diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp index b0853f0de3..ab292ec523 100644 --- a/src/gromacs/mdtypes/inputrec.cpp +++ b/src/gromacs/mdtypes/inputrec.cpp @@ -48,6 +48,7 @@ #include #include +#include "gromacs/applied_forces/awh/read_params.h" #include "gromacs/math/veccompare.h" #include "gromacs/math/vecdump.h" #include "gromacs/mdtypes/awh_params.h" @@ -617,66 +618,70 @@ static void pr_pull(FILE* fp, int indent, const pull_params_t& pull) } } -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++; } } @@ -970,7 +975,7 @@ void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir, PS("awh", EBOOL(ir->bDoAwh)); if (ir->bDoAwh) { - pr_awh(fp, indent, ir->awhParams); + pr_awh(fp, indent, ir->awhParams.get()); } /* ENFORCED ROTATION */ @@ -1103,8 +1108,8 @@ static void cmp_pull(FILE* fp) } 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) @@ -1112,87 +1117,103 @@ static void cmp_awhDimParams(FILE* fp, /* 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(fp, "inputrec->awhParams->biaseTarget", bias1->eTarget, bias2->eTarget); + cmp_int(fp, "inputrec->awhParams->ndim", biasIndex, bias1.ndim(), bias2.ndim()); + cmpEnum( + 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( - 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(fp, "inputrec->awhParams->ePotential", awh1->ePotential, awh2->ePotential); - cmp_bool(fp, "inputrec->awhParams->shareBiasMultisim", -1, awh1->shareBiasMultisim, awh2->shareBiasMultisim); + awh1.numSamplesUpdateFreeEnergy(), + awh2.numSamplesUpdateFreeEnergy()); + cmpEnum( + 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); } } } @@ -1459,7 +1480,7 @@ void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real f 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); diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h index 4ca33733b7..7b50eb0b9b 100644 --- a/src/gromacs/mdtypes/inputrec.h +++ b/src/gromacs/mdtypes/inputrec.h @@ -59,7 +59,7 @@ struct pull_params_t; namespace gmx { class Awh; -struct AwhParams; +class AwhParams; class KeyValueTreeObject; struct MtsLevel; } // namespace gmx @@ -518,7 +518,7 @@ struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding) //! Whether to use AWH biasing for PMF calculations gmx_bool bDoAwh; //! AWH biasing parameters - gmx::AwhParams* awhParams; + std::unique_ptr awhParams; /* Enforced rotation data */ //! Whether to calculate enforced rotation potential(s) -- 2.22.0