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