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