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