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