src/gromacs/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, 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/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/pull_params.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/pulling/pull.h"
  54 #include "gromacs/random/seed.h"
  55 #include "gromacs/utility/cstringutil.h"
  56 #include "gromacs/utility/fatalerror.h"
  57 #include "gromacs/utility/smalloc.h"
  58 #include "gromacs/utility/stringutil.h"
  59
  60 #include "biasparams.h"
  61 #include "biassharing.h"
  62
  63 namespace gmx
  64 {
  65
  66 const char* eawhtarget_names[eawhtargetNR + 1] = { "constant", "cutoff", "boltzmann",
  67                                                    "local-boltzmann", nullptr };
  68
  69 const char* eawhgrowth_names[eawhgrowthNR + 1] = { "exp-linear", "linear", nullptr };
  70
  71 const char* eawhpotential_names[eawhpotentialNR + 1] = { "convolved", "umbrella", nullptr };
  72
  73 const char* eawhcoordprovider_names[eawhcoordproviderNR + 1] = { "pull", nullptr };
  74
  75 /*! \brief
  76  * Read parameters of an AWH bias dimension.
  77  *
  78  * \param[in,out] inp        Input file entries.
  79  * \param[in] prefix         Prefix for dimension parameters.
  80  * \param[in,out] dimParams  AWH dimensional parameters.
  81  * \param[in] pull_params    Pull parameters.
  82  * \param[in,out] wi         Struct for bookeeping warnings.
  83  * \param[in] bComment       True if comments should be printed.
  84  */
  85 static void readDimParams(std::vector<t_inpfile>* inp,
  86                           const std::string&      prefix,
  87                           AwhDimParams*           dimParams,
  88                           const pull_params_t*    pull_params,
  89                           warninp_t               wi,
  90                           bool                    bComment)
  91 {
  92     std::string opt;
  93     if (bComment)
  94     {
  95         printStringNoNewline(
  96                 inp, "The provider of the reaction coordinate, currently only pull is supported");
  97     }
  98
  99     opt                       = prefix + "-coord-provider";
 100     dimParams->eCoordProvider = get_eeenum(inp, opt, eawhcoordprovider_names, wi);
 101
 102     if (bComment)
 103     {
 104         printStringNoNewline(inp, "The coordinate index for this dimension");
 105     }
 106     opt = prefix + "-coord-index";
 107     int coordIndexInput;
 108     coordIndexInput = get_eint(inp, opt, 1, wi);
 109     if (coordIndexInput < 1)
 110     {
 111         gmx_fatal(FARGS,
 112                   "Failed to read a valid coordinate index for %s. "
 113                   "Note that the pull coordinate indexing starts at 1.",
 114                   opt.c_str());
 115     }
 116
 117     /* The pull coordinate indices start at 1 in the input file, at 0 internally */
 118     dimParams->coordIndex = coordIndexInput - 1;
 119
 120     /* The pull settings need to be consistent with the AWH settings */
 121     if (!(pull_params->coord[dimParams->coordIndex].eType == epullEXTERNAL))
 122     {
 123         gmx_fatal(FARGS, "AWH biasing can only be  applied to pull type %s", EPULLTYPE(epullEXTERNAL));
 124     }
 125
 126     if (dimParams->coordIndex >= pull_params->ncoord)
 127     {
 128         gmx_fatal(FARGS,
 129                   "The given AWH coordinate index (%d) is larger than the number of pull "
 130                   "coordinates (%d)",
 131                   coordIndexInput, pull_params->ncoord);
 132     }
 133     if (pull_params->coord[dimParams->coordIndex].rate != 0)
 134     {
 135         auto message = formatString(
 136                 "Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate",
 137                 coordIndexInput, pull_params->coord[dimParams->coordIndex].rate);
 138         warning_error(wi, message);
 139     }
 140
 141     /* Grid params for each axis */
 142     int eGeom = pull_params->coord[dimParams->coordIndex].eGeom;
 143
 144     if (bComment)
 145     {
 146         printStringNoNewline(inp, "Start and end values for each coordinate dimension");
 147     }
 148
 149     opt               = prefix + "-start";
 150     dimParams->origin = get_ereal(inp, opt, 0., wi);
 151
 152     opt            = prefix + "-end";
 153     dimParams->end = get_ereal(inp, opt, 0., wi);
 154
 155     if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
 156     {
 157         auto message = formatString(
 158                 "The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
 159                 "This will result in only one point along this axis in the coordinate value grid.",
 160                 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end);
 161         warning(wi, message);
 162     }
 163     /* Check that the requested interval is in allowed range */
 164     if (eGeom == epullgDIST)
 165     {
 166         if (dimParams->origin < 0 || dimParams->end < 0)
 167         {
 168             gmx_fatal(FARGS,
 169                       "%s-start (%g) or %s-end (%g) set to a negative value. With pull geometry "
 170                       "distance coordinate values are non-negative. "
 171                       "Perhaps you want to use geometry %s instead?",
 172                       prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
 173                       EPULLGEOM(epullgDIR));
 174         }
 175     }
 176     else if (eGeom == epullgANGLE || eGeom == epullgANGLEAXIS)
 177     {
 178         if (dimParams->origin < 0 || dimParams->end > 180)
 179         {
 180             gmx_fatal(FARGS,
 181                       "%s-start (%g) and %s-end (%g) are outside of the allowed range 0 to 180 deg "
 182                       "for pull geometries %s and %s ",
 183                       prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
 184                       EPULLGEOM(epullgANGLE), EPULLGEOM(epullgANGLEAXIS));
 185         }
 186     }
 187     else if (eGeom == epullgDIHEDRAL)
 188     {
 189         if (dimParams->origin < -180 || dimParams->end > 180)
 190         {
 191             gmx_fatal(FARGS,
 192                       "%s-start (%g) and %s-end (%g) are outside of the allowed range -180 to 180 "
 193                       "deg for pull geometry %s. ",
 194                       prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
 195                       EPULLGEOM(epullgDIHEDRAL));
 196         }
 197     }
 198
 199     if (bComment)
 200     {
 201         printStringNoNewline(
 202                 inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
 203     }
 204     opt                      = prefix + "-force-constant";
 205     dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
 206     if (dimParams->forceConstant <= 0)
 207     {
 208         warning_error(wi, "The force AWH bias force constant should be > 0");
 209     }
 210
 211     if (bComment)
 212     {
 213         printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps)");
 214     }
 215     opt                  = prefix + "-diffusion";
 216     dimParams->diffusion = get_ereal(inp, opt, 0, wi);
 217
 218     if (dimParams->diffusion <= 0)
 219     {
 220         const double diffusion_default = 1e-5;
 221         auto         message           = formatString(
 222                 "%s not explicitly set by user. You can choose to use a default "
 223                 "value (%g nm^2/ps or rad^2/ps) but this may very well be "
 224                 "non-optimal for your system!",
 225                 opt.c_str(), diffusion_default);
 226         warning(wi, message);
 227         dimParams->diffusion = diffusion_default;
 228     }
 229
 230     if (bComment)
 231     {
 232         printStringNoNewline(inp,
 233                              "Diameter that needs to be sampled around a point before it is "
 234                              "considered covered.");
 235     }
 236     opt                      = prefix + "-cover-diameter";
 237     dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
 238
 239     if (dimParams->coverDiameter < 0)
 240     {
 241         gmx_fatal(FARGS, "%s (%g) cannot be negative.", opt.c_str(), dimParams->coverDiameter);
 242     }
 243 }
 244
 245 /*! \brief
 246  * Check consistency of input at the AWH bias level.
 247  *
 248  * \param[in]     awhBiasParams  AWH bias parameters.
 249  * \param[in,out] wi             Struct for bookkeeping warnings.
 250  */
 251 static void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
 252 {
 253     /* Covering diameter and sharing warning. */
 254     for (int d = 0; d < awhBiasParams.ndim; d++)
 255     {
 256         double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
 257         if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
 258         {
 259             warning(wi,
 260                     "The covering diameter is only relevant to set for bias sharing simulations.");
 261         }
 262     }
 263 }
 264
 265 /*! \brief
 266  * Read parameters of an AWH bias.
 267  *
 268  * \param[in,out] inp            Input file entries.
 269  * \param[in,out] awhBiasParams  AWH dimensional parameters.
 270  * \param[in]     prefix         Prefix for bias parameters.
 271  * \param[in]     ir             Input parameter struct.
 272  * \param[in,out] wi             Struct for bookeeping warnings.
 273  * \param[in]     bComment       True if comments should be printed.
 274  */
 275 static void read_bias_params(std::vector<t_inpfile>* inp,
 276                              AwhBiasParams*          awhBiasParams,
 277                              const std::string&      prefix,
 278                              const t_inputrec*       ir,
 279                              warninp_t               wi,
 280                              bool                    bComment)
 281 {
 282     if (bComment)
 283     {
 284         printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
 285     }
 286
 287     std::string opt = prefix + "-error-init";
 288     /* We allow using a default value here without warning (but warn the user if the diffusion constant is not set). */
 289     awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi);
 290     if (awhBiasParams->errorInitial <= 0)
 291     {
 292         gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
 293     }
 294
 295     if (bComment)
 296     {
 297         printStringNoNewline(inp,
 298                              "Growth rate of the reference histogram determining the bias update "
 299                              "size: exp-linear or linear");
 300     }
 301     opt                    = prefix + "-growth";
 302     awhBiasParams->eGrowth = get_eeenum(inp, opt, eawhgrowth_names, wi);
 303
 304     if (bComment)
 305     {
 306         printStringNoNewline(inp,
 307                              "Start the simulation by equilibrating histogram towards the target "
 308                              "distribution: no or yes");
 309     }
 310     opt                                 = prefix + "-equilibrate-histogram";
 311     awhBiasParams->equilibrateHistogram = (get_eeenum(inp, opt, yesno_names, wi) != 0);
 312     if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != eawhgrowthEXP_LINEAR)
 313     {
 314         auto message =
 315                 formatString("Option %s will only have an effect for histogram growth type '%s'.",
 316                              opt.c_str(), EAWHGROWTH(eawhgrowthEXP_LINEAR));
 317         warning(wi, message);
 318     }
 319
 320     if (bComment)
 321     {
 322         printStringNoNewline(
 323                 inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
 324     }
 325     opt                    = prefix + "-target";
 326     awhBiasParams->eTarget = get_eeenum(inp, opt, eawhtarget_names, wi);
 327
 328     if ((awhBiasParams->eTarget == eawhtargetLOCALBOLTZMANN)
 329         && (awhBiasParams->eGrowth == eawhgrowthEXP_LINEAR))
 330     {
 331         auto message = formatString(
 332                 "Target type '%s' combined with histogram growth type '%s' is not "
 333                 "expected to give stable bias updates. You probably want to use growth type "
 334                 "'%s' instead.",
 335                 EAWHTARGET(eawhtargetLOCALBOLTZMANN), EAWHGROWTH(eawhgrowthEXP_LINEAR),
 336                 EAWHGROWTH(eawhgrowthLINEAR));
 337         warning(wi, message);
 338     }
 339
 340     if (bComment)
 341     {
 342         printStringNoNewline(inp,
 343                              "Boltzmann beta scaling factor for target distribution types "
 344                              "'boltzmann' and 'boltzmann-local'");
 345     }
 346     opt                              = prefix + "-target-beta-scaling";
 347     awhBiasParams->targetBetaScaling = get_ereal(inp, opt, 0, wi);
 348
 349     switch (awhBiasParams->eTarget)
 350     {
 351         case eawhtargetBOLTZMANN:
 352         case eawhtargetLOCALBOLTZMANN:
 353             if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
 354             {
 355                 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
 356                           awhBiasParams->targetBetaScaling, EAWHTARGET(awhBiasParams->eTarget));
 357             }
 358             break;
 359         default:
 360             if (awhBiasParams->targetBetaScaling != 0)
 361             {
 362                 gmx_fatal(
 363                         FARGS,
 364                         "Value for %s (%g) set explicitly but will not be used for target type %s.",
 365                         opt.c_str(), awhBiasParams->targetBetaScaling,
 366                         EAWHTARGET(awhBiasParams->eTarget));
 367             }
 368             break;
 369     }
 370
 371     if (bComment)
 372     {
 373         printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
 374     }
 375     opt                         = prefix + "-target-cutoff";
 376     awhBiasParams->targetCutoff = get_ereal(inp, opt, 0, wi);
 377
 378     switch (awhBiasParams->eTarget)
 379     {
 380         case eawhtargetCUTOFF:
 381             if (awhBiasParams->targetCutoff <= 0)
 382             {
 383                 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
 384                           awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
 385             }
 386             break;
 387         default:
 388             if (awhBiasParams->targetCutoff != 0)
 389             {
 390                 gmx_fatal(
 391                         FARGS,
 392                         "Value for %s (%g) set explicitly but will not be used for target type %s.",
 393                         opt.c_str(), awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
 394             }
 395             break;
 396     }
 397
 398     if (bComment)
 399     {
 400         printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
 401     }
 402     opt                      = prefix + "-user-data";
 403     awhBiasParams->bUserData = get_eeenum(inp, opt, yesno_names, wi);
 404
 405     if (bComment)
 406     {
 407         printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
 408     }
 409     opt                       = prefix + "-share-group";
 410     awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
 411     if (awhBiasParams->shareGroup < 0)
 412     {
 413         warning_error(wi, "AWH bias share-group should be >= 0");
 414     }
 415
 416     if (bComment)
 417     {
 418         printStringNoNewline(inp, "Dimensionality of the coordinate");
 419     }
 420     opt                 = prefix + "-ndim";
 421     awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
 422
 423     if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
 424     {
 425         gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(),
 426                   awhBiasParams->ndim, c_biasMaxNumDim);
 427     }
 428     if (awhBiasParams->ndim > 2)
 429     {
 430         warning_note(wi,
 431                      "For awh-dim > 2 the estimate based on the diffusion and the initial error is "
 432                      "currently only a rough guideline."
 433                      " You should verify its usefulness for your system before production runs!");
 434     }
 435     snew(awhBiasParams->dimParams, awhBiasParams->ndim);
 436     for (int d = 0; d < awhBiasParams->ndim; d++)
 437     {
 438         bComment              = bComment && d == 0;
 439         std::string prefixdim = prefix + formatString("-dim%d", d + 1);
 440         readDimParams(inp, prefixdim, &awhBiasParams->dimParams[d], ir->pull, wi, bComment);
 441     }
 442
 443     /* Check consistencies here that cannot be checked at read time at a lower level. */
 444     checkInputConsistencyAwhBias(*awhBiasParams, wi);
 445 }
 446
 447 /*! \brief
 448  * Check consistency of input at the AWH level.
 449  *
 450  * \param[in]     awhParams  AWH parameters.
 451  * \param[in,out] wi         Struct for bookkeeping warnings.
 452  */
 453 static void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi)
 454 {
 455     /* Each pull coord can map to at most 1 AWH coord.
 456      * Check that we have a shared bias when requesting multisim sharing.
 457      */
 458     bool haveSharedBias = false;
 459     for (int k1 = 0; k1 < awhParams.numBias; k1++)
 460     {
 461         const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1];
 462
 463         if (awhBiasParams1.shareGroup > 0)
 464         {
 465             haveSharedBias = true;
 466         }
 467
 468         /* k1 is the reference AWH, k2 is the AWH we compare with (can be equal to k1) */
 469         for (int k2 = k1; k2 < awhParams.numBias; k2++)
 470         {
 471             for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
 472             {
 473                 const AwhBiasParams& awhBiasParams2 = awhParams.awhBiasParams[k2];
 474
 475                 /* d1 is the reference dimension of the reference AWH. d2 is the dim index of the AWH to compare with. */
 476                 for (int d2 = 0; d2 < awhBiasParams2.ndim; d2++)
 477                 {
 478                     /* Give an error if (d1, k1) is different from (d2, k2) but the pull coordinate is the same */
 479                     if ((d1 != d2 || k1 != k2)
 480                         && (awhBiasParams1.dimParams[d1].coordIndex == awhBiasParams2.dimParams[d2].coordIndex))
 481                     {
 482                         char errormsg[STRLEN];
 483                         sprintf(errormsg,
 484                                 "One pull coordinate (%d) cannot be mapped to two separate AWH "
 485                                 "dimensions (awh%d-dim%d and awh%d-dim%d). "
 486                                 "If this is really what you want to do you will have to duplicate "
 487                                 "this pull coordinate.",
 488                                 awhBiasParams1.dimParams[d1].coordIndex + 1, k1 + 1, d1 + 1, k2 + 1,
 489                                 d2 + 1);
 490                         gmx_fatal(FARGS, "%s", errormsg);
 491                     }
 492                 }
 493             }
 494         }
 495     }
 496
 497     if (awhParams.shareBiasMultisim && !haveSharedBias)
 498     {
 499         warning(wi,
 500                 "Sharing of biases over multiple simulations is requested, but no bias is marked "
 501                 "as shared (share-group > 0)");
 502     }
 503
 504     /* mdrun does not support this (yet), but will check again */
 505     if (haveBiasSharingWithinSimulation(awhParams))
 506     {
 507         warning(wi,
 508                 "You have shared biases within a single simulation, but mdrun does not support "
 509                 "this (yet)");
 510     }
 511 }
 512
 513 AwhParams* readAndCheckAwhParams(std::vector<t_inpfile>* inp, const t_inputrec* ir, warninp_t wi)
 514 {
 515     AwhParams* awhParams;
 516     snew(awhParams, 1);
 517     std::string opt;
 518
 519     /* Parameters common for all biases */
 520
 521     printStringNoNewline(inp, "The way to apply the biasing potential: convolved or umbrella");
 522     opt                   = "awh-potential";
 523     awhParams->ePotential = get_eeenum(inp, opt, eawhpotential_names, wi);
 524
 525     printStringNoNewline(inp,
 526                          "The random seed used for sampling the umbrella center in the case of "
 527                          "umbrella type potential");
 528     opt             = "awh-seed";
 529     awhParams->seed = get_eint(inp, opt, -1, wi);
 530     if (awhParams->seed == -1)
 531     {
 532         awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
 533         fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
 534     }
 535
 536     printStringNoNewline(inp, "Data output interval in number of steps");
 537     opt               = "awh-nstout";
 538     awhParams->nstOut = get_eint(inp, opt, 100000, wi);
 539     if (awhParams->nstOut <= 0)
 540     {
 541         auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
 542                                     opt.c_str(), awhParams->nstOut);
 543         warning_error(wi, message);
 544     }
 545     /* This restriction can be removed by changing a flag of print_ebin() */
 546     if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
 547     {
 548         auto message = formatString("%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(),
 549                                     awhParams->nstOut, ir->nstenergy);
 550         warning_error(wi, message);
 551     }
 552
 553     printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
 554     opt                       = "awh-nstsample";
 555     awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
 556
 557     printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
 558     opt                                   = "awh-nsamples-update";
 559     awhParams->numSamplesUpdateFreeEnergy = get_eint(inp, opt, 10, wi);
 560     if (awhParams->numSamplesUpdateFreeEnergy <= 0)
 561     {
 562         warning_error(wi, opt + " needs to be an integer > 0");
 563     }
 564
 565     printStringNoNewline(
 566             inp, "When true, biases with share-group>0 are shared between multiple simulations");
 567     opt                          = "awh-share-multisim";
 568     awhParams->shareBiasMultisim = (get_eeenum(inp, opt, yesno_names, wi) != 0);
 569
 570     printStringNoNewline(inp, "The number of independent AWH biases");
 571     opt                = "awh-nbias";
 572     awhParams->numBias = get_eint(inp, opt, 1, wi);
 573     if (awhParams->numBias <= 0)
 574     {
 575         gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
 576     }
 577
 578     /* Read the parameters specific to each AWH bias */
 579     snew(awhParams->awhBiasParams, awhParams->numBias);
 580
 581     for (int k = 0; k < awhParams->numBias; k++)
 582     {
 583         bool        bComment  = (k == 0);
 584         std::string prefixawh = formatString("awh%d", k + 1);
 585         read_bias_params(inp, &awhParams->awhBiasParams[k], prefixawh, ir, wi, bComment);
 586     }
 587
 588     /* Do a final consistency check before returning */
 589     checkInputConsistencyAwh(*awhParams, wi);
 590
 591     if (ir->init_step != 0)
 592     {
 593         warning_error(wi, "With AWH init-step should be 0");
 594     }
 595
 596     return awhParams;
 597 }
 598
 599 /*! \brief
 600  * Gets the period of a pull coordinate.
 601  *
 602  * \param[in] pullCoordParams   The parameters for the pull coordinate.
 603  * \param[in] pbc               The PBC setup
 604  * \param[in] intervalLength    The length of the AWH interval for this pull coordinate
 605  * \returns the period (or 0 if not periodic).
 606  */
 607 static double get_pull_coord_period(const t_pull_coord& pullCoordParams, const t_pbc& pbc, const real intervalLength)
 608 {
 609     double period = 0;
 610
 611     if (pullCoordParams.eGeom == epullgDIR)
 612     {
 613         const real margin = 0.001;
 614         // Make dims periodic when the interval covers > 95%
 615         const real periodicFraction = 0.95;
 616
 617         // Check if the pull direction is along a box vector
 618         for (int dim = 0; dim < pbc.ndim_ePBC; dim++)
 619         {
 620             const real boxLength    = norm(pbc.box[dim]);
 621             const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]);
 622             if (innerProduct >= (1 - margin) * boxLength && innerProduct <= (1 + margin) * boxLength)
 623             {
 624                 if (intervalLength > (1 + margin) * boxLength)
 625                 {
 626                     gmx_fatal(FARGS,
 627                               "The AWH interval (%f nm) for a pull coordinate is larger than the "
 628                               "box size (%f nm)",
 629                               intervalLength, boxLength);
 630                 }
 631
 632                 if (intervalLength > periodicFraction * boxLength)
 633                 {
 634                     period = boxLength;
 635                 }
 636             }
 637         }
 638     }
 639     else if (pullCoordParams.eGeom == epullgDIHEDRAL)
 640     {
 641         /* The dihedral angle is periodic in -180 to 180 deg */
 642         period = 360;
 643     }
 644
 645     return period;
 646 }
 647
 648 /*! \brief
 649  * Checks if the given interval is defined in the correct periodic interval.
 650  *
 651  * \param[in] origin      Start value of interval.
 652  * \param[in] end         End value of interval.
 653  * \param[in] period      Period (or 0 if not periodic).
 654  * \returns true if the end point values are in the correct periodic interval.
 655  */
 656 static bool intervalIsInPeriodicInterval(double origin, double end, double period)
 657 {
 658     return (period == 0) || (std::fabs(origin) <= 0.5 * period && std::fabs(end) <= 0.5 * period);
 659 }
 660
 661 /*! \brief
 662  * Checks if a value is within an interval.
 663  *
 664  * \param[in] origin      Start value of interval.
 665  * \param[in] end         End value of interval.
 666  * \param[in] period      Period (or 0 if not periodic).
 667  * \param[in] value       Value to check.
 668  * \returns true if the value is within the interval.
 669  */
 670 static bool valueIsInInterval(double origin, double end, double period, double value)
 671 {
 672     bool bIn_interval;
 673
 674     if (period > 0)
 675     {
 676         if (origin < end)
 677         {
 678             /* The interval closes within the periodic interval */
 679             bIn_interval = (value >= origin) && (value <= end);
 680         }
 681         else
 682         {
 683             /* The interval wraps around the periodic boundary */
 684             bIn_interval = ((value >= origin) && (value <= 0.5 * period))
 685                            || ((value >= -0.5 * period) && (value <= end));
 686         }
 687     }
 688     else
 689     {
 690         bIn_interval = (value >= origin) && (value <= end);
 691     }
 692
 693     return bIn_interval;
 694 }
 695
 696 /*! \brief
 697  * Check if the starting configuration is consistent with the given interval.
 698  *
 699  * \param[in]     awhParams  AWH parameters.
 700  * \param[in,out] wi         Struct for bookeeping warnings.
 701  */
 702 static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t wi)
 703 {
 704     for (int k = 0; k < awhParams->numBias; k++)
 705     {
 706         AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
 707         for (int d = 0; d < awhBiasParams->ndim; d++)
 708         {
 709             AwhDimParams* dimParams  = &awhBiasParams->dimParams[d];
 710             int           coordIndex = dimParams->coordIndex;
 711             double origin = dimParams->origin, end = dimParams->end, period = dimParams->period;
 712             double coordValueInit = dimParams->coordValueInit;
 713
 714             if ((period == 0) && (origin > end))
 715             {
 716                 gmx_fatal(FARGS,
 717                           "For the non-periodic pull coordinates awh%d-dim%d-start (%f) cannot be "
 718                           "larger than awh%d-dim%d-end (%f)",
 719                           k + 1, d + 1, origin, k + 1, d + 1, end);
 720             }
 721
 722             /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
 723                Make sure the AWH interval is within this reference interval.
 724
 725                Note: we could fairly simply allow using a  more general interval (e.g. [x, x + period]) but it complicates
 726                things slightly and I don't see that there is a great need for it. It would also mean that the interval would
 727                depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180,
 728                independent of AWH parameters.
 729              */
 730             if (!intervalIsInPeriodicInterval(origin, end, period))
 731             {
 732                 gmx_fatal(FARGS,
 733                           "When using AWH with periodic pull coordinate geometries "
 734                           "awh%d-dim%d-start (%.8g) and "
 735                           "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take "
 736                           "values in between "
 737                           "minus half a period and plus half a period, i.e. in the interval [%.8g, "
 738                           "%.8g].",
 739                           k + 1, d + 1, origin, k + 1, d + 1, end, period, -0.5 * period, 0.5 * period);
 740             }
 741
 742             /* Warn if the pull initial coordinate value is not in the grid */
 743             if (!valueIsInInterval(origin, end, period, coordValueInit))
 744             {
 745                 auto message = formatString(
 746                         "The initial coordinate value (%.8g) for pull coordinate index %d falls "
 747                         "outside "
 748                         "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end "
 749                         "(%.8g). "
 750                         "This can lead to large initial forces pulling the coordinate towards the "
 751                         "sampling interval.",
 752                         coordValueInit, coordIndex + 1, k + 1, d + 1, origin, k + 1, d + 1, end);
 753                 warning(wi, message);
 754             }
 755         }
 756     }
 757 }
 758
 759 void setStateDependentAwhParams(AwhParams*           awhParams,
 760                                 const pull_params_t* pull_params,
 761                                 pull_t*              pull_work,
 762                                 const matrix         box,
 763                                 int                  ePBC,
 764                                 const tensor&        compressibility,
 765                                 const t_grpopts*     inputrecGroupOptions,
 766                                 warninp_t            wi)
 767 {
 768     /* The temperature is not really state depenendent but is not known
 769      * when read_awhParams is called (in get ir).
 770      * It is known first after do_index has been called in grompp.cpp.
 771      */
 772     if (inputrecGroupOptions->ref_t == nullptr || inputrecGroupOptions->ref_t[0] <= 0)
 773     {
 774         gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
 775     }
 776     for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
 777     {
 778         if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
 779         {
 780             gmx_fatal(FARGS,
 781                       "AWH biasing is currently only supported for identical temperatures for all "
 782                       "temperature coupling groups");
 783         }
 784     }
 785
 786     t_pbc pbc;
 787     set_pbc(&pbc, ePBC, box);
 788
 789     for (int k = 0; k < awhParams->numBias; k++)
 790     {
 791         AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
 792         for (int d = 0; d < awhBiasParams->ndim; d++)
 793         {
 794             AwhDimParams*       dimParams       = &awhBiasParams->dimParams[d];
 795             const t_pull_coord& pullCoordParams = pull_params->coord[dimParams->coordIndex];
 796
 797             if (pullCoordParams.eGeom == epullgDIRPBC)
 798             {
 799                 gmx_fatal(FARGS,
 800                           "AWH does not support pull geometry '%s'. "
 801                           "If the maximum distance between the groups is always less than half the "
 802                           "box size, "
 803                           "you can use geometry '%s' instead.",
 804                           EPULLGEOM(epullgDIRPBC), EPULLGEOM(epullgDIR));
 805             }
 806
 807             dimParams->period =
 808                     get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin);
 809             // We would like to check for scaling, but we don't have the full inputrec available here
 810             if (dimParams->period > 0
 811                 && !(pullCoordParams.eGeom == epullgANGLE || pullCoordParams.eGeom == epullgDIHEDRAL))
 812             {
 813                 bool coordIsScaled = false;
 814                 for (int d2 = 0; d2 < DIM; d2++)
 815                 {
 816                     if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0)
 817                     {
 818                         coordIsScaled = true;
 819                     }
 820                 }
 821                 if (coordIsScaled)
 822                 {
 823                     std::string mesg = gmx::formatString(
 824                             "AWH dimension %d of bias %d is periodic with pull geometry '%s', "
 825                             "while you should are applying pressure scaling to the corresponding "
 826                             "box vector, this is not supported.",
 827                             d + 1, k + 1, EPULLGEOM(pullCoordParams.eGeom));
 828                     warning(wi, mesg.c_str());
 829                 }
 830             }
 831
 832             /* The initial coordinate value, converted to external user units. */
 833             dimParams->coordValueInit = get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
 834
 835             dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(&pullCoordParams);
 836         }
 837     }
 838     checkInputConsistencyInterval(awhParams, wi);
 839
 840     /* Register AWH as external potential with pull to check consistency. */
 841     Awh::registerAwhWithPull(*awhParams, pull_work);
 842 }
 843
 844 } // namespace gmx