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