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