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] pull_params      Pull parameters.
 557  * \param[in] coord_ind        Pull coordinate index.
 558  * \param[in] box              Box vectors.
 559  * \returns the period (or 0 if not periodic).
 560  */
 561 static double get_pull_coord_period(const pull_params_t *pull_params,
 562                                     int                  coord_ind,
 563                                     const matrix         box)
 564 {
 565     double        period;
 566     t_pull_coord *pcrd_params = &pull_params->coord[coord_ind];
 567
 568     if (pcrd_params->eGeom == epullgDIRPBC)
 569     {
 570         /* For direction periodic, we need the pull vector to be one of the box vectors
 571            (or more generally I guess it could be an integer combination of boxvectors).
 572            This boxvector should to be orthogonal to the (periodic) plane spanned by the other two box vectors.
 573            Here we assume that the pull vector is either x, y or z.
 574          * E.g. for pull vec = (1, 0, 0) the box vector tensor should look like:
 575          * | x 0 0 |
 576          * | 0 a c |
 577          * | 0 b d |
 578          *
 579            The period is then given by the box length x.
 580
 581            Note: we make these checks here for AWH and not in pull because we allow pull to be more general.
 582          */
 583         int m_pullvec = -1, count_nonzeros = 0;
 584
 585         /* Check that pull vec has only one component and which component it is. This component gives the relevant box vector */
 586         for (int m = 0; m < DIM; m++)
 587         {
 588             if (pcrd_params->vec[m] != 0)
 589             {
 590                 m_pullvec = m;
 591                 count_nonzeros++;
 592             }
 593         }
 594         if (count_nonzeros != 1)
 595         {
 596             gmx_fatal(FARGS, "For AWH biasing pull coordinate %d with pull geometry %s, the pull vector needs to be parallel to "
 597                       "a box vector that is parallel to either the x, y or z axis and is orthogonal to the other box vectors.",
 598                       coord_ind + 1, EPULLGEOM(epullgDIRPBC));
 599         }
 600
 601         /* Check that there is a box vec parallel to pull vec and that this boxvec is orthogonal to the other box vectors */
 602         for (int m = 0; m < DIM; m++)
 603         {
 604             for (int n = 0; n < DIM; n++)
 605             {
 606                 if ((n != m) && (n == m_pullvec || m == m_pullvec) && box[m][n] > 0)
 607                 {
 608                     gmx_fatal(FARGS, "For AWH biasing pull coordinate %d with pull geometry %s, there needs to be a box vector parallel to the pull vector that is "
 609                               "orthogonal to the other box vectors.",
 610                               coord_ind + 1, EPULLGEOM(epullgDIRPBC));
 611                 }
 612             }
 613         }
 614
 615         /* If this box vector only has one component as we assumed the norm should be equal to the absolute value of that component */
 616         period = static_cast<double>(norm(box[m_pullvec]));
 617     }
 618     else if (pcrd_params->eGeom == epullgDIHEDRAL)
 619     {
 620         /* The dihedral angle is periodic in -180 to 180 deg */
 621         period = 360;
 622     }
 623     else
 624     {
 625         period = 0;
 626     }
 627
 628     return period;
 629 }
 630
 631 /*! \brief
 632  * Checks if the given interval is defined in the correct periodic interval.
 633  *
 634  * \param[in] origin      Start value of interval.
 635  * \param[in] end         End value of interval.
 636  * \param[in] period      Period (or 0 if not periodic).
 637  * \returns true if the end point values are in the correct periodic interval.
 638  */
 639 static bool intervalIsInPeriodicInterval(double origin, double end, double period)
 640 {
 641     return (period == 0) || (std::fabs(origin) <= 0.5*period && std::fabs(end) <= 0.5*period);
 642 }
 643
 644 /*! \brief
 645  * Checks if a value is within an interval.
 646  *
 647  * \param[in] origin      Start value of interval.
 648  * \param[in] end         End value of interval.
 649  * \param[in] period      Period (or 0 if not periodic).
 650  * \param[in] value       Value to check.
 651  * \returns true if the value is within the interval.
 652  */
 653 static bool valueIsInInterval(double origin, double end, double period, double value)
 654 {
 655     bool bIn_interval;
 656
 657     if (period > 0)
 658     {
 659         if (origin < end)
 660         {
 661             /* The interval closes within the periodic interval */
 662             bIn_interval = (value >= origin) && (value <= end);
 663         }
 664         else
 665         {
 666             /* The interval wraps around the periodic boundary */
 667             bIn_interval = ((value >= origin) && (value <= 0.5*period)) || ((value >= -0.5*period) && (value <= end));
 668         }
 669     }
 670     else
 671     {
 672         bIn_interval = (value >= origin) && (value <= end);
 673     }
 674
 675     return bIn_interval;
 676 }
 677
 678 /*! \brief
 679  * Check if the starting configuration is consistent with the given interval.
 680  *
 681  * \param[in]     awhParams  AWH parameters.
 682  * \param[in,out] wi         Struct for bookeeping warnings.
 683  */
 684 static void checkInputConsistencyInterval(const AwhParams *awhParams, warninp_t wi)
 685 {
 686     for (int k = 0; k < awhParams->numBias; k++)
 687     {
 688         AwhBiasParams    *awhBiasParams  = &awhParams->awhBiasParams[k];
 689         for (int d = 0; d < awhBiasParams->ndim; d++)
 690         {
 691             AwhDimParams *dimParams      = &awhBiasParams->dimParams[d];
 692             int           coordIndex     = dimParams->coordIndex;
 693             double        origin         = dimParams->origin, end = dimParams->end, period = dimParams->period;
 694             double        coordValueInit = dimParams->coordValueInit;
 695
 696             if ((period == 0) && (origin > end))
 697             {
 698                 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)",
 699                           k + 1, d + 1, origin, k + 1, d + 1, end);
 700             }
 701
 702             /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
 703                Make sure the AWH interval is within this reference interval.
 704
 705                Note: we could fairly simply allow using a  more general interval (e.g. [x, x + period]) but it complicates
 706                things slightly and I don't see that there is a great need for it. It would also mean that the interval would
 707                depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180,
 708                independent of AWH parameters.
 709              */
 710             if (!intervalIsInPeriodicInterval(origin, end, period))
 711             {
 712                 gmx_fatal(FARGS, "When using AWH with periodic pull coordinate geometries awh%d-dim%d-start (%.8g) and "
 713                           "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take values in between "
 714                           "minus half a period and plus half a period, i.e. in the interval [%.8g, %.8g].",
 715                           k + 1, d + 1, origin, k + 1, d + 1, end,
 716                           period, -0.5*period, 0.5*period);
 717
 718             }
 719
 720             /* Warn if the pull initial coordinate value is not in the grid */
 721             if (!valueIsInInterval(origin, end, period, coordValueInit))
 722             {
 723                 auto message = formatString
 724                         ("The initial coordinate value (%.8g) for pull coordinate index %d falls outside "
 725                         "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end (%.8g). "
 726                         "This can lead to large initial forces pulling the coordinate towards the sampling interval.",
 727                         coordValueInit, coordIndex + 1, k + 1, d + 1, origin, k + 1, d + 1, end);
 728                 warning(wi, message);
 729             }
 730         }
 731     }
 732 }
 733
 734 void setStateDependentAwhParams(AwhParams *awhParams,
 735                                 const pull_params_t *pull_params, pull_t *pull_work,
 736                                 const matrix box,  int ePBC,
 737                                 const t_grpopts *inputrecGroupOptions, warninp_t wi)
 738 {
 739     /* The temperature is not really state depenendent but is not known
 740      * when read_awhParams is called (in get ir).
 741      * It is known first after do_index has been called in grompp.cpp.
 742      */
 743     if (inputrecGroupOptions->ref_t == nullptr ||
 744         inputrecGroupOptions->ref_t[0] <= 0)
 745     {
 746         gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
 747     }
 748     for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
 749     {
 750         if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
 751         {
 752             gmx_fatal(FARGS, "AWH biasing is currently only supported for identical temperatures for all temperature coupling groups");
 753         }
 754     }
 755
 756     t_pbc          pbc;
 757     set_pbc(&pbc, ePBC, box);
 758
 759     for (int k = 0; k < awhParams->numBias; k++)
 760     {
 761         AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
 762         for (int d = 0; d < awhBiasParams->ndim; d++)
 763         {
 764             AwhDimParams *dimParams = &awhBiasParams->dimParams[d];
 765
 766             /* The periodiciy of the AWH grid in certain cases depends on the simulation box */
 767             dimParams->period = get_pull_coord_period(pull_params, dimParams->coordIndex, box);
 768
 769             /* The initial coordinate value, converted to external user units. */
 770             dimParams->coordValueInit =
 771                 get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
 772
 773             t_pull_coord *pullCoord = &pull_params->coord[dimParams->coordIndex];
 774             dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(pullCoord);
 775         }
 776     }
 777     checkInputConsistencyInterval(awhParams, wi);
 778
 779     /* Register AWH as external potential with pull to check consistency. */
 780     Awh::registerAwhWithPull(*awhParams, pull_work);
 781 }
 782
 783 } // namespace gmx