From: Mark Abraham Date: Fri, 29 Nov 2019 12:52:13 +0000 (+0100) Subject: Merge branch release-2019 into release-2020 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=475b9a11a7bcc07d3cb20cd5d49e76ec538c9a2a;p=alexxy%2Fgromacs.git Merge branch release-2019 into release-2020 Change-Id: I138d910b0d5cb8555460f87345c4e312e528a3ce --- 475b9a11a7bcc07d3cb20cd5d49e76ec538c9a2a diff --cc src/gromacs/awh/biasparams.h index 938eefce1e,2cf134d224..812963c27a --- a/src/gromacs/awh/biasparams.h +++ b/src/gromacs/awh/biasparams.h @@@ -69,160 -69,175 +69,163 @@@ class GridAxis */ class BiasParams { - public: - /*! \brief Switch to turn off update skips, useful for testing. - */ - enum class DisableUpdateSkips - { - no, /**< Allow update skips (when supported by the method) */ - yes /**< Disable update skips */ - }; - - /*! \brief - * Check if the parameters permit skipping updates. - * - * Generally, we can skip updates of points that are non-local - * at the time of the update if we for later times, when the points - * with skipped updates have become local, know exactly how to apply - * the previous updates. The free energy updates only depend - * on local sampling, but the histogram rescaling factors - * generally depend on the histogram size (all samples). - * If the histogram size is kept constant or the scaling factors - * are trivial, this is not a problem. However, if the histogram growth - * is scaled down by some factor the size at the time of the update - * needs to be known. It would be fairly simple to, for a deterministically - * growing histogram, backtrack and calculate this value, but currently - * we just disallow this case. This is not a restriction because it - * only affects the local Boltzmann target type for which every update - * is currently anyway global because the target is always updated globally. - * - * \returns true when we can skip updates. - */ - inline bool skipUpdates() const - { - return (!disableUpdateSkips_ && localWeightScaling == 1); - } - - /*! \brief - * Returns the the radius that needs to be sampled around a point before it is considered covered. - */ - inline const awh_ivec &coverRadius() const - { - return coverRadius_; - } - - /*! \brief - * Returns whether we should sample the coordinate. - * - * \param[in] step The MD step number. - */ - inline bool isSampleCoordStep(int64_t step) const - { - return (step > 0 && step % numStepsSampleCoord_ == 0); - } - - /*! \brief - * Returns whether we should update the free energy. - * - * \param[in] step The MD step number. - */ - inline bool isUpdateFreeEnergyStep(int64_t step) const - { - int stepIntervalUpdateFreeEnergy = numSamplesUpdateFreeEnergy_*numStepsSampleCoord_; - return (step > 0 && step % stepIntervalUpdateFreeEnergy == 0); - } - - /*! \brief - * Returns whether we should update the target distribution. - * - * \param[in] step The MD step number. - */ - inline bool isUpdateTargetStep(int64_t step) const - { - return step % numStepsUpdateTarget_ == 0; - } - - /*! \brief - * Returns if to do checks for covering in the initial stage. - * - * To avoid overhead due to expensive checks, we do not check - * at every free energy update. However, if checks are - * performed too rarely the detection of coverings will be - * delayed, ultimately affecting free energy convergence. - * - * \param[in] step Time step. - * \returns true at steps where checks should be performed. - * \note Only returns true at free energy update steps. - */ - bool isCheckCoveringStep(int64_t step) const - { - return step > 0 && (step % numStepsCheckCovering_ == 0); - } - - /*! \brief - * Returns if to perform checks for anomalies in the histogram. - * - * To avoid overhead due to expensive checks, we do not check - * at every free energy update. These checks are only used for - * warning the user and can be made as infrequently as - * neccessary without affecting the algorithm itself. - * - * \param[in] step Time step. - * \returns true at steps where checks should be performed. - * \note Only returns true at free energy update steps. - * \todo Currently this function just calls isCheckCoveringStep but the checks could be done less frequently. - */ - bool isCheckHistogramForAnomaliesStep(int64_t step) const - { - return isCheckCoveringStep(step); - } - - /*! \brief Constructor. - * - * The local Boltzmann target distibution is defined by - * 1) Adding the sampled weights instead of the target weights to the reference weight histogram. - * 2) Scaling the weights of these samples by the beta scaling factor. - * 3) Setting the target distribution equal the reference weight histogram. - * This requires the following special update settings: - * localWeightScaling = targetParam - * idealWeighthistUpdate = false - * Note: these variables could in principle be set to something else also for other target distribution types. - * However, localWeightScaling < 1 is in general expected to give lower efficiency and, except for local Boltzmann, - * idealWeightHistUpdate = false gives (in my experience) unstable, non-converging results. - * - * \param[in] awhParams AWH parameters. - * \param[in] awhBiasParams Bias parameters. - * \param[in] dimParams Bias dimension parameters. - * \param[in] beta 1/(k_B T) in units of 1/(kJ/mol), should be > 0. - * \param[in] mdTimeStep The MD time step. - * \param[in] numSharingSimulations The number of simulations to share the bias across. - * \param[in] gridAxis The grid axes. - * \param[in] disableUpdateSkips If to disable update skips, useful for testing. - * \param[in] biasIndex Index of the bias. - */ - BiasParams(const AwhParams &awhParams, - const AwhBiasParams &awhBiasParams, - const std::vector &dimParams, - double beta, - double mdTimeStep, - DisableUpdateSkips disableUpdateSkips, - int numSharingSimulations, - const std::vector &gridAxis, - int biasIndex); - - /* Data members */ - const double invBeta; /**< 1/beta = kT in kJ/mol */ - private: - const int64_t numStepsSampleCoord_; /**< Number of steps per coordinate value sample. */ - public: - const int numSamplesUpdateFreeEnergy_; /**< Number of samples per free energy update. */ - private: - const int64_t numStepsUpdateTarget_; /**< Number of steps per updating the target distribution. */ - const int64_t numStepsCheckCovering_; /**< Number of steps per checking for covering. */ - public: - const int eTarget; /**< Type of target distribution. */ - const double freeEnergyCutoffInKT; /**< Free energy cut-off in kT for cut-off target distribution. */ - const double temperatureScaleFactor; /**< Temperature scaling factor for temperature scaled targed distributions. */ - const bool idealWeighthistUpdate; /**< Update reference weighthistogram using the target distribution? Otherwise use the realized distribution. */ - const int numSharedUpdate; /**< The number of (multi-)simulations sharing the bias update */ - const double updateWeight; /**< The probability weight accumulated for each update. */ - const double localWeightScaling; /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */ - const double initialErrorInKT; /**< Estimated initial free energy error in kT. */ - const double initialHistogramSize; /**< Initial reference weight histogram size. */ - private: - awh_ivec coverRadius_; /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */ - public: - const bool convolveForce; /**< True if we convolve the force, false means use MC between umbrellas. */ - const int biasIndex; /**< Index of the bias, used as a second random seed and for priting. */ - private: - const bool disableUpdateSkips_; /**< If true, we disallow update skips, even when the method supports it. */ +public: + /*! \brief Switch to turn off update skips, useful for testing. + */ + enum class DisableUpdateSkips + { + no, /**< Allow update skips (when supported by the method) */ + yes /**< Disable update skips */ + }; + + /*! \brief + * Check if the parameters permit skipping updates. + * + * Generally, we can skip updates of points that are non-local + * at the time of the update if we for later times, when the points + * with skipped updates have become local, know exactly how to apply + * the previous updates. The free energy updates only depend + * on local sampling, but the histogram rescaling factors + * generally depend on the histogram size (all samples). + * If the histogram size is kept constant or the scaling factors + * are trivial, this is not a problem. However, if the histogram growth + * is scaled down by some factor the size at the time of the update + * needs to be known. It would be fairly simple to, for a deterministically + * growing histogram, backtrack and calculate this value, but currently + * we just disallow this case. This is not a restriction because it + * only affects the local Boltzmann target type for which every update + * is currently anyway global because the target is always updated globally. + * + * \returns true when we can skip updates. + */ + inline bool skipUpdates() const { return (!disableUpdateSkips_ && localWeightScaling == 1); } + + /*! \brief + * Returns the radius that needs to be sampled around a point before it is considered covered. + */ + inline const awh_ivec& coverRadius() const { return coverRadius_; } + + /*! \brief + * Returns whether we should sample the coordinate. + * + * \param[in] step The MD step number. + */ + inline bool isSampleCoordStep(int64_t step) const + { + return (step > 0 && step % numStepsSampleCoord_ == 0); + } + + /*! \brief + * Returns whether we should update the free energy. + * + * \param[in] step The MD step number. + */ + inline bool isUpdateFreeEnergyStep(int64_t step) const + { + int stepIntervalUpdateFreeEnergy = numSamplesUpdateFreeEnergy_ * numStepsSampleCoord_; + return (step > 0 && step % stepIntervalUpdateFreeEnergy == 0); + } + + /*! \brief + * Returns whether we should update the target distribution. + * + * \param[in] step The MD step number. + */ + inline bool isUpdateTargetStep(int64_t step) const { return step % numStepsUpdateTarget_ == 0; } + + /*! \brief + * Returns if to do checks for covering in the initial stage. + * + * To avoid overhead due to expensive checks, we do not check + * at every free energy update. However, if checks are + * performed too rarely the detection of coverings will be + * delayed, ultimately affecting free energy convergence. + * + * \param[in] step Time step. + * \returns true at steps where checks should be performed. + * \note Only returns true at free energy update steps. + */ - bool isCheckCoveringStep(int64_t step) const { return step % numStepsCheckCovering_ == 0; } ++ bool isCheckCoveringStep(int64_t step) const ++ { ++ return step > 0 && (step % numStepsCheckCovering_ == 0); ++ } + + /*! \brief + * Returns if to perform checks for anomalies in the histogram. + * + * To avoid overhead due to expensive checks, we do not check + * at every free energy update. These checks are only used for + * warning the user and can be made as infrequently as + * neccessary without affecting the algorithm itself. + * + * \param[in] step Time step. + * \returns true at steps where checks should be performed. + * \note Only returns true at free energy update steps. + * \todo Currently this function just calls isCheckCoveringStep but the checks could be done less frequently. + */ + bool isCheckHistogramForAnomaliesStep(int64_t step) const { return isCheckCoveringStep(step); } + + /*! \brief Constructor. + * + * The local Boltzmann target distibution is defined by + * 1) Adding the sampled weights instead of the target weights to the reference weight histogram. + * 2) Scaling the weights of these samples by the beta scaling factor. + * 3) Setting the target distribution equal the reference weight histogram. + * This requires the following special update settings: + * localWeightScaling = targetParam + * idealWeighthistUpdate = false + * Note: these variables could in principle be set to something else also for other target distribution types. + * However, localWeightScaling < 1 is in general expected to give lower efficiency and, except for local Boltzmann, + * idealWeightHistUpdate = false gives (in my experience) unstable, non-converging results. + * + * \param[in] awhParams AWH parameters. + * \param[in] awhBiasParams Bias parameters. + * \param[in] dimParams Bias dimension parameters. + * \param[in] beta 1/(k_B T) in units of 1/(kJ/mol), should be > 0. + * \param[in] mdTimeStep The MD time step. + * \param[in] numSharingSimulations The number of simulations to share the bias across. + * \param[in] gridAxis The grid axes. + * \param[in] disableUpdateSkips If to disable update skips, useful for testing. + * \param[in] biasIndex Index of the bias. + */ + BiasParams(const AwhParams& awhParams, + const AwhBiasParams& awhBiasParams, + const std::vector& dimParams, + double beta, + double mdTimeStep, + DisableUpdateSkips disableUpdateSkips, + int numSharingSimulations, + const std::vector& gridAxis, + int biasIndex); + + /* Data members */ + const double invBeta; /**< 1/beta = kT in kJ/mol */ +private: + const int64_t numStepsSampleCoord_; /**< Number of steps per coordinate value sample. */ +public: + const int numSamplesUpdateFreeEnergy_; /**< Number of samples per free energy update. */ +private: + const int64_t numStepsUpdateTarget_; /**< Number of steps per updating the target distribution. */ + const int64_t numStepsCheckCovering_; /**< Number of steps per checking for covering. */ +public: + const int eTarget; /**< Type of target distribution. */ + const double freeEnergyCutoffInKT; /**< Free energy cut-off in kT for cut-off target distribution. */ + const double temperatureScaleFactor; /**< Temperature scaling factor for temperature scaled targed distributions. */ + const bool idealWeighthistUpdate; /**< Update reference weighthistogram using the target distribution? Otherwise use the realized distribution. */ + const int numSharedUpdate; /**< The number of (multi-)simulations sharing the bias update */ + const double updateWeight; /**< The probability weight accumulated for each update. */ + const double localWeightScaling; /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */ + const double initialErrorInKT; /**< Estimated initial free energy error in kT. */ + const double initialHistogramSize; /**< Initial reference weight histogram size. */ +private: + awh_ivec coverRadius_; /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */ +public: + const bool convolveForce; /**< True if we convolve the force, false means use MC between umbrellas. */ + const int biasIndex; /**< Index of the bias, used as a second random seed and for priting. */ +private: + const bool disableUpdateSkips_; /**< If true, we disallow update skips, even when the method supports it. */ }; -} // namespace gmx +} // namespace gmx #endif /* GMX_AWH_BIASPARAMS_H */ diff --cc src/gromacs/awh/read_params.cpp index 55837e1f88,0000000000..2b27717731 mode 100644,000000..100644 --- a/src/gromacs/awh/read_params.cpp +++ b/src/gromacs/awh/read_params.cpp @@@ -1,837 -1,0 +1,843 @@@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 1991-2000, University of Groningen, The Netherlands. + * Copyright (c) 2001-2004, The GROMACS development team. + * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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 "read_params.h" + +#include "gromacs/awh/awh.h" +#include "gromacs/fileio/readinp.h" +#include "gromacs/fileio/warninp.h" +#include "gromacs/math/units.h" +#include "gromacs/math/utilities.h" +#include "gromacs/math/vec.h" +#include "gromacs/mdtypes/awh_params.h" +#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/mdtypes/md_enums.h" +#include "gromacs/mdtypes/pull_params.h" +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/pulling/pull.h" +#include "gromacs/random/seed.h" +#include "gromacs/utility/cstringutil.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/smalloc.h" +#include "gromacs/utility/stringutil.h" + +#include "biasparams.h" +#include "biassharing.h" + +namespace gmx +{ + +const char* eawhtarget_names[eawhtargetNR + 1] = { "constant", "cutoff", "boltzmann", + "local-boltzmann", nullptr }; + +const char* eawhgrowth_names[eawhgrowthNR + 1] = { "exp-linear", "linear", nullptr }; + +const char* eawhpotential_names[eawhpotentialNR + 1] = { "convolved", "umbrella", nullptr }; + +const char* eawhcoordprovider_names[eawhcoordproviderNR + 1] = { "pull", nullptr }; + +/*! \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] pull_params Pull parameters. + * \param[in,out] wi Struct for bookeeping warnings. + * \param[in] bComment True if comments should be printed. + */ +static void readDimParams(std::vector* inp, + const std::string& prefix, + AwhDimParams* dimParams, + const pull_params_t* pull_params, + warninp_t wi, + bool bComment) +{ + std::string opt; + if (bComment) + { + printStringNoNewline( + inp, "The provider of the reaction coordinate, currently only pull 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); + 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 */ + dimParams->coordIndex = coordIndexInput - 1; + + /* The pull settings need to be consistent with the AWH settings */ + if (!(pull_params->coord[dimParams->coordIndex].eType == epullEXTERNAL)) + { + gmx_fatal(FARGS, "AWH biasing can only be applied to pull type %s", EPULLTYPE(epullEXTERNAL)); + } + + if (dimParams->coordIndex >= pull_params->ncoord) + { + gmx_fatal(FARGS, + "The given AWH coordinate index (%d) is larger than the number of pull " + "coordinates (%d)", + coordIndexInput, pull_params->ncoord); + } + if (pull_params->coord[dimParams->coordIndex].rate != 0) + { + auto message = formatString( + "Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate", + coordIndexInput, pull_params->coord[dimParams->coordIndex].rate); + warning_error(wi, message); + } + + /* Grid params for each axis */ + int eGeom = pull_params->coord[dimParams->coordIndex].eGeom; + + 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 (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); + warning(wi, message); + } + /* Check that the requested interval is in allowed range */ + if (eGeom == epullgDIST) + { + 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)); + } + } + else if (eGeom == epullgANGLE || eGeom == epullgANGLEAXIS) + { + 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)); + } + } + else if (eGeom == epullgDIHEDRAL) + { + 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)); + } + } + + 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 (dimParams->forceConstant <= 0) + { + warning_error(wi, "The force AWH bias force constant should be > 0"); + } + + if (bComment) + { + printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps)"); + } + opt = prefix + "-diffusion"; + dimParams->diffusion = get_ereal(inp, opt, 0, wi); + + if (dimParams->diffusion <= 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); + dimParams->diffusion = diffusion_default; + } + + if (bComment) + { + printStringNoNewline(inp, + "Diameter that needs to be sampled around a point before it is " + "considered covered."); + } + opt = prefix + "-cover-diameter"; + dimParams->coverDiameter = get_ereal(inp, opt, 0, wi); + + if (dimParams->coverDiameter < 0) + { + gmx_fatal(FARGS, "%s (%g) cannot be negative.", opt.c_str(), dimParams->coverDiameter); + } +} + +/*! \brief + * Check consistency of input at the AWH bias level. + * + * \param[in] awhBiasParams AWH bias parameters. + * \param[in,out] wi Struct for bookkeeping warnings. + */ +static 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] ir Input parameter struct. + * \param[in,out] wi Struct for bookeeping warnings. + * \param[in] bComment True if comments should be printed. + */ +static void read_bias_params(std::vector* inp, + AwhBiasParams* awhBiasParams, + 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"; + /* We allow using a default value here without warning (but warn the user if the diffusion constant is not set). */ + awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi); + if (awhBiasParams->errorInitial <= 0) + { + gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str()); + } + + 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 (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != eawhgrowthEXP_LINEAR) + { + auto message = + formatString("Option %s will only have an effect for histogram growth type '%s'.", + opt.c_str(), EAWHGROWTH(eawhgrowthEXP_LINEAR)); + warning(wi, message); + } + + 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 ((awhBiasParams->eTarget == eawhtargetLOCALBOLTZMANN) + && (awhBiasParams->eGrowth == eawhgrowthEXP_LINEAR)) + { + 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)); + warning(wi, message); + } + + 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); + + switch (awhBiasParams->eTarget) + { + case eawhtargetBOLTZMANN: + case eawhtargetLOCALBOLTZMANN: + 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)); + } + break; + default: + 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)); + } + break; + } + + 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); + + switch (awhBiasParams->eTarget) + { + case eawhtargetCUTOFF: + if (awhBiasParams->targetCutoff <= 0) + { + gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(), + awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget)); + } + break; + default: + 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)); + } + break; + } + + 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 (awhBiasParams->shareGroup < 0) + { + warning_error(wi, "AWH bias share-group should be >= 0"); + } + + if (bComment) + { + printStringNoNewline(inp, "Dimensionality of the coordinate"); + } + opt = prefix + "-ndim"; + awhBiasParams->ndim = get_eint(inp, opt, 0, wi); + + 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); + } + 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!"); + } + 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], ir->pull, wi, bComment); + } + + /* 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 level. + * + * \param[in] awhParams AWH parameters. + * \param[in,out] wi Struct for bookkeeping warnings. + */ +static void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi) +{ + /* Each pull coord can map to at most 1 AWH coord. + * Check that we have a shared bias when requesting multisim sharing. + */ + bool haveSharedBias = false; + for (int k1 = 0; k1 < awhParams.numBias; k1++) + { + const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1]; + + 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 d1 = 0; d1 < awhBiasParams1.ndim; d1++) + { + 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++) + { + /* 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)) + { + char errormsg[STRLEN]; + sprintf(errormsg, + "One pull coordinate (%d) cannot be mapped to two separate AWH " + "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, + d2 + 1); + gmx_fatal(FARGS, "%s", errormsg); + } + } + } + } + } + + if (awhParams.shareBiasMultisim && !haveSharedBias) + { + warning(wi, + "Sharing of biases over multiple simulations is requested, but no bias is marked " + "as shared (share-group > 0)"); + } + + /* mdrun does not support this (yet), but will check again */ + if (haveBiasSharingWithinSimulation(awhParams)) + { + warning(wi, + "You have shared biases within a single simulation, but mdrun does not support " + "this (yet)"); + } +} + +AwhParams* readAndCheckAwhParams(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 = get_eeenum(inp, opt, eawhpotential_names, 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) + { + awhParams->seed = static_cast(gmx::makeRandomSeed()); + fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed); + } + + printStringNoNewline(inp, "Data output interval in number of steps"); + opt = "awh-nstout"; + awhParams->nstOut = get_eint(inp, opt, 100000, wi); + 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); + } + + printStringNoNewline(inp, "Coordinate sampling interval in number of steps"); + opt = "awh-nstsample"; + awhParams->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); + if (awhParams->numSamplesUpdateFreeEnergy <= 0) + { + warning_error(wi, opt + " needs to be an integer > 0"); + } + + 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); + + printStringNoNewline(inp, "The number of independent AWH biases"); + opt = "awh-nbias"; + awhParams->numBias = get_eint(inp, opt, 1, wi); + if (awhParams->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++) + { + bool bComment = (k == 0); + std::string prefixawh = formatString("awh%d", k + 1); + read_bias_params(inp, &awhParams->awhBiasParams[k], prefixawh, ir, wi, bComment); + } + + /* 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"); + } + + return awhParams; +} + +/*! \brief + * Gets the period of a pull coordinate. + * + * \param[in] pullCoordParams The parameters for the pull coordinate. + * \param[in] pbc The PBC setup + * \param[in] intervalLength The length of the AWH interval for this pull coordinate + * \returns the period (or 0 if not periodic). + */ +static double get_pull_coord_period(const t_pull_coord& pullCoordParams, const t_pbc& pbc, const real intervalLength) +{ + double period = 0; + + if (pullCoordParams.eGeom == epullgDIR) + { + const real margin = 0.001; + // Make dims periodic when the interval covers > 95% + const real periodicFraction = 0.95; + + // Check if the pull direction is along a box vector + for (int dim = 0; dim < pbc.ndim_ePBC; dim++) + { + const real boxLength = norm(pbc.box[dim]); + const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]); + if (innerProduct >= (1 - margin) * boxLength && innerProduct <= (1 + margin) * boxLength) + { - GMX_RELEASE_ASSERT(intervalLength < (1 + margin) * boxLength, - "We have checked before that interval <= period"); ++ if (intervalLength > (1 + margin) * boxLength) ++ { ++ gmx_fatal(FARGS, ++ "The AWH interval (%f nm) for a pull coordinate is larger than the " ++ "box size (%f nm)", ++ intervalLength, boxLength); ++ } ++ + if (intervalLength > periodicFraction * boxLength) + { + period = boxLength; + } + } + } + } + else if (pullCoordParams.eGeom == epullgDIHEDRAL) + { + /* The dihedral angle is periodic in -180 to 180 deg */ + period = 360; + } + + return period; +} + +/*! \brief + * Checks if the given interval is defined in the correct periodic interval. + * + * \param[in] origin Start value of interval. + * \param[in] end End value of interval. + * \param[in] period Period (or 0 if not periodic). + * \returns true if the end point values are in the correct periodic interval. + */ +static bool intervalIsInPeriodicInterval(double origin, double end, double period) +{ + return (period == 0) || (std::fabs(origin) <= 0.5 * period && std::fabs(end) <= 0.5 * period); +} + +/*! \brief + * Checks if a value is within an interval. + * + * \param[in] origin Start value of interval. + * \param[in] end End value of interval. + * \param[in] period Period (or 0 if not periodic). + * \param[in] value Value to check. + * \returns true if the value is within the interval. + */ +static bool valueIsInInterval(double origin, double end, double period, double value) +{ + bool bIn_interval; + + if (period > 0) + { + if (origin < end) + { + /* The interval closes within the periodic interval */ + bIn_interval = (value >= origin) && (value <= end); + } + else + { + /* The interval wraps around the periodic boundary */ + bIn_interval = ((value >= origin) && (value <= 0.5 * period)) + || ((value >= -0.5 * period) && (value <= end)); + } + } + else + { + bIn_interval = (value >= origin) && (value <= end); + } + + return bIn_interval; +} + +/*! \brief + * Check if the starting configuration is consistent with the given interval. + * + * \param[in] awhParams AWH parameters. + * \param[in,out] wi Struct for bookeeping warnings. + */ +static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t wi) +{ + for (int k = 0; k < awhParams->numBias; k++) + { + AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k]; + for (int d = 0; d < awhBiasParams->ndim; d++) + { + AwhDimParams* dimParams = &awhBiasParams->dimParams[d]; + int coordIndex = dimParams->coordIndex; + double origin = dimParams->origin, end = dimParams->end, period = dimParams->period; + double coordValueInit = dimParams->coordValueInit; + + 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); + } + + /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval. + Make sure the AWH interval is within this reference interval. + + Note: we could fairly simply allow using a more general interval (e.g. [x, x + period]) but it complicates + things slightly and I don't see that there is a great need for it. It would also mean that the interval would + depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180, + independent of AWH parameters. + */ + if (!intervalIsInPeriodicInterval(origin, end, period)) + { + gmx_fatal(FARGS, + "When using AWH with periodic pull coordinate geometries " + "awh%d-dim%d-start (%.8g) and " + "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take " + "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); + } + + /* Warn if the pull initial coordinate value is not in the grid */ + if (!valueIsInInterval(origin, end, period, coordValueInit)) + { + auto message = formatString( + "The initial coordinate value (%.8g) for pull coordinate index %d falls " + "outside " + "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end " + "(%.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); + warning(wi, message); + } + } + } +} + +void setStateDependentAwhParams(AwhParams* awhParams, + const pull_params_t* pull_params, + pull_t* pull_work, + const matrix box, + int ePBC, + const tensor& compressibility, + const t_grpopts* inputrecGroupOptions, + warninp_t wi) +{ + /* The temperature is not really state depenendent but is not known + * when read_awhParams is called (in get ir). + * It is known first after do_index has been called in grompp.cpp. + */ + if (inputrecGroupOptions->ref_t == nullptr || inputrecGroupOptions->ref_t[0] <= 0) + { + gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0"); + } + for (int i = 1; i < inputrecGroupOptions->ngtc; i++) + { + if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0]) + { + gmx_fatal(FARGS, + "AWH biasing is currently only supported for identical temperatures for all " + "temperature coupling groups"); + } + } + + t_pbc pbc; + set_pbc(&pbc, ePBC, box); + + for (int k = 0; k < awhParams->numBias; k++) + { + AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k]; + for (int d = 0; d < awhBiasParams->ndim; d++) + { + AwhDimParams* dimParams = &awhBiasParams->dimParams[d]; + const t_pull_coord& pullCoordParams = pull_params->coord[dimParams->coordIndex]; + + if (pullCoordParams.eGeom == epullgDIRPBC) + { + 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)); + } + + dimParams->period = + 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)) + { + bool coordIsScaled = false; + for (int d2 = 0; d2 < DIM; d2++) + { + if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0) + { + coordIsScaled = true; + } + } + if (coordIsScaled) + { + 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 corresponding " + "box vector, this is not supported.", + d + 1, k + 1, EPULLGEOM(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); + } + } + checkInputConsistencyInterval(awhParams, wi); + + /* Register AWH as external potential with pull to check consistency. */ + Awh::registerAwhWithPull(*awhParams, pull_work); +} + +} // namespace gmx diff --cc src/gromacs/fileio/pdbio.cpp index fa4be6e34b,e475e66a2c..3466d8d4f1 --- a/src/gromacs/fileio/pdbio.cpp +++ b/src/gromacs/fileio/pdbio.cpp @@@ -765,11 -824,10 +765,9 @@@ static void gmx_conect_addline(gmx_cone n = sscanf(line, format.c_str(), &aj); if (n == 1) { - srenew(con->conect, ++con->nconect); - con->conect[con->nconect - 1].ai = ai - 1; - con->conect[con->nconect - 1].aj = aj - 1; + gmx_conect_add(con, ai - 1, aj - 1); /* to prevent duplicated records */ } - } - while (n == 1); + } while (n == 1); } } diff --cc src/gromacs/pulling/pull.cpp index 4e159aa760,aaa86cc981..c20deec924 --- a/src/gromacs/pulling/pull.cpp +++ b/src/gromacs/pulling/pull.cpp @@@ -510,13 -512,19 +510,17 @@@ static void low_get_pull_coord_dr(cons /* This function returns the distance based on the contents of the pull struct. * pull->coord[coord_ind].dr, and potentially vector, are updated. */ -static void get_pull_coord_dr(struct pull_t *pull, - int coord_ind, - const t_pbc *pbc) +static void get_pull_coord_dr(struct pull_t* pull, int coord_ind, const t_pbc* pbc) { - pull_coord_work_t *pcrd = &pull->coord[coord_ind]; - PullCoordSpatialData &spatialData = pcrd->spatialData; + pull_coord_work_t* pcrd = &pull->coord[coord_ind]; + PullCoordSpatialData& spatialData = pcrd->spatialData; - double md2; + double md2; - if (pcrd->params.eGeom == epullgDIRPBC) + /* With AWH pulling we allow for periodic pulling with geometry=direction. + * TODO: Store a periodicity flag instead of checking for external pull provider. + */ - if (pcrd->params.eGeom == epullgDIRPBC || (pcrd->params.eGeom == epullgDIR && - pcrd->params.eType == epullEXTERNAL)) ++ if (pcrd->params.eGeom == epullgDIRPBC ++ || (pcrd->params.eGeom == epullgDIR && pcrd->params.eType == epullEXTERNAL)) { md2 = -1; }