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,
 630                            const std::string&      prefix,
 631                            const t_inputrec&       ir,
 632                            warninp_t               wi,
 633                            bool                    bComment)
 634 {
 635     std::string opt;
 636     if (bComment)
 637     {
 638         printStringNoNewline(
 639                 inp,
 640                 "The provider of the reaction coordinate, "
 641                 "currently only 'pull' and 'fep-lambda' (free energy lambda state) is supported");
 642     }
 643
 644     opt             = prefix + "-coord-provider";
 645     eCoordProvider_ = getEnum<AwhCoordinateProviderType>(inp, opt.c_str(), wi);
 646
 647     if (bComment)
 648     {
 649         printStringNoNewline(inp, "The coordinate index for this dimension");
 650     }
 651     opt = prefix + "-coord-index";
 652     int coordIndexInput;
 653     coordIndexInput = get_eint(inp, opt, 1, wi);
 654     if (coordIndexInput < 1)
 655     {
 656         gmx_fatal(FARGS,
 657                   "Failed to read a valid coordinate index for %s. "
 658                   "Note that the pull coordinate indexing starts at 1.",
 659                   opt.c_str());
 660     }
 661
 662     /* The pull coordinate indices start at 1 in the input file, at 0 internally */
 663     coordIndex_ = coordIndexInput - 1;
 664
 665     if (bComment)
 666     {
 667         printStringNoNewline(inp, "Start and end values for each coordinate dimension");
 668     }
 669
 670     opt     = prefix + "-start";
 671     origin_ = get_ereal(inp, opt, 0., wi);
 672
 673     opt  = prefix + "-end";
 674     end_ = get_ereal(inp, opt, 0., wi);
 675
 676     if (bComment)
 677     {
 678         printStringNoNewline(
 679                 inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
 680     }
 681     opt            = prefix + "-force-constant";
 682     forceConstant_ = get_ereal(inp, opt, 0, wi);
 683
 684     if (bComment)
 685     {
 686         printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps or ps^-1)");
 687     }
 688     opt                   = prefix + "-diffusion";
 689     double diffusionValue = get_ereal(inp, opt, 0, wi);
 690     if (diffusionValue <= 0)
 691     {
 692         const double diffusion_default = 1e-5;
 693         auto         message           = formatString(
 694                 "%s not explicitly set by user. You can choose to use a default "
 695                 "value (%g nm^2/ps or rad^2/ps) but this may very well be "
 696                 "non-optimal for your system!",
 697                 opt.c_str(),
 698                 diffusion_default);
 699         warning(wi, message);
 700         diffusionValue = diffusion_default;
 701     }
 702     diffusion_ = diffusionValue;
 703
 704     if (bComment)
 705     {
 706         printStringNoNewline(inp,
 707                              "Diameter that needs to be sampled around a point before it is "
 708                              "considered covered. In FEP dimensions the cover diameter is "
 709                              "specified in lambda states.");
 710     }
 711     opt            = prefix + "-cover-diameter";
 712     coverDiameter_ = get_ereal(inp, opt, 0, wi);
 713
 714     checkDimParams(prefix, *this, ir, wi);
 715 }
 716
 717 AwhDimParams::AwhDimParams(ISerializer* serializer)
 718 {
 719     GMX_RELEASE_ASSERT(serializer->reading(),
 720                        "Can not use writing serializer for creating datastructure");
 721     serializer->doEnumAsInt(&eCoordProvider_);
 722     serializer->doInt(&coordIndex_);
 723     serializer->doDouble(&origin_);
 724     serializer->doDouble(&end_);
 725     serializer->doDouble(&period_);
 726     serializer->doDouble(&forceConstant_);
 727     serializer->doDouble(&diffusion_);
 728     serializer->doDouble(&coordValueInit_);
 729     serializer->doDouble(&coverDiameter_);
 730 }
 731
 732 void AwhDimParams::serialize(ISerializer* serializer)
 733 {
 734     GMX_RELEASE_ASSERT(!serializer->reading(),
 735                        "Can not use reading serializer for writing datastructure");
 736     serializer->doEnumAsInt(&eCoordProvider_);
 737     serializer->doInt(&coordIndex_);
 738     serializer->doDouble(&origin_);
 739     serializer->doDouble(&end_);
 740     serializer->doDouble(&period_);
 741     serializer->doDouble(&forceConstant_);
 742     serializer->doDouble(&diffusion_);
 743     serializer->doDouble(&coordValueInit_);
 744     serializer->doDouble(&coverDiameter_);
 745 }
 746
 747 AwhBiasParams::AwhBiasParams(std::vector<t_inpfile>* inp,
 748                              const std::string&      prefix,
 749                              const t_inputrec&       ir,
 750                              warninp_t               wi,
 751                              bool                    bComment)
 752 {
 753     if (bComment)
 754     {
 755         printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
 756     }
 757
 758     std::string opt = prefix + "-error-init";
 759     errorInitial_   = get_ereal(inp, opt, 10, wi);
 760
 761     if (bComment)
 762     {
 763         printStringNoNewline(inp,
 764                              "Growth rate of the reference histogram determining the bias update "
 765                              "size: exp-linear or linear");
 766     }
 767     opt      = prefix + "-growth";
 768     eGrowth_ = getEnum<AwhHistogramGrowthType>(inp, opt.c_str(), wi);
 769
 770     if (bComment)
 771     {
 772         printStringNoNewline(inp,
 773                              "Start the simulation by equilibrating histogram towards the target "
 774                              "distribution: no or yes");
 775     }
 776     opt                   = prefix + "-equilibrate-histogram";
 777     equilibrateHistogram_ = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
 778
 779     if (bComment)
 780     {
 781         printStringNoNewline(
 782                 inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
 783     }
 784     opt      = prefix + "-target";
 785     eTarget_ = getEnum<AwhTargetType>(inp, opt.c_str(), wi);
 786
 787     if (bComment)
 788     {
 789         printStringNoNewline(inp,
 790                              "Boltzmann beta scaling factor for target distribution types "
 791                              "'boltzmann' and 'boltzmann-local'");
 792     }
 793     opt                = prefix + "-target-beta-scaling";
 794     targetBetaScaling_ = get_ereal(inp, opt, 0, wi);
 795
 796     if (bComment)
 797     {
 798         printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
 799     }
 800     opt           = prefix + "-target-cutoff";
 801     targetCutoff_ = get_ereal(inp, opt, 0, wi);
 802
 803     if (bComment)
 804     {
 805         printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
 806     }
 807     opt        = prefix + "-user-data";
 808     bUserData_ = getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No;
 809
 810     if (bComment)
 811     {
 812         printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
 813     }
 814     opt         = prefix + "-share-group";
 815     shareGroup_ = get_eint(inp, opt, 0, wi);
 816
 817     if (bComment)
 818     {
 819         printStringNoNewline(inp, "Dimensionality of the coordinate");
 820     }
 821     opt      = prefix + "-ndim";
 822     int ndim = get_eint(inp, opt, 0, wi);
 823
 824     /* Check this before starting to read the AWH dimension parameters. */
 825     if (ndim <= 0 || ndim > c_biasMaxNumDim)
 826     {
 827         gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), ndim, c_biasMaxNumDim);
 828     }
 829     for (int d = 0; d < ndim; d++)
 830     {
 831         bComment              = bComment && d == 0;
 832         std::string prefixdim = prefix + formatString("-dim%d", d + 1);
 833         dimParams_.emplace_back(inp, prefixdim, ir, wi, bComment);
 834     }
 835     checkBiasParams(*this, prefix, ir, wi);
 836 }
 837
 838 AwhBiasParams::AwhBiasParams(ISerializer* serializer)
 839 {
 840     GMX_RELEASE_ASSERT(serializer->reading(),
 841                        "Can not use writing serializer to create datastructure");
 842     serializer->doEnumAsInt(&eTarget_);
 843     serializer->doDouble(&targetBetaScaling_);
 844     serializer->doDouble(&targetCutoff_);
 845     serializer->doEnumAsInt(&eGrowth_);
 846     int temp = 0;
 847     serializer->doInt(&temp);
 848     bUserData_ = static_cast<bool>(temp);
 849     serializer->doDouble(&errorInitial_);
 850     int numDimensions = dimParams_.size();
 851     serializer->doInt(&numDimensions);
 852     serializer->doInt(&shareGroup_);
 853     serializer->doBool(&equilibrateHistogram_);
 854
 855     for (int k = 0; k < numDimensions; k++)
 856     {
 857         dimParams_.emplace_back(serializer);
 858     }
 859     /* Check consistencies here that cannot be checked at read time at a lower level. */
 860     checkInputConsistencyAwhBias(*this, nullptr);
 861 }
 862
 863 void AwhBiasParams::serialize(ISerializer* serializer)
 864 {
 865     GMX_RELEASE_ASSERT(!serializer->reading(),
 866                        "Can not use reading serializer to write datastructure");
 867     serializer->doEnumAsInt(&eTarget_);
 868     serializer->doDouble(&targetBetaScaling_);
 869     serializer->doDouble(&targetCutoff_);
 870     serializer->doEnumAsInt(&eGrowth_);
 871     int temp = static_cast<int>(bUserData_);
 872     serializer->doInt(&temp);
 873     serializer->doDouble(&errorInitial_);
 874     int numDimensions = ndim();
 875     serializer->doInt(&numDimensions);
 876     serializer->doInt(&shareGroup_);
 877     serializer->doBool(&equilibrateHistogram_);
 878
 879     for (int k = 0; k < numDimensions; k++)
 880     {
 881         dimParams_[k].serialize(serializer);
 882     }
 883 }
 884
 885 AwhParams::AwhParams(std::vector<t_inpfile>* inp, const t_inputrec& ir, warninp_t wi)
 886 {
 887     std::string opt;
 888
 889     /* Parameters common for all biases */
 890
 891     printStringNoNewline(inp, "The way to apply the biasing potential: convolved or umbrella");
 892     opt            = "awh-potential";
 893     potentialEnum_ = getEnum<AwhPotentialType>(inp, opt.c_str(), wi);
 894
 895     printStringNoNewline(inp,
 896                          "The random seed used for sampling the umbrella center in the case of "
 897                          "umbrella type potential");
 898     opt   = "awh-seed";
 899     seed_ = get_eint(inp, opt, -1, wi);
 900     if (seed_ == -1)
 901     {
 902         seed_ = static_cast<int>(gmx::makeRandomSeed());
 903         fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", seed_);
 904     }
 905
 906     printStringNoNewline(inp, "Data output interval in number of steps");
 907     opt     = "awh-nstout";
 908     nstOut_ = get_eint(inp, opt, 100000, wi);
 909
 910     printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
 911     opt             = "awh-nstsample";
 912     nstSampleCoord_ = get_eint(inp, opt, 10, wi);
 913
 914     printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
 915     opt                         = "awh-nsamples-update";
 916     numSamplesUpdateFreeEnergy_ = get_eint(inp, opt, 10, wi);
 917
 918     printStringNoNewline(
 919             inp, "When true, biases with share-group>0 are shared between multiple simulations");
 920     opt                = "awh-share-multisim";
 921     shareBiasMultisim_ = (getEnum<Boolean>(inp, opt.c_str(), wi) != Boolean::No);
 922
 923     printStringNoNewline(inp, "The number of independent AWH biases");
 924     opt         = "awh-nbias";
 925     int numBias = get_eint(inp, opt, 1, wi);
 926     /* Check this before starting to read the AWH biases. */
 927     if (numBias <= 0)
 928     {
 929         gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
 930     }
 931
 932     /* Read the parameters specific to each AWH bias */
 933     for (int k = 0; k < numBias; k++)
 934     {
 935         bool        bComment  = (k == 0);
 936         std::string prefixawh = formatString("awh%d", k + 1);
 937         awhBiasParams_.emplace_back(inp, prefixawh, ir, wi, bComment);
 938     }
 939     checkAwhParams(*this, ir, wi);
 940     checkInputConsistencyAwh(*this, wi);
 941 }
 942
 943 AwhParams::AwhParams(ISerializer* serializer)
 944 {
 945     GMX_RELEASE_ASSERT(serializer->reading(),
 946                        "Can not use writing serializer to read AWH parameters");
 947     int numberOfBiases = awhBiasParams_.size();
 948     serializer->doInt(&numberOfBiases);
 949     serializer->doInt(&nstOut_);
 950     serializer->doInt64(&seed_);
 951     serializer->doInt(&nstSampleCoord_);
 952     serializer->doInt(&numSamplesUpdateFreeEnergy_);
 953     serializer->doEnumAsInt(&potentialEnum_);
 954     serializer->doBool(&shareBiasMultisim_);
 955
 956     if (numberOfBiases > 0)
 957     {
 958         for (int k = 0; k < numberOfBiases; k++)
 959         {
 960             awhBiasParams_.emplace_back(serializer);
 961         }
 962     }
 963     checkInputConsistencyAwh(*this, nullptr);
 964 }
 965
 966 void AwhParams::serialize(ISerializer* serializer)
 967 {
 968     GMX_RELEASE_ASSERT(!serializer->reading(),
 969                        "Can not use reading serializer to write AWH parameters");
 970     int numberOfBiases = numBias();
 971     serializer->doInt(&numberOfBiases);
 972     serializer->doInt(&nstOut_);
 973     serializer->doInt64(&seed_);
 974     serializer->doInt(&nstSampleCoord_);
 975     serializer->doInt(&numSamplesUpdateFreeEnergy_);
 976     serializer->doEnumAsInt(&potentialEnum_);
 977     serializer->doBool(&shareBiasMultisim_);
 978
 979     if (numberOfBiases > 0)
 980     {
 981         for (int k = 0; k < numberOfBiases; k++)
 982         {
 983             awhBiasParams_[k].serialize(serializer);
 984         }
 985     }
 986 }
 987
 988 /*! \brief
 989  * Gets the period of a pull coordinate.
 990  *
 991  * \param[in] pullCoordParams   The parameters for the pull coordinate.
 992  * \param[in] pbc               The PBC setup
 993  * \param[in] intervalLength    The length of the AWH interval for this pull coordinate
 994  * \returns the period (or 0 if not periodic).
 995  */
 996 static double get_pull_coord_period(const t_pull_coord& pullCoordParams, const t_pbc& pbc, const real intervalLength)
 997 {
 998     double period = 0;
 999
1000     if (pullCoordParams.eGeom == PullGroupGeometry::Direction)
1001     {
1002         const real margin = 0.001;
1003         // Make dims periodic when the interval covers > 95%
1004         const real periodicFraction = 0.95;
1005
1006         // Check if the pull direction is along a box vector
1007         for (int dim = 0; dim < pbc.ndim_ePBC; dim++)
1008         {
1009             const real boxLength    = norm(pbc.box[dim]);
1010             const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]);
1011             if (innerProduct >= (1 - margin) * boxLength && innerProduct <= (1 + margin) * boxLength)
1012             {
1013                 if (intervalLength > (1 + margin) * boxLength)
1014                 {
1015                     gmx_fatal(FARGS,
1016                               "The AWH interval (%f nm) for a pull coordinate is larger than the "
1017                               "box size (%f nm)",
1018                               intervalLength,
1019                               boxLength);
1020                 }
1021
1022                 if (intervalLength > periodicFraction * boxLength)
1023                 {
1024                     period = boxLength;
1025                 }
1026             }
1027         }
1028     }
1029     else if (pullCoordParams.eGeom == PullGroupGeometry::Dihedral)
1030     {
1031         /* The dihedral angle is periodic in -180 to 180 deg */
1032         period = 360;
1033     }
1034
1035     return period;
1036 }
1037
1038 /*! \brief
1039  * Checks if the given interval is defined in the correct periodic interval.
1040  *
1041  * \param[in] origin      Start value of interval.
1042  * \param[in] end         End value of interval.
1043  * \param[in] period      Period (or 0 if not periodic).
1044  * \returns true if the end point values are in the correct periodic interval.
1045  */
1046 static bool intervalIsInPeriodicInterval(double origin, double end, double period)
1047 {
1048     return (period == 0) || (std::fabs(origin) <= 0.5 * period && std::fabs(end) <= 0.5 * period);
1049 }
1050
1051 /*! \brief
1052  * Checks if a value is within an interval.
1053  *
1054  * \param[in] origin      Start value of interval.
1055  * \param[in] end         End value of interval.
1056  * \param[in] period      Period (or 0 if not periodic).
1057  * \param[in] value       Value to check.
1058  * \returns true if the value is within the interval.
1059  */
1060 static bool valueIsInInterval(double origin, double end, double period, double value)
1061 {
1062     bool bIn_interval;
1063
1064     if (period > 0)
1065     {
1066         if (origin < end)
1067         {
1068             /* The interval closes within the periodic interval */
1069             bIn_interval = (value >= origin) && (value <= end);
1070         }
1071         else
1072         {
1073             /* The interval wraps around the periodic boundary */
1074             bIn_interval = ((value >= origin) && (value <= 0.5 * period))
1075                            || ((value >= -0.5 * period) && (value <= end));
1076         }
1077     }
1078     else
1079     {
1080         bIn_interval = (value >= origin) && (value <= end);
1081     }
1082
1083     return bIn_interval;
1084 }
1085
1086 /*! \brief
1087  * Check if the starting configuration is consistent with the given interval.
1088  *
1089  * \param[in]     awhParams  AWH parameters.
1090  * \param[in,out] wi         Struct for bookeeping warnings.
1091  */
1092 static void checkInputConsistencyInterval(const AwhParams& awhParams, warninp_t wi)
1093 {
1094     const auto& awhBiasParams = awhParams.awhBiasParams();
1095     for (int k = 0; k < gmx::ssize(awhBiasParams); k++)
1096     {
1097         const auto& dimParams = awhBiasParams[k].dimParams();
1098         for (int d = 0; d < gmx::ssize(dimParams); d++)
1099         {
1100             int    coordIndex = dimParams[d].coordinateIndex();
1101             double origin = dimParams[d].origin(), end = dimParams[d].end(),
1102                    period         = dimParams[d].period();
1103             double coordValueInit = dimParams[d].initialCoordinate();
1104
1105             if ((period == 0) && (origin > end))
1106             {
1107                 gmx_fatal(FARGS,
1108                           "For the non-periodic pull coordinates awh%d-dim%d-start (%f) cannot be "
1109                           "larger than awh%d-dim%d-end (%f)",
1110                           k + 1,
1111                           d + 1,
1112                           origin,
1113                           k + 1,
1114                           d + 1,
1115                           end);
1116             }
1117
1118             /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
1119                Make sure the AWH interval is within this reference interval.
1120
1121                Note: we could fairly simply allow using a  more general interval (e.g. [x, x + period]) but it complicates
1122                things slightly and I don't see that there is a great need for it. It would also mean that the interval would
1123                depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180,
1124                independent of AWH parameters.
1125              */
1126             if (!intervalIsInPeriodicInterval(origin, end, period))
1127             {
1128                 gmx_fatal(FARGS,
1129                           "When using AWH with periodic pull coordinate geometries "
1130                           "awh%d-dim%d-start (%.8g) and "
1131                           "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take "
1132                           "values in between "
1133                           "minus half a period and plus half a period, i.e. in the interval [%.8g, "
1134                           "%.8g].",
1135                           k + 1,
1136                           d + 1,
1137                           origin,
1138                           k + 1,
1139                           d + 1,
1140                           end,
1141                           period,
1142                           -0.5 * period,
1143                           0.5 * period);
1144             }
1145
1146             /* Warn if the pull initial coordinate value is not in the grid */
1147             if (!valueIsInInterval(origin, end, period, coordValueInit))
1148             {
1149                 auto message = formatString(
1150                         "The initial coordinate value (%.8g) for pull coordinate index %d falls "
1151                         "outside "
1152                         "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end "
1153                         "(%.8g). "
1154                         "This can lead to large initial forces pulling the coordinate towards the "
1155                         "sampling interval.",
1156                         coordValueInit,
1157                         coordIndex + 1,
1158                         k + 1,
1159                         d + 1,
1160                         origin,
1161                         k + 1,
1162                         d + 1,
1163                         end);
1164                 warning(wi, message);
1165             }
1166         }
1167     }
1168 }
1169
1170 /*! \brief
1171  * Sets AWH parameters, for one AWH pull dimension.
1172  *
1173  * \param[in,out] dimParams          AWH dimension parameters.
1174  * \param[in]     biasIndex             The index of the bias containing this AWH pull dimension.
1175  * \param[in]     dimIndex              The index of this AWH pull dimension.
1176  * \param[in]     pull_params           Pull parameters.
1177  * \param[in,out] pull_work             Pull working struct to register AWH bias in.
1178  * \param[in]     pbc                   A pbc information structure.
1179  * \param[in]     compressibility       Compressibility matrix for pressure coupling,
1180  * pass all 0 without pressure coupling.
1181  * \param[in,out] wi                    Struct for bookeeping warnings.
1182  *
1183  */
1184 static void setStateDependentAwhPullDimParams(AwhDimParams*        dimParams,
1185                                               const int            biasIndex,
1186                                               const int            dimIndex,
1187                                               const pull_params_t& pull_params,
1188                                               pull_t*              pull_work,
1189                                               const t_pbc&         pbc,
1190                                               const tensor&        compressibility,
1191                                               warninp_t            wi)
1192 {
1193     const t_pull_coord& pullCoordParams = pull_params.coord[dimParams->coordinateIndex()];
1194
1195     if (pullCoordParams.eGeom == PullGroupGeometry::DirectionPBC)
1196     {
1197         gmx_fatal(FARGS,
1198                   "AWH does not support pull geometry '%s'. "
1199                   "If the maximum distance between the groups is always "
1200                   "less than half the box size, "
1201                   "you can use geometry '%s' instead.",
1202                   enumValueToString(PullGroupGeometry::DirectionPBC),
1203                   enumValueToString(PullGroupGeometry::Direction));
1204     }
1205
1206     dimParams->setPeriod(
1207             get_pull_coord_period(pullCoordParams, pbc, dimParams->end() - dimParams->origin()));
1208     // We would like to check for scaling, but we don't have the full inputrec available here
1209     if (dimParams->period() > 0
1210         && !(pullCoordParams.eGeom == PullGroupGeometry::Angle
1211              || pullCoordParams.eGeom == PullGroupGeometry::Dihedral))
1212     {
1213         bool coordIsScaled = false;
1214         for (int d2 = 0; d2 < DIM; d2++)
1215         {
1216             if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0)
1217             {
1218                 coordIsScaled = true;
1219             }
1220         }
1221         if (coordIsScaled)
1222         {
1223             std::string mesg = gmx::formatString(
1224                     "AWH dimension %d of bias %d is periodic with pull geometry '%s', "
1225                     "while you should be applying pressure scaling to the "
1226                     "corresponding box vector, this is not supported.",
1227                     dimIndex + 1,
1228                     biasIndex + 1,
1229                     enumValueToString(pullCoordParams.eGeom));
1230             warning(wi, mesg.c_str());
1231         }
1232     }
1233
1234     /* The initial coordinate value, converted to external user units. */
1235     double initialCoordinate = get_pull_coord_value(pull_work, dimParams->coordinateIndex(), &pbc);
1236     initialCoordinate *= pull_conversion_factor_internal2userinput(&pullCoordParams);
1237     dimParams->setInitialCoordinate(initialCoordinate);
1238 }
1239
1240 void setStateDependentAwhParams(AwhParams*           awhParams,
1241                                 const pull_params_t& pull_params,
1242                                 pull_t*              pull_work,
1243                                 const matrix         box,
1244                                 PbcType              pbcType,
1245                                 const tensor&        compressibility,
1246                                 const t_grpopts*     inputrecGroupOptions,
1247                                 const real           initLambda,
1248                                 const gmx_mtop_t&    mtop,
1249                                 warninp_t            wi)
1250 {
1251     /* The temperature is not really state depenendent but is not known
1252      * when read_awhParams is called (in get ir).
1253      * It is known first after do_index has been called in grompp.cpp.
1254      */
1255     if (inputrecGroupOptions->ref_t == nullptr || inputrecGroupOptions->ref_t[0] <= 0)
1256     {
1257         gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
1258     }
1259     for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
1260     {
1261         if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
1262         {
1263             gmx_fatal(FARGS,
1264                       "AWH biasing is currently only supported for identical temperatures for all "
1265                       "temperature coupling groups");
1266         }
1267     }
1268
1269     t_pbc pbc;
1270     set_pbc(&pbc, pbcType, box);
1271
1272     auto awhBiasParams = awhParams->awhBiasParams();
1273     for (int k = 0; k < awhParams->numBias(); k++)
1274     {
1275         auto awhBiasDimensionParams = awhBiasParams[k].dimParams();
1276         for (int d = 0; d < awhBiasParams[k].ndim(); d++)
1277         {
1278             AwhDimParams* dimParams = &awhBiasDimensionParams[d];
1279             if (dimParams->coordinateProvider() == AwhCoordinateProviderType::Pull)
1280             {
1281                 setStateDependentAwhPullDimParams(
1282                         dimParams, k, d, pull_params, pull_work, pbc, compressibility, wi);
1283             }
1284             else
1285             {
1286                 dimParams->setInitialCoordinate(initLambda);
1287                 checkFepLambdaDimDecouplingConsistency(mtop, wi);
1288             }
1289         }
1290     }
1291     checkInputConsistencyInterval(*awhParams, wi);
1292
1293     /* Register AWH as external potential with pull (for AWH dimensions that use pull as
1294      * reaction coordinate provider) to check consistency. */
1295     Awh::registerAwhWithPull(*awhParams, pull_work);
1296 }
1297
1298 void checkAwhParams(const AwhParams& awhParams, const t_inputrec& ir, warninp_t wi)
1299 {
1300     std::string opt;
1301     checkMtsConsistency(ir, wi);
1302
1303     opt = "awh-nstout";
1304     if (awhParams.nstout() <= 0)
1305     {
1306         auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
1307                                     opt.c_str(),
1308                                     awhParams.nstout());
1309         warning_error(wi, message);
1310     }
1311     /* This restriction can be removed by changing a flag of print_ebin() */
1312     if (ir.nstenergy == 0 || awhParams.nstout() % ir.nstenergy != 0)
1313     {
1314         auto message = formatString(
1315                 "%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(), awhParams.nstout(), ir.nstenergy);
1316         warning_error(wi, message);
1317     }
1318
1319     opt = "awh-nsamples-update";
1320     if (awhParams.numSamplesUpdateFreeEnergy() <= 0)
1321     {
1322         warning_error(wi, opt + " needs to be an integer > 0");
1323     }
1324
1325     for (int k = 0; k < awhParams.numBias(); k++)
1326     {
1327         std::string prefixawh = formatString("awh%d", k + 1);
1328         checkBiasParams(awhParams.awhBiasParams()[k], prefixawh, ir, wi);
1329     }
1330
1331     if (ir.init_step != 0)
1332     {
1333         warning_error(wi, "With AWH init-step should be 0");
1334     }
1335 }
1336
1337 } // namespace gmx