src/gromacs/applied_forces/awh/read_params.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "read_params.h"
  41
  42 #include "gromacs/applied_forces/awh/awh.h"
  43 #include "gromacs/fileio/readinp.h"
  44 #include "gromacs/fileio/warninp.h"
  45 #include "gromacs/math/units.h"
  46 #include "gromacs/math/utilities.h"
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/mdtypes/awh_params.h"
  49 #include "gromacs/mdtypes/inputrec.h"
  50 #include "gromacs/mdtypes/md_enums.h"
  51 #include "gromacs/mdtypes/multipletimestepping.h"
  52 #include "gromacs/mdtypes/pull_params.h"
  53 #include "gromacs/pbcutil/pbc.h"
  54 #include "gromacs/pulling/pull.h"
  55 #include "gromacs/random/seed.h"
  56 #include "gromacs/topology/mtop_util.h"
  57 #include "gromacs/utility/cstringutil.h"
  58 #include "gromacs/utility/enumerationhelpers.h"
  59 #include "gromacs/utility/fatalerror.h"
  60 #include "gromacs/utility/smalloc.h"
  61 #include "gromacs/utility/stringutil.h"
  62
  63 #include "biasparams.h"
  64 #include "biassharing.h"
  65
  66 namespace gmx
  67 {
  68
  69 const char* enumValueToString(AwhTargetType enumValue)
  70 {
  71     constexpr gmx::EnumerationArray<AwhTargetType, const char*> awhTargetTypeNames = {
  72         "constant", "cutoff", "boltzmann", "local-boltzmann"
  73     };
  74     return awhTargetTypeNames[enumValue];
  75 }
  76
  77 const char* enumValueToString(AwhHistogramGrowthType enumValue)
  78 {
  79     constexpr gmx::EnumerationArray<AwhHistogramGrowthType, const char*> awhHistogramGrowthTypeNames = {
  80         "exp-linear", "linear"
  81     };
  82     return awhHistogramGrowthTypeNames[enumValue];
  83 }
  84
  85 const char* enumValueToString(AwhPotentialType enumValue)
  86 {
  87     constexpr gmx::EnumerationArray<AwhPotentialType, const char*> awhPotentialTypeNames = {
  88         "convolved", "umbrella"
  89     };
  90     return awhPotentialTypeNames[enumValue];
  91 }
  92
  93 const char* enumValueToString(AwhCoordinateProviderType enumValue)
  94 {
  95     constexpr gmx::EnumerationArray<AwhCoordinateProviderType, const char*> awhCoordinateProviderTypeNames = {
  96         "pull", "fep-lambda"
  97     };
  98     return awhCoordinateProviderTypeNames[enumValue];
  99 }
 100
 101 namespace
 102 {
 103
 104 /*! \brief
 105  * Check multiple time-stepping consistency between AWH and pull and/or FEP
 106  *
 107  * \param[in] inputrec  Inputput parameter struct.
 108  * \param[in,out] wi    Struct for bookeeping warnings.
 109  */
 110 void checkMtsConsistency(const t_inputrec& inputrec, warninp_t wi)
 111 {
 112     if (!inputrec.useMts)
 113     {
 114         return;
 115     }
 116
 117     GMX_RELEASE_ASSERT(inputrec.mtsLevels.size() == 2, "Only 2 MTS levels supported here");
 118
 119     bool usesPull = false;
 120     bool usesFep  = false;
 121     for (int b = 0; b < inputrec.awhParams->numBias; b++)
 122     {
 123         const auto& biasParams = inputrec.awhParams->awhBiasParams[b];
 124         for (int d = 0; d < biasParams.ndim; d++)
 125         {
 126             switch (biasParams.dimParams[d].eCoordProvider)
 127             {
 128                 case AwhCoordinateProviderType::Pull: usesPull = true; break;
 129                 case AwhCoordinateProviderType::FreeEnergyLambda: usesFep = true; break;
 130                 default: GMX_RELEASE_ASSERT(false, "Unsupported coord provider");
 131             }
 132         }
 133     }
 134     const int awhMtsLevel = forceGroupMtsLevel(inputrec.mtsLevels, MtsForceGroups::Awh);
 135     if (usesPull && forceGroupMtsLevel(inputrec.mtsLevels, MtsForceGroups::Pull) != awhMtsLevel)
 136     {
 137         warning_error(wi,
 138                       "When AWH is applied to pull coordinates, pull and AWH should be computed at "
 139                       "the same MTS level");
 140     }
 141     if (usesFep && awhMtsLevel != ssize(inputrec.mtsLevels) - 1)
 142     {
 143         warning_error(wi,
 144                       "When AWH is applied to the free-energy lambda with MTS, AWH should be "
 145                       "computed at the slow MTS level");
 146     }
 147
 148     if (inputrec.awhParams->nstSampleCoord % inputrec.mtsLevels[awhMtsLevel].stepFactor != 0)
 149     {
 150         warning_error(wi,
 151                       "With MTS applied to AWH, awh-nstsample should be a multiple of mts-factor");
 152     }
 153 }
 154
 155 /*! \brief
 156  * Check the parameters of an AWH bias pull dimension.
 157  *
 158  * \param[in] prefix         Prefix for dimension parameters.
 159  * \param[in,out] dimParams  AWH dimensional parameters.
 160  * \param[in] pull_params    Pull parameters.
 161  * \param[in,out] wi         Struct for bookeeping warnings.
 162  */
 163 void checkPullDimParams(const std::string&   prefix,
 164                         AwhDimParams*        dimParams,
 165                         const pull_params_t& pull_params,
 166                         warninp_t            wi)
 167 {
 168     if (dimParams->coordIndex < 0)
 169     {
 170         gmx_fatal(FARGS,
 171                   "Failed to read a valid coordinate index for %s-coord-index. "
 172                   "Note that the pull coordinate indexing starts at 1.",
 173                   prefix.c_str());
 174     }
 175     if (dimParams->coordIndex >= pull_params.ncoord)
 176     {
 177         gmx_fatal(FARGS,
 178                   "The given AWH coordinate index (%d) is larger than the number of pull "
 179                   "coordinates (%d)",
 180                   dimParams->coordIndex + 1,
 181                   pull_params.ncoord);
 182     }
 183     if (pull_params.coord[dimParams->coordIndex].rate != 0)
 184     {
 185         auto message = formatString(
 186                 "Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate",
 187                 dimParams->coordIndex + 1,
 188                 pull_params.coord[dimParams->coordIndex].rate);
 189         warning_error(wi, message);
 190     }
 191
 192     if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
 193     {
 194         auto message = formatString(
 195                 "The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
 196                 "This will result in only one point along this axis in the coordinate value grid.",
 197                 prefix.c_str(),
 198                 dimParams->origin,
 199                 prefix.c_str(),
 200                 dimParams->end);
 201         warning(wi, message);
 202     }
 203
 204     if (dimParams->forceConstant <= 0)
 205     {
 206         warning_error(wi, "The force AWH bias force constant should be > 0");
 207     }
 208
 209     /* Grid params for each axis */
 210     PullGroupGeometry eGeom = pull_params.coord[dimParams->coordIndex].eGeom;
 211
 212     /* Check that the requested interval is in allowed range */
 213     if (eGeom == PullGroupGeometry::Distance)
 214     {
 215         if (dimParams->origin < 0 || dimParams->end < 0)
 216         {
 217             gmx_fatal(FARGS,
 218                       "%s-start (%g) or %s-end (%g) set to a negative value. With pull "
 219                       "geometry distance coordinate values are non-negative. "
 220                       "Perhaps you want to use geometry %s instead?",
 221                       prefix.c_str(),
 222                       dimParams->origin,
 223                       prefix.c_str(),
 224                       dimParams->end,
 225                       enumValueToString(PullGroupGeometry::Direction));
 226         }
 227     }
 228     else if (eGeom == PullGroupGeometry::Angle || eGeom == PullGroupGeometry::AngleAxis)
 229     {
 230         if (dimParams->origin < 0 || dimParams->end > 180)
 231         {
 232             gmx_fatal(FARGS,
 233                       "%s-start (%g) and %s-end (%g) are outside of the allowed range "
 234                       "0 to 180 deg for pull geometries %s and %s ",
 235                       prefix.c_str(),
 236                       dimParams->origin,
 237                       prefix.c_str(),
 238                       dimParams->end,
 239                       enumValueToString(PullGroupGeometry::Angle),
 240                       enumValueToString(PullGroupGeometry::AngleAxis));
 241         }
 242     }
 243     else if (eGeom == PullGroupGeometry::Dihedral)
 244     {
 245         if (dimParams->origin < -180 || dimParams->end > 180)
 246         {
 247             gmx_fatal(FARGS,
 248                       "%s-start (%g) and %s-end (%g) are outside of the allowed range "
 249                       "-180 to 180 deg for pull geometry %s. ",
 250                       prefix.c_str(),
 251                       dimParams->origin,
 252                       prefix.c_str(),
 253                       dimParams->end,
 254                       enumValueToString(PullGroupGeometry::Dihedral));
 255         }
 256     }
 257 }
 258
 259 /*! \brief
 260  * Check parameters of an AWH bias in a free energy lambda state dimension.
 261  *
 262  * \param[in] prefix         Prefix for dimension parameters.
 263  * \param[in,out] dimParams  AWH dimensional parameters.
 264  * \param[in] lambdaParams   The free energy lambda related parameters.
 265  * \param[in] efep           This is the type of FEP calculation (efep enumerator).
 266  * \param[in,out] wi         Struct for bookeeping warnings.
 267  */
 268 void checkFepLambdaDimParams(const std::string&               prefix,
 269                              const AwhDimParams*              dimParams,
 270                              const t_lambda*                  lambdaParams,
 271                              const FreeEnergyPerturbationType efep,
 272                              warninp_t                        wi)
 273 {
 274     std::string opt;
 275
 276     if (!lambdaParams)
 277     {
 278         gmx_fatal(FARGS,
 279                   "There must be free energy input if using AWH to steer the free energy lambda "
 280                   "state.");
 281     }
 282
 283     if (lambdaParams->lambda_neighbors != -1)
 284     {
 285         gmx_fatal(FARGS,
 286                   "When running AWH coupled to the free energy lambda state all lambda states "
 287                   "should be used as neighbors in order to get correct probabilities, i.e. "
 288                   "calc-lambda-neighbors (%d) must be %d.",
 289                   lambdaParams->lambda_neighbors,
 290                   -1);
 291     }
 292
 293     if (efep == FreeEnergyPerturbationType::SlowGrowth || lambdaParams->delta_lambda != 0)
 294     {
 295         gmx_fatal(FARGS,
 296                   "AWH coupled to the free energy lambda state is not compatible with slow-growth "
 297                   "and delta-lambda must be 0.");
 298     }
 299
 300     if (efep == FreeEnergyPerturbationType::Expanded)
 301     {
 302         gmx_fatal(FARGS,
 303                   "AWH is not treated like other expanded ensemble methods. Do not use expanded.");
 304     }
 305
 306     if (dimParams->origin < 0)
 307     {
 308         opt = prefix + "-start";
 309         gmx_fatal(FARGS,
 310                   "When running AWH coupled to the free energy lambda state the lower lambda state "
 311                   "for AWH, %s (%.0f), must be >= 0.",
 312                   opt.c_str(),
 313                   dimParams->origin);
 314     }
 315     if (dimParams->end >= lambdaParams->n_lambda)
 316     {
 317         opt = prefix + "-end";
 318         gmx_fatal(FARGS,
 319                   "When running AWH coupled to the free energy lambda state the upper lambda state "
 320                   "for AWH, %s (%.0f), must be < n_lambda (%d).",
 321                   opt.c_str(),
 322                   dimParams->origin,
 323                   lambdaParams->n_lambda);
 324     }
 325     if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
 326     {
 327         auto message = formatString(
 328                 "The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
 329                 "This will result in only one lambda point along this free energy lambda state "
 330                 "axis in the coordinate value grid.",
 331                 prefix.c_str(),
 332                 dimParams->origin,
 333                 prefix.c_str(),
 334                 dimParams->end);
 335         warning(wi, message);
 336     }
 337
 338     if (dimParams->forceConstant != 0)
 339     {
 340         warning_error(
 341                 wi,
 342                 "The force AWH bias force constant is not used with free energy lambda state as "
 343                 "coordinate provider.");
 344     }
 345 }
 346
 347 /*! \brief
 348  * Check that AWH FEP is not combined with incompatible decoupling.
 349  *
 350  * \param[in] mtop      System topology.
 351  * \param[in,out] wi    Struct for bookeeping warnings.
 352  */
 353 void checkFepLambdaDimDecouplingConsistency(const gmx_mtop_t& mtop, warninp_t wi)
 354 {
 355     if (haveFepPerturbedMasses(mtop))
 356     {
 357         warning(wi,
 358                 "Masses may not be perturbed when using the free energy lambda state as AWH "
 359                 "coordinate provider. If you are using fep-lambdas to specify lambda states "
 360                 "make sure that you also specify mass-lambdas without perturbation.");
 361     }
 362     if (havePerturbedConstraints(mtop))
 363     {
 364         warning(wi,
 365                 "Constraints may not be perturbed when using the free energy lambda state as AWH "
 366                 "coordinate provider. If you are using fep-lambdas to specify lambda states "
 367                 "make sure that you also specify mass-lambdas without perturbation.");
 368     }
 369 }
 370
 371 /*! \brief
 372  * Read parameters of an AWH bias dimension.
 373  *
 374  * \param[in,out] inp        Input file entries.
 375  * \param[in] prefix         Prefix for dimension parameters.
 376  * \param[in,out] dimParams  AWH dimensional parameters.
 377  * \param[in,out] wi         Struct for bookeeping warnings.
 378  * \param[in] bComment       True if comments should be printed.
 379  */
 380 void readDimParams(std::vector<t_inpfile>* inp,
 381                    const std::string&      prefix,
 382                    AwhDimParams*           dimParams,
 383                    warninp_t               wi,
 384                    bool                    bComment)
 385 {
 386     std::string opt;
 387     if (bComment)
 388     {
 389         printStringNoNewline(
 390                 inp,
 391                 "The provider of the reaction coordinate, "
 392                 "currently only 'pull' and 'fep-lambda' (free energy lambda state) is supported");
 393     }
 394     opt                       = prefix + "-coord-provider";
 395     dimParams->eCoordProvider = getEnum<AwhCoordinateProviderType>(inp, opt.c_str(), wi);
 396
 397     if (bComment)
 398     {
 399         printStringNoNewline(inp, "The coordinate index for this dimension");
 400     }
 401     opt = prefix + "-coord-index";
 402     int coordIndexInput;
 403     coordIndexInput = get_eint(inp, opt, 1, wi);
 404
 405     /* The pull coordinate indices start at 1 in the input file, at 0 internally */
 406     dimParams->coordIndex = coordIndexInput - 1;
 407
 408     if (bComment)
 409     {
 410         printStringNoNewline(inp, "Start and end values for each coordinate dimension");
 411     }
 412     opt               = prefix + "-start";
 413     dimParams->origin = get_ereal(inp, opt, 0., wi);
 414     opt               = prefix + "-end";
 415     dimParams->end    = get_ereal(inp, opt, 0., wi);
 416
 417     if (bComment)
 418     {
 419         printStringNoNewline(
 420                 inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
 421     }
 422     opt                      = prefix + "-force-constant";
 423     dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
 424
 425     if (bComment)
 426     {
 427         printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps, rad^2/ps or ps^-1)");
 428     }
 429     opt                  = prefix + "-diffusion";
 430     dimParams->diffusion = get_ereal(inp, opt, 0, wi);
 431
 432     if (bComment)
 433     {
 434         printStringNoNewline(inp,
 435                              "Diameter that needs to be sampled around a point before it is "
 436                              "considered covered. In FEP dimensions the cover diameter is "
 437                              "specified in lambda states.");
 438     }
 439     opt                      = prefix + "-cover-diameter";
 440     dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
 441 }
 442
 443 /*! \brief
 444  * Check the parameters of an AWH bias dimension.
 445  *
 446  * \param[in] prefix         Prefix for dimension parameters.
 447  * \param[in,out] dimParams  AWH dimensional parameters.
 448  * \param[in] ir             Input parameter struct.
 449  * \param[in,out] wi         Struct for bookeeping warnings.
 450  */
 451 void checkDimParams(const std::string& prefix, AwhDimParams* dimParams, const t_inputrec* ir, warninp_t wi)
 452 {
 453     if (dimParams->eCoordProvider == AwhCoordinateProviderType::Pull)
 454     {
 455         if (!ir->bPull)
 456         {
 457             gmx_fatal(FARGS,
 458                       "AWH biasing along a pull dimension is only compatible with COM pulling "
 459                       "turned on");
 460         }
 461         checkPullDimParams(prefix, dimParams, *ir->pull, wi);
 462     }
 463     else if (dimParams->eCoordProvider == AwhCoordinateProviderType::FreeEnergyLambda)
 464     {
 465         if (ir->efep == FreeEnergyPerturbationType::No)
 466         {
 467             gmx_fatal(FARGS,
 468                       "AWH biasing along a free energy lambda state dimension is only compatible "
 469                       "with free energy turned on");
 470         }
 471         checkFepLambdaDimParams(prefix, dimParams, ir->fepvals.get(), ir->efep, wi);
 472     }
 473     else
 474     {
 475         gmx_fatal(FARGS,
 476                   "AWH biasing can only be  applied to pull and free energy lambda state "
 477                   "coordinates");
 478     }
 479 }
 480
 481 /*! \brief
 482  * Check consistency of input at the AWH bias level.
 483  *
 484  * \param[in]     awhBiasParams  AWH bias parameters.
 485  * \param[in,out] wi             Struct for bookkeeping warnings.
 486  */
 487 void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
 488 {
 489     /* Covering diameter and sharing warning. */
 490     for (int d = 0; d < awhBiasParams.ndim; d++)
 491     {
 492         double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
 493         if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
 494         {
 495             warning(wi,
 496                     "The covering diameter is only relevant to set for bias sharing simulations.");
 497         }
 498     }
 499 }
 500
 501 /*! \brief
 502  * Read parameters of an AWH bias.
 503  *
 504  * \param[in,out] inp            Input file entries.
 505  * \param[in,out] awhBiasParams  AWH dimensional parameters.
 506  * \param[in]     prefix         Prefix for bias parameters.
 507  * \param[in,out] wi             Struct for bookeeping warnings.
 508  * \param[in]     bComment       True if comments should be printed.
 509  */
 510 void readBiasParams(std::vector<t_inpfile>* inp,
 511                     AwhBiasParams*          awhBiasParams,
 512                     const std::string&      prefix,
 513                     warninp_t               wi,
 514                     bool                    bComment)
 515 {
 516     if (bComment)
 517     {
 518         printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
 519     }
 520
 521     std::string opt             = prefix + "-error-init";
 522     awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi);
 523
 524     if (bComment)
 525     {
 526         printStringNoNewline(inp,
 527                              "Growth rate of the reference histogram determining the bias update "
 528                              "size: exp-linear or linear");
 529     }
 530     opt                    = prefix + "-growth";
 531     awhBiasParams->eGrowth = getEnum<AwhHistogramGrowthType>(inp, opt.c_str(), wi);
 532
 533     if (bComment)
 534     {
 535         printStringNoNewline(inp,
 536                              "Start the simulation by equilibrating histogram towards the target "
 537                              "distribution: no or yes");
 538     }
 539     opt                                 = prefix + "-equilibrate-histogram";
 540     awhBiasParams->equilibrateHistogram = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
 541
 542     if (bComment)
 543     {
 544         printStringNoNewline(
 545                 inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
 546     }
 547     opt                    = prefix + "-target";
 548     awhBiasParams->eTarget = getEnum<AwhTargetType>(inp, opt.c_str(), wi);
 549
 550     if (bComment)
 551     {
 552         printStringNoNewline(inp,
 553                              "Boltzmann beta scaling factor for target distribution types "
 554                              "'boltzmann' and 'boltzmann-local'");
 555     }
 556     opt                              = prefix + "-target-beta-scaling";
 557     awhBiasParams->targetBetaScaling = get_ereal(inp, opt, 0, wi);
 558
 559     if (bComment)
 560     {
 561         printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
 562     }
 563     opt                         = prefix + "-target-cutoff";
 564     awhBiasParams->targetCutoff = get_ereal(inp, opt, 0, wi);
 565
 566     if (bComment)
 567     {
 568         printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
 569     }
 570     opt                      = prefix + "-user-data";
 571     awhBiasParams->bUserData = getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No;
 572
 573     if (bComment)
 574     {
 575         printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
 576     }
 577     opt                       = prefix + "-share-group";
 578     awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
 579
 580     if (bComment)
 581     {
 582         printStringNoNewline(inp, "Dimensionality of the coordinate");
 583     }
 584     opt                 = prefix + "-ndim";
 585     awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
 586
 587     /* Check this before starting to read the AWH dimension parameters. */
 588     if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
 589     {
 590         gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), awhBiasParams->ndim, c_biasMaxNumDim);
 591     }
 592     snew(awhBiasParams->dimParams, awhBiasParams->ndim);
 593     for (int d = 0; d < awhBiasParams->ndim; d++)
 594     {
 595         bComment              = bComment && d == 0;
 596         std::string prefixdim = prefix + formatString("-dim%d", d + 1);
 597         readDimParams(inp, prefixdim, &awhBiasParams->dimParams[d], wi, bComment);
 598     }
 599 }
 600
 601 /*! \brief
 602  * Check the parameters of an AWH bias.
 603  *
 604  * \param[in]     awhBiasParams  AWH dimensional parameters.
 605  * \param[in]     prefix         Prefix for bias parameters.
 606  * \param[in]     ir             Input parameter struct.
 607  * \param[in,out] wi             Struct for bookeeping warnings.
 608  */
 609 void checkBiasParams(const AwhBiasParams* awhBiasParams, const std::string& prefix, const t_inputrec* ir, warninp_t wi)
 610 {
 611     std::string opt = prefix + "-error-init";
 612     if (awhBiasParams->errorInitial <= 0)
 613     {
 614         gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
 615     }
 616
 617     opt = prefix + "-equilibrate-histogram";
 618     if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != AwhHistogramGrowthType::ExponentialLinear)
 619     {
 620         auto message =
 621                 formatString("Option %s will only have an effect for histogram growth type '%s'.",
 622                              opt.c_str(),
 623                              enumValueToString(AwhHistogramGrowthType::ExponentialLinear));
 624         warning(wi, message);
 625     }
 626
 627     if ((awhBiasParams->eTarget == AwhTargetType::LocalBoltzmann)
 628         && (awhBiasParams->eGrowth == AwhHistogramGrowthType::ExponentialLinear))
 629     {
 630         auto message = formatString(
 631                 "Target type '%s' combined with histogram growth type '%s' is not "
 632                 "expected to give stable bias updates. You probably want to use growth type "
 633                 "'%s' instead.",
 634                 enumValueToString(AwhTargetType::LocalBoltzmann),
 635                 enumValueToString(AwhHistogramGrowthType::ExponentialLinear),
 636                 enumValueToString(AwhHistogramGrowthType::Linear));
 637         warning(wi, message);
 638     }
 639
 640     opt = prefix + "-target-beta-scaling";
 641     switch (awhBiasParams->eTarget)
 642     {
 643         case AwhTargetType::Boltzmann:
 644         case AwhTargetType::LocalBoltzmann:
 645             if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
 646             {
 647                 gmx_fatal(FARGS,
 648                           "%s = %g is not useful for target type %s.",
 649                           opt.c_str(),
 650                           awhBiasParams->targetBetaScaling,
 651                           enumValueToString(awhBiasParams->eTarget));
 652             }
 653             break;
 654         default:
 655             if (awhBiasParams->targetBetaScaling != 0)
 656             {
 657                 gmx_fatal(
 658                         FARGS,
 659                         "Value for %s (%g) set explicitly but will not be used for target type %s.",
 660                         opt.c_str(),
 661                         awhBiasParams->targetBetaScaling,
 662                         enumValueToString(awhBiasParams->eTarget));
 663             }
 664             break;
 665     }
 666
 667     opt = prefix + "-target-cutoff";
 668     switch (awhBiasParams->eTarget)
 669     {
 670         case AwhTargetType::Cutoff:
 671             if (awhBiasParams->targetCutoff <= 0)
 672             {
 673                 gmx_fatal(FARGS,
 674                           "%s = %g is not useful for target type %s.",
 675                           opt.c_str(),
 676                           awhBiasParams->targetCutoff,
 677                           enumValueToString(awhBiasParams->eTarget));
 678             }
 679             break;
 680         default:
 681             if (awhBiasParams->targetCutoff != 0)
 682             {
 683                 gmx_fatal(
 684                         FARGS,
 685                         "Value for %s (%g) set explicitly but will not be used for target type %s.",
 686                         opt.c_str(),
 687                         awhBiasParams->targetCutoff,
 688                         enumValueToString(awhBiasParams->eTarget));
 689             }
 690             break;
 691     }
 692
 693     opt = prefix + "-share-group";
 694     if (awhBiasParams->shareGroup < 0)
 695     {
 696         warning_error(wi, "AWH bias share-group should be >= 0");
 697     }
 698
 699     opt = prefix + "-ndim";
 700     if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
 701     {
 702         gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), awhBiasParams->ndim, c_biasMaxNumDim);
 703     }
 704     if (awhBiasParams->ndim > 2)
 705     {
 706         warning_note(wi,
 707                      "For awh-dim > 2 the estimate based on the diffusion and the initial error is "
 708                      "currently only a rough guideline."
 709                      " You should verify its usefulness for your system before production runs!");
 710     }
 711     for (int d = 0; d < awhBiasParams->ndim; d++)
 712     {
 713         std::string prefixdim = prefix + formatString("-dim%d", d + 1);
 714         checkDimParams(prefixdim, &awhBiasParams->dimParams[d], ir, wi);
 715     }
 716
 717     /* Check consistencies here that cannot be checked at read time at a lower level. */
 718     checkInputConsistencyAwhBias(*awhBiasParams, wi);
 719 }
 720
 721 /*! \brief
 722  * Check consistency of input at the AWH level.
 723  *
 724  * \param[in]     awhParams  AWH parameters.
 725  * \param[in,out] wi         Struct for bookkeeping warnings.
 726  */
 727 void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi)
 728 {
 729     /* Each pull coord can map to at most 1 AWH coord.
 730      * Check that we have a shared bias when requesting multisim sharing.
 731      */
 732     bool haveSharedBias = false;
 733     for (int k1 = 0; k1 < awhParams.numBias; k1++)
 734     {
 735         const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1];
 736
 737         if (awhBiasParams1.shareGroup > 0)
 738         {
 739             haveSharedBias = true;
 740         }
 741
 742         /* k1 is the reference AWH, k2 is the AWH we compare with (can be equal to k1) */
 743         for (int k2 = k1; k2 < awhParams.numBias; k2++)
 744         {
 745             for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
 746             {
 747                 if (awhBiasParams1.dimParams[d1].eCoordProvider == AwhCoordinateProviderType::FreeEnergyLambda)
 748                 {
 749                     continue;
 750                 }
 751                 const AwhBiasParams& awhBiasParams2 = awhParams.awhBiasParams[k2];
 752
 753                 /* d1 is the reference dimension of the reference AWH. d2 is the dim index of the AWH to compare with. */
 754                 for (int d2 = 0; d2 < awhBiasParams2.ndim; d2++)
 755                 {
 756                     if (awhBiasParams2.dimParams[d2].eCoordProvider == AwhCoordinateProviderType::FreeEnergyLambda)
 757                     {
 758                         continue;
 759                     }
 760                     /* Give an error if (d1, k1) is different from (d2, k2) but the pull coordinate is the same */
 761                     if ((d1 != d2 || k1 != k2)
 762                         && (awhBiasParams1.dimParams[d1].coordIndex == awhBiasParams2.dimParams[d2].coordIndex))
 763                     {
 764                         char errormsg[STRLEN];
 765                         sprintf(errormsg,
 766                                 "One pull coordinate (%d) cannot be mapped to two separate AWH "
 767                                 "dimensions (awh%d-dim%d and awh%d-dim%d). "
 768                                 "If this is really what you want to do you will have to duplicate "
 769                                 "this pull coordinate.",
 770                                 awhBiasParams1.dimParams[d1].coordIndex + 1,
 771                                 k1 + 1,
 772                                 d1 + 1,
 773                                 k2 + 1,
 774                                 d2 + 1);
 775                         gmx_fatal(FARGS, "%s", errormsg);
 776                     }
 777                 }
 778             }
 779         }
 780     }
 781
 782     if (awhParams.shareBiasMultisim && !haveSharedBias)
 783     {
 784         warning(wi,
 785                 "Sharing of biases over multiple simulations is requested, but no bias is marked "
 786                 "as shared (share-group > 0)");
 787     }
 788
 789     /* mdrun does not support this (yet), but will check again */
 790     if (haveBiasSharingWithinSimulation(awhParams))
 791     {
 792         warning(wi,
 793                 "You have shared biases within a single simulation, but mdrun does not support "
 794                 "this (yet)");
 795     }
 796 }
 797 } // namespace
 798
 799 AwhParams* readAwhParams(std::vector<t_inpfile>* inp, warninp_t wi)
 800 {
 801     AwhParams* awhParams;
 802     snew(awhParams, 1);
 803     std::string opt;
 804
 805     /* Parameters common for all biases */
 806
 807     printStringNoNewline(inp, "The way to apply the biasing potential: convolved or umbrella");
 808     opt                   = "awh-potential";
 809     awhParams->ePotential = getEnum<AwhPotentialType>(inp, opt.c_str(), wi);
 810
 811     printStringNoNewline(inp,
 812                          "The random seed used for sampling the umbrella center in the case of "
 813                          "umbrella type potential");
 814     opt             = "awh-seed";
 815     awhParams->seed = get_eint(inp, opt, -1, wi);
 816     if (awhParams->seed == -1)
 817     {
 818         awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
 819         fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
 820     }
 821
 822     printStringNoNewline(inp, "Data output interval in number of steps");
 823     opt               = "awh-nstout";
 824     awhParams->nstOut = get_eint(inp, opt, 100000, wi);
 825
 826     printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
 827     opt                       = "awh-nstsample";
 828     awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
 829
 830     printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
 831     opt                                   = "awh-nsamples-update";
 832     awhParams->numSamplesUpdateFreeEnergy = get_eint(inp, opt, 10, wi);
 833
 834     printStringNoNewline(
 835             inp, "When true, biases with share-group>0 are shared between multiple simulations");
 836     opt                          = "awh-share-multisim";
 837     awhParams->shareBiasMultisim = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
 838
 839     printStringNoNewline(inp, "The number of independent AWH biases");
 840     opt                = "awh-nbias";
 841     awhParams->numBias = get_eint(inp, opt, 1, wi);
 842     /* Check this before starting to read the AWH biases. */
 843     if (awhParams->numBias <= 0)
 844     {
 845         gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
 846     }
 847
 848     /* Read the parameters specific to each AWH bias */
 849     snew(awhParams->awhBiasParams, awhParams->numBias);
 850
 851     for (int k = 0; k < awhParams->numBias; k++)
 852     {
 853         bool        bComment  = (k == 0);
 854         std::string prefixawh = formatString("awh%d", k + 1);
 855         readBiasParams(inp, &awhParams->awhBiasParams[k], prefixawh, wi, bComment);
 856     }
 857
 858     return awhParams;
 859 }
 860
 861 void checkAwhParams(const AwhParams* awhParams, const t_inputrec* ir, warninp_t wi)
 862 {
 863     std::string opt;
 864
 865     checkMtsConsistency(*ir, wi);
 866
 867     opt = "awh-nstout";
 868     if (awhParams->nstOut <= 0)
 869     {
 870         auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
 871                                     opt.c_str(),
 872                                     awhParams->nstOut);
 873         warning_error(wi, message);
 874     }
 875     /* This restriction can be removed by changing a flag of print_ebin() */
 876     if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
 877     {
 878         auto message = formatString(
 879                 "%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(), awhParams->nstOut, ir->nstenergy);
 880         warning_error(wi, message);
 881     }
 882
 883     opt = "awh-nsamples-update";
 884     if (awhParams->numSamplesUpdateFreeEnergy <= 0)
 885     {
 886         warning_error(wi, opt + " needs to be an integer > 0");
 887     }
 888
 889     for (int k = 0; k < awhParams->numBias; k++)
 890     {
 891         std::string prefixawh = formatString("awh%d", k + 1);
 892         checkBiasParams(&awhParams->awhBiasParams[k], prefixawh, ir, wi);
 893     }
 894
 895     /* Do a final consistency check before returning */
 896     checkInputConsistencyAwh(*awhParams, wi);
 897
 898     if (ir->init_step != 0)
 899     {
 900         warning_error(wi, "With AWH init-step should be 0");
 901     }
 902 }
 903
 904 /*! \brief
 905  * Gets the period of a pull coordinate.
 906  *
 907  * \param[in] pullCoordParams   The parameters for the pull coordinate.
 908  * \param[in] pbc               The PBC setup
 909  * \param[in] intervalLength    The length of the AWH interval for this pull coordinate
 910  * \returns the period (or 0 if not periodic).
 911  */
 912 static double get_pull_coord_period(const t_pull_coord& pullCoordParams, const t_pbc& pbc, const real intervalLength)
 913 {
 914     double period = 0;
 915
 916     if (pullCoordParams.eGeom == PullGroupGeometry::Direction)
 917     {
 918         const real margin = 0.001;
 919         // Make dims periodic when the interval covers > 95%
 920         const real periodicFraction = 0.95;
 921
 922         // Check if the pull direction is along a box vector
 923         for (int dim = 0; dim < pbc.ndim_ePBC; dim++)
 924         {
 925             const real boxLength    = norm(pbc.box[dim]);
 926             const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]);
 927             if (innerProduct >= (1 - margin) * boxLength && innerProduct <= (1 + margin) * boxLength)
 928             {
 929                 if (intervalLength > (1 + margin) * boxLength)
 930                 {
 931                     gmx_fatal(FARGS,
 932                               "The AWH interval (%f nm) for a pull coordinate is larger than the "
 933                               "box size (%f nm)",
 934                               intervalLength,
 935                               boxLength);
 936                 }
 937
 938                 if (intervalLength > periodicFraction * boxLength)
 939                 {
 940                     period = boxLength;
 941                 }
 942             }
 943         }
 944     }
 945     else if (pullCoordParams.eGeom == PullGroupGeometry::Dihedral)
 946     {
 947         /* The dihedral angle is periodic in -180 to 180 deg */
 948         period = 360;
 949     }
 950
 951     return period;
 952 }
 953
 954 /*! \brief
 955  * Checks if the given interval is defined in the correct periodic interval.
 956  *
 957  * \param[in] origin      Start value of interval.
 958  * \param[in] end         End value of interval.
 959  * \param[in] period      Period (or 0 if not periodic).
 960  * \returns true if the end point values are in the correct periodic interval.
 961  */
 962 static bool intervalIsInPeriodicInterval(double origin, double end, double period)
 963 {
 964     return (period == 0) || (std::fabs(origin) <= 0.5 * period && std::fabs(end) <= 0.5 * period);
 965 }
 966
 967 /*! \brief
 968  * Checks if a value is within an interval.
 969  *
 970  * \param[in] origin      Start value of interval.
 971  * \param[in] end         End value of interval.
 972  * \param[in] period      Period (or 0 if not periodic).
 973  * \param[in] value       Value to check.
 974  * \returns true if the value is within the interval.
 975  */
 976 static bool valueIsInInterval(double origin, double end, double period, double value)
 977 {
 978     bool bIn_interval;
 979
 980     if (period > 0)
 981     {
 982         if (origin < end)
 983         {
 984             /* The interval closes within the periodic interval */
 985             bIn_interval = (value >= origin) && (value <= end);
 986         }
 987         else
 988         {
 989             /* The interval wraps around the periodic boundary */
 990             bIn_interval = ((value >= origin) && (value <= 0.5 * period))
 991                            || ((value >= -0.5 * period) && (value <= end));
 992         }
 993     }
 994     else
 995     {
 996         bIn_interval = (value >= origin) && (value <= end);
 997     }
 998
 999     return bIn_interval;
1000 }
1001
1002 /*! \brief
1003  * Check if the starting configuration is consistent with the given interval.
1004  *
1005  * \param[in]     awhParams  AWH parameters.
1006  * \param[in,out] wi         Struct for bookeeping warnings.
1007  */
1008 static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t wi)
1009 {
1010     for (int k = 0; k < awhParams->numBias; k++)
1011     {
1012         AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
1013         for (int d = 0; d < awhBiasParams->ndim; d++)
1014         {
1015             AwhDimParams* dimParams  = &awhBiasParams->dimParams[d];
1016             int           coordIndex = dimParams->coordIndex;
1017             double origin = dimParams->origin, end = dimParams->end, period = dimParams->period;
1018             double coordValueInit = dimParams->coordValueInit;
1019
1020             if ((period == 0) && (origin > end))
1021             {
1022                 gmx_fatal(FARGS,
1023                           "For the non-periodic pull coordinates awh%d-dim%d-start (%f) cannot be "
1024                           "larger than awh%d-dim%d-end (%f)",
1025                           k + 1,
1026                           d + 1,
1027                           origin,
1028                           k + 1,
1029                           d + 1,
1030                           end);
1031             }
1032
1033             /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
1034                Make sure the AWH interval is within this reference interval.
1035
1036                Note: we could fairly simply allow using a  more general interval (e.g. [x, x + period]) but it complicates
1037                things slightly and I don't see that there is a great need for it. It would also mean that the interval would
1038                depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180,
1039                independent of AWH parameters.
1040              */
1041             if (!intervalIsInPeriodicInterval(origin, end, period))
1042             {
1043                 gmx_fatal(FARGS,
1044                           "When using AWH with periodic pull coordinate geometries "
1045                           "awh%d-dim%d-start (%.8g) and "
1046                           "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take "
1047                           "values in between "
1048                           "minus half a period and plus half a period, i.e. in the interval [%.8g, "
1049                           "%.8g].",
1050                           k + 1,
1051                           d + 1,
1052                           origin,
1053                           k + 1,
1054                           d + 1,
1055                           end,
1056                           period,
1057                           -0.5 * period,
1058                           0.5 * period);
1059             }
1060
1061             /* Warn if the pull initial coordinate value is not in the grid */
1062             if (!valueIsInInterval(origin, end, period, coordValueInit))
1063             {
1064                 auto message = formatString(
1065                         "The initial coordinate value (%.8g) for pull coordinate index %d falls "
1066                         "outside "
1067                         "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end "
1068                         "(%.8g). "
1069                         "This can lead to large initial forces pulling the coordinate towards the "
1070                         "sampling interval.",
1071                         coordValueInit,
1072                         coordIndex + 1,
1073                         k + 1,
1074                         d + 1,
1075                         origin,
1076                         k + 1,
1077                         d + 1,
1078                         end);
1079                 warning(wi, message);
1080             }
1081         }
1082     }
1083 }
1084
1085 /*! \brief
1086  * Sets AWH parameters, for one AWH pull dimension.
1087  *
1088  * \param[in,out] dimParams          AWH dimension parameters.
1089  * \param[in]     biasIndex             The index of the bias containing this AWH pull dimension.
1090  * \param[in]     dimIndex              The index of this AWH pull dimension.
1091  * \param[in]     pull_params           Pull parameters.
1092  * \param[in,out] pull_work             Pull working struct to register AWH bias in.
1093  * \param[in]     pbc                   A pbc information structure.
1094  * \param[in]     compressibility       Compressibility matrix for pressure coupling,
1095  * pass all 0 without pressure coupling.
1096  * \param[in,out] wi                    Struct for bookeeping warnings.
1097  *
1098  */
1099 static void setStateDependentAwhPullDimParams(AwhDimParams*        dimParams,
1100                                               const int            biasIndex,
1101                                               const int            dimIndex,
1102                                               const pull_params_t* pull_params,
1103                                               pull_t*              pull_work,
1104                                               const t_pbc&         pbc,
1105                                               const tensor&        compressibility,
1106                                               warninp_t            wi)
1107 {
1108     const t_pull_coord& pullCoordParams = pull_params->coord[dimParams->coordIndex];
1109
1110     if (pullCoordParams.eGeom == PullGroupGeometry::DirectionPBC)
1111     {
1112         gmx_fatal(FARGS,
1113                   "AWH does not support pull geometry '%s'. "
1114                   "If the maximum distance between the groups is always "
1115                   "less than half the box size, "
1116                   "you can use geometry '%s' instead.",
1117                   enumValueToString(PullGroupGeometry::DirectionPBC),
1118                   enumValueToString(PullGroupGeometry::Direction));
1119     }
1120
1121     dimParams->period = get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin);
1122     // We would like to check for scaling, but we don't have the full inputrec available here
1123     if (dimParams->period > 0
1124         && !(pullCoordParams.eGeom == PullGroupGeometry::Angle
1125              || pullCoordParams.eGeom == PullGroupGeometry::Dihedral))
1126     {
1127         bool coordIsScaled = false;
1128         for (int d2 = 0; d2 < DIM; d2++)
1129         {
1130             if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0)
1131             {
1132                 coordIsScaled = true;
1133             }
1134         }
1135         if (coordIsScaled)
1136         {
1137             std::string mesg = gmx::formatString(
1138                     "AWH dimension %d of bias %d is periodic with pull geometry '%s', "
1139                     "while you should are applying pressure scaling to the "
1140                     "corresponding box vector, this is not supported.",
1141                     dimIndex + 1,
1142                     biasIndex + 1,
1143                     enumValueToString(pullCoordParams.eGeom));
1144             warning(wi, mesg.c_str());
1145         }
1146     }
1147
1148     /* The initial coordinate value, converted to external user units. */
1149     dimParams->coordValueInit = get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
1150
1151     dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(&pullCoordParams);
1152 }
1153
1154 void setStateDependentAwhParams(AwhParams*           awhParams,
1155                                 const pull_params_t& pull_params,
1156                                 pull_t*              pull_work,
1157                                 const matrix         box,
1158                                 PbcType              pbcType,
1159                                 const tensor&        compressibility,
1160                                 const t_grpopts*     inputrecGroupOptions,
1161                                 const real           initLambda,
1162                                 const gmx_mtop_t&    mtop,
1163                                 warninp_t            wi)
1164 {
1165     /* The temperature is not really state depenendent but is not known
1166      * when read_awhParams is called (in get ir).
1167      * It is known first after do_index has been called in grompp.cpp.
1168      */
1169     if (inputrecGroupOptions->ref_t == nullptr || inputrecGroupOptions->ref_t[0] <= 0)
1170     {
1171         gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
1172     }
1173     for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
1174     {
1175         if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
1176         {
1177             gmx_fatal(FARGS,
1178                       "AWH biasing is currently only supported for identical temperatures for all "
1179                       "temperature coupling groups");
1180         }
1181     }
1182
1183     t_pbc pbc;
1184     set_pbc(&pbc, pbcType, box);
1185
1186     for (int k = 0; k < awhParams->numBias; k++)
1187     {
1188         AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
1189         for (int d = 0; d < awhBiasParams->ndim; d++)
1190         {
1191             AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
1192             if (dimParams->eCoordProvider == AwhCoordinateProviderType::Pull)
1193             {
1194                 setStateDependentAwhPullDimParams(
1195                         dimParams, k, d, &pull_params, pull_work, pbc, compressibility, wi);
1196             }
1197             else
1198             {
1199                 dimParams->coordValueInit = initLambda;
1200                 checkFepLambdaDimDecouplingConsistency(mtop, wi);
1201             }
1202         }
1203     }
1204     checkInputConsistencyInterval(awhParams, wi);
1205
1206     /* Register AWH as external potential with pull (for AWH dimensions that use pull as
1207      * reaction coordinate provider) to check consistency. */
1208     Awh::registerAwhWithPull(*awhParams, pull_work);
1209 }
1210
1211 } // namespace gmx