2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
39 #include "read-params.h"
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"
59 #include "biasparams.h"
60 #include "biassharing.h"
65 const char *eawhtarget_names[eawhtargetNR+1] = {
66 "constant", "cutoff", "boltzmann", "local-boltzmann", nullptr
69 const char *eawhgrowth_names[eawhgrowthNR+1] = {
70 "exp-linear", "linear", nullptr
73 const char *eawhpotential_names[eawhpotentialNR+1] = {
74 "convolved", "umbrella", nullptr
77 const char *eawhcoordprovider_names[eawhcoordproviderNR+1] = {
82 * Read parameters of an AWH bias dimension.
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.
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)
98 printStringNoNewline(inp, "The provider of the reaction coordinate, currently only pull is supported");
101 opt = prefix + "-coord-provider";
102 dimParams->eCoordProvider = get_eeenum(inp, opt, eawhcoordprovider_names, wi);
106 printStringNoNewline(inp, "The coordinate index for this dimension");
108 opt = prefix + "-coord-index";
110 coordIndexInput = get_eint(inp, opt, 1, wi);
111 if (coordIndexInput < 1)
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());
117 /* The pull coordinate indices start at 1 in the input file, at 0 internally */
118 dimParams->coordIndex = coordIndexInput - 1;
120 /* The pull settings need to be consistent with the AWH settings */
121 if (!(pull_params->coord[dimParams->coordIndex].eType == epullEXTERNAL) )
123 gmx_fatal(FARGS, "AWH biasing can only be applied to pull type %s",
124 EPULLTYPE(epullEXTERNAL));
127 if (dimParams->coordIndex >= pull_params->ncoord)
129 gmx_fatal(FARGS, "The given AWH coordinate index (%d) is larger than the number of pull coordinates (%d)",
130 coordIndexInput, pull_params->ncoord);
132 if (pull_params->coord[dimParams->coordIndex].rate != 0)
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);
139 /* Grid params for each axis */
140 int eGeom = pull_params->coord[dimParams->coordIndex].eGeom;
144 printStringNoNewline(inp, "Start and end values for each coordinate dimension");
147 opt = prefix + "-start";
148 dimParams->origin = get_ereal(inp, opt, 0., wi);
150 opt = prefix + "-end";
151 dimParams->end = get_ereal(inp, opt, 0., wi);
153 if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
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);
160 /* Check that the requested interval is in allowed range */
161 if (eGeom == epullgDIST)
163 if (dimParams->origin < 0 || dimParams->end < 0)
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));
170 else if (eGeom == epullgANGLE || eGeom == epullgANGLEAXIS)
172 if (dimParams->origin < 0 || dimParams->end > 180)
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));
178 else if (eGeom == epullgDIHEDRAL)
180 if (dimParams->origin < -180 || dimParams->end > 180)
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));
189 printStringNoNewline(inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
191 opt = prefix + "-force-constant";
192 dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
193 if (dimParams->forceConstant <= 0)
195 warning_error(wi, "The force AWH bias force constant should be > 0");
200 printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps)");
202 opt = prefix + "-diffusion";
203 dimParams->diffusion = get_ereal(inp, opt, 0, wi);
205 if (dimParams->diffusion <= 0)
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;
218 printStringNoNewline(inp, "Diameter that needs to be sampled around a point before it is considered covered.");
220 opt = prefix + "-cover-diameter";
221 dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
223 if (dimParams->coverDiameter < 0)
225 gmx_fatal(FARGS, "%s (%g) cannot be negative.",
226 opt.c_str(), dimParams->coverDiameter);
231 * Check consistency of input at the AWH bias level.
233 * \param[in] awhBiasParams AWH bias parameters.
234 * \param[in,out] wi Struct for bookkeeping warnings.
236 static void checkInputConsistencyAwhBias(const AwhBiasParams &awhBiasParams,
239 /* Covering diameter and sharing warning. */
240 for (int d = 0; d < awhBiasParams.ndim; d++)
242 double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
243 if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
245 warning(wi, "The covering diameter is only relevant to set for bias sharing simulations.");
251 * Read parameters of an AWH bias.
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.
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)
265 printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
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)
273 gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
278 printStringNoNewline(inp, "Growth rate of the reference histogram determining the bias update size: exp-linear or linear");
280 opt = prefix + "-growth";
281 awhBiasParams->eGrowth = get_eeenum(inp, opt, eawhgrowth_names, wi);
285 printStringNoNewline(inp, "Start the simulation by equilibrating histogram towards the target distribution: no or yes");
287 opt = prefix + "-equilibrate-histogram";
288 awhBiasParams->equilibrateHistogram = (get_eeenum(inp, opt, yesno_names, wi) != 0);
289 if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != eawhgrowthEXP_LINEAR)
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);
298 printStringNoNewline(inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
300 opt = prefix + "-target";
301 awhBiasParams->eTarget = get_eeenum(inp, opt, eawhtarget_names, wi);
303 if ((awhBiasParams->eTarget == eawhtargetLOCALBOLTZMANN) &&
304 (awhBiasParams->eGrowth == eawhgrowthEXP_LINEAR))
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 "
309 EAWHTARGET(eawhtargetLOCALBOLTZMANN), EAWHGROWTH(eawhgrowthEXP_LINEAR),
310 EAWHGROWTH(eawhgrowthLINEAR));
311 warning(wi, message);
316 printStringNoNewline(inp, "Boltzmann beta scaling factor for target distribution types 'boltzmann' and 'boltzmann-local'");
318 opt = prefix + "-target-beta-scaling";
319 awhBiasParams->targetBetaScaling = get_ereal(inp, opt, 0, wi);
321 switch (awhBiasParams->eTarget)
323 case eawhtargetBOLTZMANN:
324 case eawhtargetLOCALBOLTZMANN:
325 if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
327 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.",
328 opt.c_str(), awhBiasParams->targetBetaScaling, EAWHTARGET(awhBiasParams->eTarget));
332 if (awhBiasParams->targetBetaScaling != 0)
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));
342 printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
344 opt = prefix + "-target-cutoff";
345 awhBiasParams->targetCutoff = get_ereal(inp, opt, 0, wi);
347 switch (awhBiasParams->eTarget)
349 case eawhtargetCUTOFF:
350 if (awhBiasParams->targetCutoff <= 0)
352 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.",
353 opt.c_str(), awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
357 if (awhBiasParams->targetCutoff != 0)
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));
367 printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
369 opt = prefix + "-user-data";
370 awhBiasParams->bUserData = get_eeenum(inp, opt, yesno_names, wi);
374 printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
376 opt = prefix + "-share-group";
377 awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
378 if (awhBiasParams->shareGroup < 0)
380 warning_error(wi, "AWH bias share-group should be >= 0");
385 printStringNoNewline(inp, "Dimensionality of the coordinate");
387 opt = prefix + "-ndim";
388 awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
390 if (awhBiasParams->ndim <= 0 ||
391 awhBiasParams->ndim > c_biasMaxNumDim)
393 gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(), awhBiasParams->ndim, c_biasMaxNumDim);
395 if (awhBiasParams->ndim > 2)
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!");
400 snew(awhBiasParams->dimParams, awhBiasParams->ndim);
401 for (int d = 0; d < awhBiasParams->ndim; d++)
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);
408 /* Check consistencies here that cannot be checked at read time at a lower level. */
409 checkInputConsistencyAwhBias(*awhBiasParams, wi);
413 * Check consistency of input at the AWH level.
415 * \param[in] awhParams AWH parameters.
416 * \param[in,out] wi Struct for bookkeeping warnings.
418 static void checkInputConsistencyAwh(const AwhParams &awhParams,
421 /* Each pull coord can map to at most 1 AWH coord.
422 * Check that we have a shared bias when requesting multisim sharing.
424 bool haveSharedBias = false;
425 for (int k1 = 0; k1 < awhParams.numBias; k1++)
427 const AwhBiasParams &awhBiasParams1 = awhParams.awhBiasParams[k1];
429 if (awhBiasParams1.shareGroup > 0)
431 haveSharedBias = true;
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++)
437 for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
439 const AwhBiasParams &awhBiasParams2 = awhParams.awhBiasParams[k2];
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++)
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) )
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);
458 if (awhParams.shareBiasMultisim && !haveSharedBias)
460 warning(wi, "Sharing of biases over multiple simulations is requested, but no bias is marked as shared (share-group > 0)");
463 /* mdrun does not support this (yet), but will check again */
464 if (haveBiasSharingWithinSimulation(awhParams))
466 warning(wi, "You have shared biases within a single simulation, but mdrun does not support this (yet)");
470 AwhParams *readAndCheckAwhParams(std::vector<t_inpfile> *inp, const t_inputrec *ir, warninp_t wi)
472 AwhParams *awhParams;
476 /* Parameters common for all biases */
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);
482 printStringNoNewline(inp, "The random seed used for sampling the umbrella center in the case of umbrella type potential");
484 awhParams->seed = get_eint(inp, opt, -1, wi);
485 if (awhParams->seed == -1)
487 awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
488 fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
491 printStringNoNewline(inp, "Data output interval in number of steps");
493 awhParams->nstOut = get_eint(inp, opt, 100000, wi);
494 if (awhParams->nstOut <= 0)
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);
500 /* This restriction can be removed by changing a flag of print_ebin() */
501 if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
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);
508 printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
509 opt = "awh-nstsample";
510 awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
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)
517 warning_error(wi, opt + " needs to be an integer > 0");
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);
524 printStringNoNewline(inp, "The number of independent AWH biases");
526 awhParams->numBias = get_eint(inp, opt, 1, wi);
527 if (awhParams->numBias <= 0)
529 gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
532 /* Read the parameters specific to each AWH bias */
533 snew(awhParams->awhBiasParams, awhParams->numBias);
535 for (int k = 0; k < awhParams->numBias; k++)
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);
542 /* Do a final consistency check before returning */
543 checkInputConsistencyAwh(*awhParams, wi);
545 if (ir->init_step != 0)
547 warning_error(wi, "With AWH init-step should be 0");
554 * Gets the period of a pull coordinate.
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).
561 static double get_pull_coord_period(const pull_params_t *pull_params,
566 t_pull_coord *pcrd_params = &pull_params->coord[coord_ind];
568 if (pcrd_params->eGeom == epullgDIRPBC)
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:
579 The period is then given by the box length x.
581 Note: we make these checks here for AWH and not in pull because we allow pull to be more general.
583 int m_pullvec = -1, count_nonzeros = 0;
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++)
588 if (pcrd_params->vec[m] != 0)
594 if (count_nonzeros != 1)
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));
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++)
604 for (int n = 0; n < DIM; n++)
606 if ((n != m) && (n == m_pullvec || m == m_pullvec) && box[m][n] > 0)
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));
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]));
618 else if (pcrd_params->eGeom == epullgDIHEDRAL)
620 /* The dihedral angle is periodic in -180 to 180 deg */
632 * Checks if the given interval is defined in the correct periodic interval.
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.
639 static bool intervalIsInPeriodicInterval(double origin, double end, double period)
641 return (period == 0) || (std::fabs(origin) <= 0.5*period && std::fabs(end) <= 0.5*period);
645 * Checks if a value is within an interval.
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.
653 static bool valueIsInInterval(double origin, double end, double period, double value)
661 /* The interval closes within the periodic interval */
662 bIn_interval = (value >= origin) && (value <= end);
666 /* The interval wraps around the periodic boundary */
667 bIn_interval = ((value >= origin) && (value <= 0.5*period)) || ((value >= -0.5*period) && (value <= end));
672 bIn_interval = (value >= origin) && (value <= end);
679 * Check if the starting configuration is consistent with the given interval.
681 * \param[in] awhParams AWH parameters.
682 * \param[in,out] wi Struct for bookeeping warnings.
684 static void checkInputConsistencyInterval(const AwhParams *awhParams, warninp_t wi)
686 for (int k = 0; k < awhParams->numBias; k++)
688 AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
689 for (int d = 0; d < awhBiasParams->ndim; d++)
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;
696 if ((period == 0) && (origin > end))
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);
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.
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.
710 if (!intervalIsInPeriodicInterval(origin, end, period))
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);
720 /* Warn if the pull initial coordinate value is not in the grid */
721 if (!valueIsInInterval(origin, end, period, coordValueInit))
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);
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)
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.
743 if (inputrecGroupOptions->ref_t == nullptr ||
744 inputrecGroupOptions->ref_t[0] <= 0)
746 gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
748 for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
750 if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
752 gmx_fatal(FARGS, "AWH biasing is currently only supported for identical temperatures for all temperature coupling groups");
757 set_pbc(&pbc, ePBC, box);
759 for (int k = 0; k < awhParams->numBias; k++)
761 AwhBiasParams *awhBiasParams = &awhParams->awhBiasParams[k];
762 for (int d = 0; d < awhBiasParams->ndim; d++)
764 AwhDimParams *dimParams = &awhBiasParams->dimParams[d];
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);
769 /* The initial coordinate value, converted to external user units. */
770 dimParams->coordValueInit =
771 get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
773 t_pull_coord *pullCoord = &pull_params->coord[dimParams->coordIndex];
774 dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(pullCoord);
777 checkInputConsistencyInterval(awhParams, wi);
779 /* Register AWH as external potential with pull to check consistency. */
780 Awh::registerAwhWithPull(*awhParams, pull_work);