* 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 "read_params.h"
+#include <algorithm>
+
#include "gromacs/applied_forces/awh/awh.h"
#include "gromacs/fileio/readinp.h"
#include "gromacs/fileio/warninp.h"
#include "gromacs/random/seed.h"
#include "gromacs/topology/mtop_util.h"
#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"
namespace gmx
{
-const char* eawhtarget_names[eawhtargetNR + 1] = { "constant", "cutoff", "boltzmann",
- "local-boltzmann", nullptr };
+const char* enumValueToString(AwhTargetType enumValue)
+{
+ constexpr gmx::EnumerationArray<AwhTargetType, const char*> awhTargetTypeNames = {
+ "constant", "cutoff", "boltzmann", "local-boltzmann"
+ };
+ return awhTargetTypeNames[enumValue];
+}
-const char* eawhgrowth_names[eawhgrowthNR + 1] = { "exp-linear", "linear", nullptr };
+const char* enumValueToString(AwhHistogramGrowthType enumValue)
+{
+ constexpr gmx::EnumerationArray<AwhHistogramGrowthType, const char*> awhHistogramGrowthTypeNames = {
+ "exp-linear", "linear"
+ };
+ return awhHistogramGrowthTypeNames[enumValue];
+}
-const char* eawhpotential_names[eawhpotentialNR + 1] = { "convolved", "umbrella", nullptr };
+const char* enumValueToString(AwhPotentialType enumValue)
+{
+ constexpr gmx::EnumerationArray<AwhPotentialType, const char*> awhPotentialTypeNames = {
+ "convolved", "umbrella"
+ };
+ return awhPotentialTypeNames[enumValue];
+}
-const char* eawhcoordprovider_names[eawhcoordproviderNR + 1] = { "pull", "fep-lambda", nullptr };
+const char* enumValueToString(AwhCoordinateProviderType enumValue)
+{
+ constexpr gmx::EnumerationArray<AwhCoordinateProviderType, const char*> awhCoordinateProviderTypeNames = {
+ "pull", "fep-lambda"
+ };
+ return awhCoordinateProviderTypeNames[enumValue];
+}
namespace
{
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 eawhcoordproviderPULL: usesPull = true; break;
- case eawhcoordproviderFREE_ENERGY_LAMBDA: usesFep = true; break;
+ case AwhCoordinateProviderType::Pull: usesPull = true; break;
+ case AwhCoordinateProviderType::FreeEnergyLambda: usesFep = true; break;
default: GMX_RELEASE_ASSERT(false, "Unsupported coord provider");
}
}
"computed at the slow MTS level");
}
- if (inputrec.awhParams->nstSampleCoord % inputrec.mtsLevels[awhMtsLevel].stepFactor != 0)
+ if (inputrec.awhParams->nstSampleCoord() % inputrec.mtsLevels[awhMtsLevel].stepFactor != 0)
{
warning_error(wi,
"With MTS applied to AWH, awh-nstsample should be a multiple of mts-factor");
* \param[in,out] wi Struct for bookeeping warnings.
*/
void checkPullDimParams(const std::string& prefix,
- AwhDimParams* dimParams,
+ const AwhDimParams& dimParams,
const pull_params_t& pull_params,
warninp_t wi)
{
- if (dimParams->coordIndex < 0)
+ if (dimParams.coordinateIndex() < 0)
{
gmx_fatal(FARGS,
"Failed to read a valid coordinate index for %s-coord-index. "
"Note that the pull coordinate indexing starts at 1.",
prefix.c_str());
}
- if (dimParams->coordIndex >= pull_params.ncoord)
+ if (dimParams.coordinateIndex() >= pull_params.ncoord)
{
gmx_fatal(FARGS,
"The given AWH coordinate index (%d) is larger than the number of pull "
"coordinates (%d)",
- dimParams->coordIndex + 1, pull_params.ncoord);
+ 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, prefix.c_str(), dimParams->end);
+ prefix.c_str(),
+ dimParams.origin(),
+ prefix.c_str(),
+ 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 */
- int 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 == epullgDIST)
+ 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, prefix.c_str(), dimParams->end,
- EPULLGEOM(epullgDIR));
+ prefix.c_str(),
+ dimParams.origin(),
+ prefix.c_str(),
+ dimParams.end(),
+ enumValueToString(PullGroupGeometry::Direction));
}
}
- else if (eGeom == epullgANGLE || eGeom == epullgANGLEAXIS)
+ 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, prefix.c_str(), dimParams->end,
- EPULLGEOM(epullgANGLE), EPULLGEOM(epullgANGLEAXIS));
+ prefix.c_str(),
+ dimParams.origin(),
+ prefix.c_str(),
+ dimParams.end(),
+ enumValueToString(PullGroupGeometry::Angle),
+ enumValueToString(PullGroupGeometry::AngleAxis));
}
}
- else if (eGeom == epullgDIHEDRAL)
+ 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, prefix.c_str(), dimParams->end,
- EPULLGEOM(epullgDIHEDRAL));
+ prefix.c_str(),
+ dimParams.origin(),
+ prefix.c_str(),
+ dimParams.end(),
+ enumValueToString(PullGroupGeometry::Dihedral));
}
}
}
* \param[in] efep This is the type of FEP calculation (efep enumerator).
* \param[in,out] wi Struct for bookeeping warnings.
*/
-void checkFepLambdaDimParams(const std::string& prefix,
- const AwhDimParams* dimParams,
- const t_lambda* lambdaParams,
- const int efep,
- warninp_t wi)
+void checkFepLambdaDimParams(const std::string& prefix,
+ const AwhDimParams& dimParams,
+ const t_lambda* lambdaParams,
+ const FreeEnergyPerturbationType efep,
+ warninp_t wi)
{
std::string opt;
"When running AWH coupled to the free energy lambda state all lambda states "
"should be used as neighbors in order to get correct probabilities, i.e. "
"calc-lambda-neighbors (%d) must be %d.",
- lambdaParams->lambda_neighbors, -1);
+ lambdaParams->lambda_neighbors,
+ -1);
}
- if (efep == efepSLOWGROWTH || lambdaParams->delta_lambda != 0)
+ if (efep == FreeEnergyPerturbationType::SlowGrowth || lambdaParams->delta_lambda != 0)
{
gmx_fatal(FARGS,
"AWH coupled to the free energy lambda state is not compatible with slow-growth "
"and delta-lambda must be 0.");
}
- if (efep == efepEXPANDED)
+ if (efep == FreeEnergyPerturbationType::Expanded)
{
gmx_fatal(FARGS,
"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);
+ opt.c_str(),
+ 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, lambdaParams->n_lambda);
+ opt.c_str(),
+ 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, prefix.c_str(), dimParams->end);
+ prefix.c_str(),
+ dimParams.origin(),
+ prefix.c_str(),
+ dimParams.end());
warning(wi, message);
}
- if (dimParams->forceConstant != 0)
+ if (dimParams.forceConstant() != 0)
{
warning_error(
wi,
}
}
-/*! \brief
- * Read parameters of an AWH bias dimension.
- *
- * \param[in,out] inp Input file entries.
- * \param[in] prefix Prefix for dimension parameters.
- * \param[in,out] dimParams AWH dimensional parameters.
- * \param[in,out] wi Struct for bookeeping warnings.
- * \param[in] bComment True if comments should be printed.
- */
-void readDimParams(std::vector<t_inpfile>* inp,
- const std::string& prefix,
- AwhDimParams* dimParams,
- warninp_t wi,
- bool bComment)
-{
- std::string opt;
- if (bComment)
- {
- printStringNoNewline(
- inp,
- "The provider of the reaction coordinate, "
- "currently only 'pull' and 'fep-lambda' (free energy lambda state) is supported");
- }
- opt = prefix + "-coord-provider";
- dimParams->eCoordProvider = get_eeenum(inp, opt, eawhcoordprovider_names, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "The coordinate index for this dimension");
- }
- opt = prefix + "-coord-index";
- int coordIndexInput;
- coordIndexInput = get_eint(inp, opt, 1, wi);
-
- /* The pull coordinate indices start at 1 in the input file, at 0 internally */
- dimParams->coordIndex = coordIndexInput - 1;
-
- if (bComment)
- {
- printStringNoNewline(inp, "Start and end values for each coordinate dimension");
- }
- opt = prefix + "-start";
- dimParams->origin = get_ereal(inp, opt, 0., wi);
- opt = prefix + "-end";
- dimParams->end = get_ereal(inp, opt, 0., wi);
-
- if (bComment)
- {
- printStringNoNewline(
- inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
- }
- opt = prefix + "-force-constant";
- dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps, rad^2/ps or ps^-1)");
- }
- opt = prefix + "-diffusion";
- dimParams->diffusion = get_ereal(inp, opt, 0, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp,
- "Diameter that needs to be sampled around a point before it is "
- "considered covered. In FEP dimensions the cover diameter is "
- "specified in lambda states.");
- }
- opt = prefix + "-cover-diameter";
- dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
-}
-
/*! \brief
* Check the parameters of an AWH bias dimension.
*
* \param[in] ir Input parameter struct.
* \param[in,out] wi Struct for bookeeping warnings.
*/
-void checkDimParams(const std::string& prefix, AwhDimParams* dimParams, const t_inputrec* ir, warninp_t wi)
+void checkDimParams(const std::string& prefix, const AwhDimParams& dimParams, const t_inputrec& ir, warninp_t wi)
{
- if (dimParams->eCoordProvider == eawhcoordproviderPULL)
+ 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 == eawhcoordproviderFREE_ENERGY_LAMBDA)
+ else if (dimParams.coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda)
{
- if (ir->efep == efepNO)
+ 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, ir->efep, wi);
+ checkFepLambdaDimParams(prefix, dimParams, ir.fepvals.get(), ir.efep, wi);
}
else
{
}
}
-/*! \brief
- * Check consistency of input at the AWH bias level.
- *
- * \param[in] awhBiasParams AWH bias parameters.
- * \param[in,out] wi Struct for bookkeeping warnings.
- */
-void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
-{
- /* Covering diameter and sharing warning. */
- for (int d = 0; d < awhBiasParams.ndim; d++)
- {
- double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
- if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
- {
- warning(wi,
- "The covering diameter is only relevant to set for bias sharing simulations.");
- }
- }
-}
-
-/*! \brief
- * Read parameters of an AWH bias.
- *
- * \param[in,out] inp Input file entries.
- * \param[in,out] awhBiasParams AWH dimensional parameters.
- * \param[in] prefix Prefix for bias parameters.
- * \param[in,out] wi Struct for bookeeping warnings.
- * \param[in] bComment True if comments should be printed.
- */
-void readBiasParams(std::vector<t_inpfile>* inp,
- AwhBiasParams* awhBiasParams,
- const std::string& prefix,
- warninp_t wi,
- bool bComment)
-{
- if (bComment)
- {
- printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
- }
-
- std::string opt = prefix + "-error-init";
- awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp,
- "Growth rate of the reference histogram determining the bias update "
- "size: exp-linear or linear");
- }
- opt = prefix + "-growth";
- awhBiasParams->eGrowth = get_eeenum(inp, opt, eawhgrowth_names, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp,
- "Start the simulation by equilibrating histogram towards the target "
- "distribution: no or yes");
- }
- opt = prefix + "-equilibrate-histogram";
- awhBiasParams->equilibrateHistogram = (get_eeenum(inp, opt, yesno_names, wi) != 0);
-
- if (bComment)
- {
- printStringNoNewline(
- inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
- }
- opt = prefix + "-target";
- awhBiasParams->eTarget = get_eeenum(inp, opt, eawhtarget_names, 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 = get_eeenum(inp, opt, yesno_names, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
- }
- opt = prefix + "-share-group";
- awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
-
- if (bComment)
- {
- printStringNoNewline(inp, "Dimensionality of the coordinate");
- }
- opt = prefix + "-ndim";
- awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
-
- /* Check this before starting to read the AWH dimension parameters. */
- if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
- {
- gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(),
- awhBiasParams->ndim, c_biasMaxNumDim);
- }
- snew(awhBiasParams->dimParams, awhBiasParams->ndim);
- for (int d = 0; d < awhBiasParams->ndim; d++)
- {
- bComment = bComment && d == 0;
- std::string prefixdim = prefix + formatString("-dim%d", d + 1);
- readDimParams(inp, prefixdim, &awhBiasParams->dimParams[d], wi, bComment);
- }
-}
-
/*! \brief
* Check the parameters of an AWH bias.
*
* \param[in] ir Input parameter struct.
* \param[in,out] wi Struct for bookeeping warnings.
*/
-void checkBiasParams(const AwhBiasParams* awhBiasParams, const std::string& prefix, const t_inputrec* ir, warninp_t wi)
+void checkBiasParams(const AwhBiasParams& awhBiasParams, const std::string& prefix, const t_inputrec& ir, warninp_t wi)
{
std::string opt = prefix + "-error-init";
- if (awhBiasParams->errorInitial <= 0)
+ if (awhBiasParams.initialErrorEstimate() <= 0)
{
gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
}
opt = prefix + "-equilibrate-histogram";
- if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != eawhgrowthEXP_LINEAR)
+ if (awhBiasParams.equilibrateHistogram()
+ && awhBiasParams.growthType() != AwhHistogramGrowthType::ExponentialLinear)
{
auto message =
formatString("Option %s will only have an effect for histogram growth type '%s'.",
- opt.c_str(), EAWHGROWTH(eawhgrowthEXP_LINEAR));
+ opt.c_str(),
+ enumValueToString(AwhHistogramGrowthType::ExponentialLinear));
warning(wi, message);
}
- if ((awhBiasParams->eTarget == eawhtargetLOCALBOLTZMANN)
- && (awhBiasParams->eGrowth == eawhgrowthEXP_LINEAR))
+ if ((awhBiasParams.targetDistribution() == AwhTargetType::LocalBoltzmann)
+ && (awhBiasParams.growthType() == AwhHistogramGrowthType::ExponentialLinear))
{
auto message = formatString(
"Target type '%s' combined with histogram growth type '%s' is not "
"expected to give stable bias updates. You probably want to use growth type "
"'%s' instead.",
- EAWHTARGET(eawhtargetLOCALBOLTZMANN), EAWHGROWTH(eawhgrowthEXP_LINEAR),
- EAWHGROWTH(eawhgrowthLINEAR));
+ enumValueToString(AwhTargetType::LocalBoltzmann),
+ enumValueToString(AwhHistogramGrowthType::ExponentialLinear),
+ enumValueToString(AwhHistogramGrowthType::Linear));
warning(wi, message);
}
opt = prefix + "-target-beta-scaling";
- switch (awhBiasParams->eTarget)
+ switch (awhBiasParams.targetDistribution())
{
- case eawhtargetBOLTZMANN:
- case eawhtargetLOCALBOLTZMANN:
- if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
+ case AwhTargetType::Boltzmann:
+ case AwhTargetType::LocalBoltzmann:
+ if (awhBiasParams.targetBetaScaling() < 0 || awhBiasParams.targetBetaScaling() > 1)
{
- gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
- awhBiasParams->targetBetaScaling, EAWHTARGET(awhBiasParams->eTarget));
+ gmx_fatal(FARGS,
+ "%s = %g is not useful for target type %s.",
+ opt.c_str(),
+ 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,
- EAWHTARGET(awhBiasParams->eTarget));
+ opt.c_str(),
+ awhBiasParams.targetBetaScaling(),
+ enumValueToString(awhBiasParams.targetDistribution()));
}
break;
}
opt = prefix + "-target-cutoff";
- switch (awhBiasParams->eTarget)
+ switch (awhBiasParams.targetDistribution())
{
- case eawhtargetCUTOFF:
- if (awhBiasParams->targetCutoff <= 0)
+ case AwhTargetType::Cutoff:
+ if (awhBiasParams.targetCutoff() <= 0)
{
- gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
- awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
+ gmx_fatal(FARGS,
+ "%s = %g is not useful for target type %s.",
+ opt.c_str(),
+ 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, EAWHTARGET(awhBiasParams->eTarget));
+ opt.c_str(),
+ awhBiasParams.targetCutoff(),
+ enumValueToString(awhBiasParams.targetDistribution()));
}
break;
}
opt = prefix + "-share-group";
- if (awhBiasParams->shareGroup < 0)
+ if (awhBiasParams.shareGroup() < 0)
{
warning_error(wi, "AWH bias share-group should be >= 0");
}
opt = prefix + "-ndim";
- if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
+ if (awhBiasParams.ndim() <= 0 || awhBiasParams.ndim() > c_biasMaxNumDim)
{
- gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(),
- awhBiasParams->ndim, c_biasMaxNumDim);
+ gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), awhBiasParams.ndim(), c_biasMaxNumDim);
}
- if (awhBiasParams->ndim > 2)
+ if (awhBiasParams.ndim() > 2)
{
warning_note(wi,
"For awh-dim > 2 the estimate based on the diffusion and the initial error is "
"currently only a rough guideline."
" You should verify its usefulness for your system before production runs!");
}
- for (int d = 0; d < awhBiasParams->ndim; d++)
+ for (int d = 0; d < awhBiasParams.ndim(); d++)
{
std::string prefixdim = prefix + formatString("-dim%d", d + 1);
- checkDimParams(prefixdim, &awhBiasParams->dimParams[d], ir, wi);
+ checkDimParams(prefixdim, awhBiasParams.dimParams()[d], ir, wi);
}
+}
- /* Check consistencies here that cannot be checked at read time at a lower level. */
- checkInputConsistencyAwhBias(*awhBiasParams, wi);
+/*! \brief
+ * Check consistency of input at the AWH bias level.
+ *
+ * \param[in] awhBiasParams AWH bias parameters.
+ * \param[in,out] wi Struct for bookkeeping warnings.
+ */
+void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
+{
+ /* Covering diameter and sharing warning. */
+ auto awhBiasDimensionParams = awhBiasParams.dimParams();
+ for (const auto& dimensionParam : awhBiasDimensionParams)
+ {
+ double coverDiameter = dimensionParam.coverDiameter();
+ if (awhBiasParams.shareGroup() <= 0 && coverDiameter > 0)
+ {
+ warning(wi,
+ "The covering diameter is only relevant to set for bias sharing simulations.");
+ }
+ }
}
/*! \brief
* Check that we have a shared bias when requesting multisim sharing.
*/
bool haveSharedBias = false;
- for (int k1 = 0; k1 < awhParams.numBias; k1++)
+ auto awhBiasParams = awhParams.awhBiasParams();
+ for (int k1 = 0; k1 < awhParams.numBias(); k1++)
{
- const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1];
+ const AwhBiasParams& awhBiasParams1 = awhBiasParams[k1];
- if (awhBiasParams1.shareGroup > 0)
+ if (awhBiasParams1.shareGroup() > 0)
{
haveSharedBias = true;
}
/* k1 is the reference AWH, k2 is the AWH we compare with (can be equal to k1) */
- for (int k2 = k1; k2 < awhParams.numBias; k2++)
+ for (int k2 = k1; k2 < awhParams.numBias(); k2++)
{
- for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
+ const AwhBiasParams& awhBiasParams2 = awhBiasParams[k2];
+ const auto& dimParams1 = awhBiasParams1.dimParams();
+ const auto& dimParams2 = awhBiasParams2.dimParams();
+ for (int d1 = 0; d1 < gmx::ssize(dimParams1); d1++)
{
- if (awhBiasParams1.dimParams[d1].eCoordProvider == eawhcoordproviderFREE_ENERGY_LAMBDA)
+ 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 == eawhcoordproviderFREE_ENERGY_LAMBDA)
+ if (dimParams2[d2].coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda)
{
continue;
}
/* Give an error if (d1, k1) is different from (d2, k2) but the pull coordinate is the same */
if ((d1 != d2 || k1 != k2)
- && (awhBiasParams1.dimParams[d1].coordIndex == awhBiasParams2.dimParams[d2].coordIndex))
+ && (dimParams1[d1].coordinateIndex() == dimParams2[d2].coordinateIndex()))
{
char errormsg[STRLEN];
sprintf(errormsg,
"dimensions (awh%d-dim%d and awh%d-dim%d). "
"If this is really what you want to do you will have to duplicate "
"this pull coordinate.",
- awhBiasParams1.dimParams[d1].coordIndex + 1, k1 + 1, d1 + 1, k2 + 1,
+ dimParams1[d1].coordinateIndex() + 1,
+ k1 + 1,
+ d1 + 1,
+ k2 + 1,
d2 + 1);
gmx_fatal(FARGS, "%s", errormsg);
}
}
}
- if (awhParams.shareBiasMultisim && !haveSharedBias)
+ if (awhParams.shareBiasMultisim() && !haveSharedBias)
{
warning(wi,
"Sharing of biases over multiple simulations is requested, but no bias is marked "
"this (yet)");
}
}
+
} // namespace
-AwhParams* readAwhParams(std::vector<t_inpfile>* inp, warninp_t wi)
+AwhDimParams::AwhDimParams(std::vector<t_inpfile>* inp, const std::string& prefix, warninp_t wi, bool bComment)
+{
+ std::string opt;
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp,
+ "The provider of the reaction coordinate, "
+ "currently only 'pull' and 'fep-lambda' (free energy lambda state) is supported");
+ }
+
+ opt = prefix + "-coord-provider";
+ eCoordProvider_ = getEnum<AwhCoordinateProviderType>(inp, opt.c_str(), wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "The coordinate index for this dimension");
+ }
+ opt = prefix + "-coord-index";
+ int coordIndexInput;
+ coordIndexInput = get_eint(inp, opt, 1, wi);
+ if (coordIndexInput < 1)
+ {
+ gmx_fatal(FARGS,
+ "Failed to read a valid coordinate index for %s. "
+ "Note that the pull coordinate indexing starts at 1.",
+ opt.c_str());
+ }
+
+ /* The pull coordinate indices start at 1 in the input file, at 0 internally */
+ coordIndex_ = coordIndexInput - 1;
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Start and end values for each coordinate dimension");
+ }
+
+ opt = prefix + "-start";
+ origin_ = get_ereal(inp, opt, 0., wi);
+
+ opt = prefix + "-end";
+ end_ = get_ereal(inp, opt, 0., wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
+ }
+ opt = prefix + "-force-constant";
+ forceConstant_ = get_ereal(inp, opt, 0, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps or ps^-1)");
+ }
+ opt = prefix + "-diffusion";
+ double diffusionValue = get_ereal(inp, opt, 0, wi);
+ if (diffusionValue <= 0)
+ {
+ const double diffusion_default = 1e-5;
+ auto message = formatString(
+ "%s not explicitly set by user. You can choose to use a default "
+ "value (%g nm^2/ps or rad^2/ps) but this may very well be "
+ "non-optimal for your system!",
+ opt.c_str(),
+ diffusion_default);
+ warning(wi, message);
+ diffusionValue = diffusion_default;
+ }
+ diffusion_ = diffusionValue;
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Diameter that needs to be sampled around a point before it is "
+ "considered covered. In FEP dimensions the cover diameter is "
+ "specified in lambda states.");
+ }
+ opt = prefix + "-cover-diameter";
+ coverDiameter_ = get_ereal(inp, opt, 0, wi);
+}
+
+AwhDimParams::AwhDimParams(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(serializer->reading(),
+ "Can not use writing serializer for creating datastructure");
+ serializer->doEnumAsInt(&eCoordProvider_);
+ serializer->doInt(&coordIndex_);
+ serializer->doDouble(&origin_);
+ serializer->doDouble(&end_);
+ serializer->doDouble(&period_);
+ serializer->doDouble(&forceConstant_);
+ serializer->doDouble(&diffusion_);
+ serializer->doDouble(&coordValueInit_);
+ serializer->doDouble(&coverDiameter_);
+}
+
+void AwhDimParams::serialize(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(!serializer->reading(),
+ "Can not use reading serializer for writing datastructure");
+ serializer->doEnumAsInt(&eCoordProvider_);
+ serializer->doInt(&coordIndex_);
+ serializer->doDouble(&origin_);
+ serializer->doDouble(&end_);
+ serializer->doDouble(&period_);
+ serializer->doDouble(&forceConstant_);
+ serializer->doDouble(&diffusion_);
+ serializer->doDouble(&coordValueInit_);
+ serializer->doDouble(&coverDiameter_);
+}
+
+AwhBiasParams::AwhBiasParams(std::vector<t_inpfile>* inp, const std::string& prefix, warninp_t wi, bool bComment)
+{
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
+ }
+
+ std::string opt = prefix + "-error-init";
+ errorInitial_ = get_ereal(inp, opt, 10, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Growth rate of the reference histogram determining the bias update "
+ "size: exp-linear or linear");
+ }
+ opt = prefix + "-growth";
+ eGrowth_ = getEnum<AwhHistogramGrowthType>(inp, opt.c_str(), wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Start the simulation by equilibrating histogram towards the target "
+ "distribution: no or yes");
+ }
+ opt = prefix + "-equilibrate-histogram";
+ equilibrateHistogram_ = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
+
+ if (bComment)
+ {
+ printStringNoNewline(
+ inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
+ }
+ opt = prefix + "-target";
+ eTarget_ = getEnum<AwhTargetType>(inp, opt.c_str(), wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp,
+ "Boltzmann beta scaling factor for target distribution types "
+ "'boltzmann' and 'boltzmann-local'");
+ }
+ opt = prefix + "-target-beta-scaling";
+ targetBetaScaling_ = get_ereal(inp, opt, 0, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
+ }
+ opt = prefix + "-target-cutoff";
+ targetCutoff_ = get_ereal(inp, opt, 0, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
+ }
+ opt = prefix + "-user-data";
+ bUserData_ = getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No;
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
+ }
+ opt = prefix + "-share-group";
+ shareGroup_ = get_eint(inp, opt, 0, wi);
+
+ if (bComment)
+ {
+ printStringNoNewline(inp, "Dimensionality of the coordinate");
+ }
+ opt = prefix + "-ndim";
+ int ndim = get_eint(inp, opt, 0, wi);
+
+ /* Check this before starting to read the AWH dimension parameters. */
+ if (ndim <= 0 || ndim > c_biasMaxNumDim)
+ {
+ gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), ndim, c_biasMaxNumDim);
+ }
+ for (int d = 0; d < ndim; d++)
+ {
+ bComment = bComment && d == 0;
+ std::string prefixdim = prefix + formatString("-dim%d", d + 1);
+ dimParams_.emplace_back(inp, prefixdim, wi, bComment);
+ }
+}
+
+AwhBiasParams::AwhBiasParams(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(serializer->reading(),
+ "Can not use writing serializer to create datastructure");
+ serializer->doEnumAsInt(&eTarget_);
+ serializer->doDouble(&targetBetaScaling_);
+ serializer->doDouble(&targetCutoff_);
+ serializer->doEnumAsInt(&eGrowth_);
+ int temp = 0;
+ serializer->doInt(&temp);
+ bUserData_ = static_cast<bool>(temp);
+ serializer->doDouble(&errorInitial_);
+ int numDimensions = dimParams_.size();
+ serializer->doInt(&numDimensions);
+ serializer->doInt(&shareGroup_);
+ serializer->doBool(&equilibrateHistogram_);
+
+ for (int k = 0; k < numDimensions; k++)
+ {
+ dimParams_.emplace_back(serializer);
+ }
+ /* Check consistencies here that cannot be checked at read time at a lower level. */
+ checkInputConsistencyAwhBias(*this, nullptr);
+}
+
+void AwhBiasParams::serialize(ISerializer* serializer)
+{
+ GMX_RELEASE_ASSERT(!serializer->reading(),
+ "Can not use reading serializer to write datastructure");
+ serializer->doEnumAsInt(&eTarget_);
+ serializer->doDouble(&targetBetaScaling_);
+ serializer->doDouble(&targetCutoff_);
+ serializer->doEnumAsInt(&eGrowth_);
+ int temp = static_cast<int>(bUserData_);
+ serializer->doInt(&temp);
+ serializer->doDouble(&errorInitial_);
+ int numDimensions = ndim();
+ serializer->doInt(&numDimensions);
+ serializer->doInt(&shareGroup_);
+ serializer->doBool(&equilibrateHistogram_);
+
+ for (int k = 0; k < numDimensions; k++)
+ {
+ dimParams_[k].serialize(serializer);
+ }
+}
+
+AwhParams::AwhParams(std::vector<t_inpfile>* inp, 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 = get_eeenum(inp, opt, eawhpotential_names, wi);
+ opt = "awh-potential";
+ potentialEnum_ = getEnum<AwhPotentialType>(inp, opt.c_str(), wi);
printStringNoNewline(inp,
"The random seed used for sampling the umbrella center in the case of "
"umbrella type potential");
- opt = "awh-seed";
- awhParams->seed = get_eint(inp, opt, -1, wi);
- if (awhParams->seed == -1)
+ opt = "awh-seed";
+ seed_ = get_eint(inp, opt, -1, wi);
+ if (seed_ == -1)
{
- awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
- fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
+ seed_ = static_cast<int>(gmx::makeRandomSeed());
+ fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", seed_);
}
printStringNoNewline(inp, "Data output interval in number of steps");
- opt = "awh-nstout";
- awhParams->nstOut = get_eint(inp, opt, 100000, wi);
+ opt = "awh-nstout";
+ nstOut_ = get_eint(inp, opt, 100000, wi);
printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
- opt = "awh-nstsample";
- awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
+ opt = "awh-nstsample";
+ nstSampleCoord_ = get_eint(inp, opt, 10, wi);
printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
- opt = "awh-nsamples-update";
- awhParams->numSamplesUpdateFreeEnergy = get_eint(inp, opt, 10, wi);
+ opt = "awh-nsamples-update";
+ numSamplesUpdateFreeEnergy_ = get_eint(inp, opt, 10, wi);
printStringNoNewline(
inp, "When true, biases with share-group>0 are shared between multiple simulations");
- opt = "awh-share-multisim";
- awhParams->shareBiasMultisim = (get_eeenum(inp, opt, yesno_names, wi) != 0);
+ opt = "awh-share-multisim";
+ shareBiasMultisim_ = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
printStringNoNewline(inp, "The number of independent AWH biases");
- opt = "awh-nbias";
- awhParams->numBias = get_eint(inp, opt, 1, wi);
+ opt = "awh-nbias";
+ int numBias = get_eint(inp, opt, 1, wi);
/* Check this before starting to read the AWH biases. */
- if (awhParams->numBias <= 0)
+ if (numBias <= 0)
{
gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
}
/* Read the parameters specific to each AWH bias */
- snew(awhParams->awhBiasParams, awhParams->numBias);
-
- for (int k = 0; k < awhParams->numBias; k++)
+ for (int k = 0; k < numBias; k++)
{
bool bComment = (k == 0);
std::string prefixawh = formatString("awh%d", k + 1);
- readBiasParams(inp, &awhParams->awhBiasParams[k], prefixawh, wi, bComment);
+ awhBiasParams_.emplace_back(inp, prefixawh, wi, bComment);
}
-
- return awhParams;
+ 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");
- }
-
- bool haveFepLambdaDim = false;
- for (int k = 0; k < awhParams->numBias; k++)
- {
- std::string prefixawh = formatString("awh%d", k + 1);
- checkBiasParams(&awhParams->awhBiasParams[k], prefixawh, ir, wi);
- /* Check if there is a FEP lambda dimension. */
- for (int l = 0; l < awhParams->awhBiasParams[k].ndim; l++)
+ 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++)
{
- if (awhParams->awhBiasParams[k].dimParams[l].eCoordProvider == eawhcoordproviderFREE_ENERGY_LAMBDA)
- {
- haveFepLambdaDim = true;
- break;
- }
+ awhBiasParams_.emplace_back(serializer);
}
}
+ checkInputConsistencyAwh(*this, nullptr);
+}
- if (haveFepLambdaDim)
- {
- if (awhParams->nstSampleCoord % ir->nstcalcenergy != 0)
- {
- opt = "awh-nstsample";
- auto message = formatString(
- "%s (%d) should be a multiple of nstcalcenergy (%d) when using AWH for "
- "sampling an FEP lambda dimension",
- opt.c_str(), awhParams->nstSampleCoord, ir->nstcalcenergy);
- warning_error(wi, message);
- }
- if (awhParams->ePotential != eawhpotentialUMBRELLA)
+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++)
{
- opt = "awh-potential";
- auto message = formatString(
- "%s (%s) must be set to %s when using AWH for sampling an FEP lambda dimension",
- opt.c_str(), eawhpotential_names[awhParams->ePotential],
- eawhpotential_names[eawhpotentialUMBRELLA]);
- warning_error(wi, message);
+ awhBiasParams_[k].serialize(serializer);
}
}
-
- /* 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");
- }
}
/*! \brief
{
double period = 0;
- if (pullCoordParams.eGeom == epullgDIR)
+ if (pullCoordParams.eGeom == PullGroupGeometry::Direction)
{
const real margin = 0.001;
// Make dims periodic when the interval covers > 95%
gmx_fatal(FARGS,
"The AWH interval (%f nm) for a pull coordinate is larger than the "
"box size (%f nm)",
- intervalLength, boxLength);
+ intervalLength,
+ boxLength);
}
if (intervalLength > periodicFraction * boxLength)
}
}
}
- else if (pullCoordParams.eGeom == epullgDIHEDRAL)
+ else if (pullCoordParams.eGeom == PullGroupGeometry::Dihedral)
{
/* The dihedral angle is periodic in -180 to 180 deg */
period = 360;
* \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))
{
gmx_fatal(FARGS,
"For the non-periodic pull coordinates awh%d-dim%d-start (%f) cannot be "
"larger than awh%d-dim%d-end (%f)",
- k + 1, d + 1, origin, k + 1, d + 1, end);
+ k + 1,
+ d + 1,
+ origin,
+ k + 1,
+ d + 1,
+ end);
}
/* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
"values in between "
"minus half a period and plus half a period, i.e. in the interval [%.8g, "
"%.8g].",
- k + 1, d + 1, origin, k + 1, d + 1, end, period, -0.5 * period, 0.5 * period);
+ k + 1,
+ d + 1,
+ origin,
+ k + 1,
+ d + 1,
+ end,
+ period,
+ -0.5 * period,
+ 0.5 * period);
}
/* Warn if the pull initial coordinate value is not in the grid */
"(%.8g). "
"This can lead to large initial forces pulling the coordinate towards the "
"sampling interval.",
- coordValueInit, coordIndex + 1, k + 1, d + 1, origin, k + 1, d + 1, end);
+ coordValueInit,
+ coordIndex + 1,
+ k + 1,
+ d + 1,
+ origin,
+ k + 1,
+ d + 1,
+ end);
warning(wi, message);
}
}
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 == epullgDIRPBC)
+ if (pullCoordParams.eGeom == PullGroupGeometry::DirectionPBC)
{
gmx_fatal(FARGS,
"AWH does not support pull geometry '%s'. "
"If the maximum distance between the groups is always "
"less than half the box size, "
"you can use geometry '%s' instead.",
- EPULLGEOM(epullgDIRPBC), EPULLGEOM(epullgDIR));
+ enumValueToString(PullGroupGeometry::DirectionPBC),
+ enumValueToString(PullGroupGeometry::Direction));
}
- dimParams->period = get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin);
+ dimParams->setPeriod(
+ get_pull_coord_period(pullCoordParams, pbc, dimParams->end() - dimParams->origin()));
// We would like to check for scaling, but we don't have the full inputrec available here
- if (dimParams->period > 0
- && !(pullCoordParams.eGeom == epullgANGLE || pullCoordParams.eGeom == epullgDIHEDRAL))
+ if (dimParams->period() > 0
+ && !(pullCoordParams.eGeom == PullGroupGeometry::Angle
+ || pullCoordParams.eGeom == PullGroupGeometry::Dihedral))
{
bool coordIsScaled = false;
for (int d2 = 0; d2 < DIM; d2++)
{
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, EPULLGEOM(pullCoordParams.eGeom));
+ dimIndex + 1,
+ biasIndex + 1,
+ enumValueToString(pullCoordParams.eGeom));
warning(wi, mesg.c_str());
}
}
/* The initial coordinate value, converted to external user units. */
- dimParams->coordValueInit = get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
-
- dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(&pullCoordParams);
+ double initialCoordinate = get_pull_coord_value(pull_work, dimParams->coordinateIndex(), pbc);
+ initialCoordinate *= pull_conversion_factor_internal2userinput(pullCoordParams);
+ dimParams->setInitialCoordinate(initialCoordinate);
}
void setStateDependentAwhParams(AwhParams* awhParams,
t_pbc pbc;
set_pbc(&pbc, pbcType, box);
- for (int k = 0; k < awhParams->numBias; k++)
+ auto awhBiasParams = awhParams->awhBiasParams();
+ for (int k = 0; k < awhParams->numBias(); k++)
{
- AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
- for (int d = 0; d < awhBiasParams->ndim; d++)
+ auto awhBiasDimensionParams = awhBiasParams[k].dimParams();
+ for (int d = 0; d < awhBiasParams[k].ndim(); d++)
{
- AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
- if (dimParams->eCoordProvider == eawhcoordproviderPULL)
+ AwhDimParams* dimParams = &awhBiasDimensionParams[d];
+ if (dimParams->coordinateProvider() == AwhCoordinateProviderType::Pull)
{
- setStateDependentAwhPullDimParams(dimParams, k, d, &pull_params, pull_work, pbc,
- compressibility, wi);
+ setStateDependentAwhPullDimParams(
+ 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");
+ }
+
+ bool haveFepLambdaDim = false;
+ const auto awhBiasParams = awhParams.awhBiasParams();
+ for (int k = 0; k < awhParams.numBias() && !haveFepLambdaDim; k++)
+ {
+ std::string prefixawh = formatString("awh%d", k + 1);
+ checkBiasParams(awhBiasParams[k], prefixawh, ir, wi);
+ /* Check if there is a FEP lambda dimension. */
+ const auto dimParams = awhBiasParams[k].dimParams();
+ haveFepLambdaDim = std::any_of(dimParams.begin(), dimParams.end(), [](const auto& dimParam) {
+ return dimParam.coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda;
+ });
+ }
+
+ if (haveFepLambdaDim)
+ {
+ if (awhParams.nstSampleCoord() % ir.nstcalcenergy != 0)
+ {
+ opt = "awh-nstsample";
+ auto message = formatString(
+ "%s (%d) should be a multiple of nstcalcenergy (%d) when using AWH for "
+ "sampling an FEP lambda dimension",
+ opt.c_str(),
+ awhParams.nstSampleCoord(),
+ ir.nstcalcenergy);
+ warning_error(wi, message);
+ }
+ if (awhParams.potential() != AwhPotentialType::Umbrella)
+ {
+ opt = "awh-potential";
+ auto message = formatString(
+ "%s (%s) must be set to %s when using AWH for sampling an FEP lambda dimension",
+ opt.c_str(),
+ enumValueToString(awhParams.potential()),
+ enumValueToString(AwhPotentialType::Umbrella));
+ warning_error(wi, message);
+ }
+ }
+
+ if (ir.init_step != 0)
+ {
+ warning_error(wi, "With AWH init-step should be 0");
+ }
+}
+
+bool awhHasFepLambdaDimension(const AwhParams& awhParams)
+{
+ for (const auto& biasParams : awhParams.awhBiasParams())
+ {
+ for (const auto& dimParams : biasParams.dimParams())
+ {
+ if (dimParams.coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda)
+ {
+ return true;
+ }
+ }
+ }
+
+ return false;
+}
+
} // namespace gmx
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff