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