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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
40 #include "read_params.h"
42 #include "gromacs/awh/awh.h"
43 #include "gromacs/fileio/readinp.h"
44 #include "gromacs/fileio/warninp.h"
45 #include "gromacs/math/units.h"
46 #include "gromacs/math/utilities.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/mdtypes/awh_params.h"
49 #include "gromacs/mdtypes/inputrec.h"
50 #include "gromacs/mdtypes/md_enums.h"
51 #include "gromacs/mdtypes/pull_params.h"
52 #include "gromacs/pbcutil/pbc.h"
53 #include "gromacs/pulling/pull.h"
54 #include "gromacs/random/seed.h"
55 #include "gromacs/utility/cstringutil.h"
56 #include "gromacs/utility/fatalerror.h"
57 #include "gromacs/utility/smalloc.h"
58 #include "gromacs/utility/stringutil.h"
60 #include "biasparams.h"
61 #include "biassharing.h"
66 const char* eawhtarget_names[eawhtargetNR + 1] = { "constant", "cutoff", "boltzmann",
67 "local-boltzmann", nullptr };
69 const char* eawhgrowth_names[eawhgrowthNR + 1] = { "exp-linear", "linear", nullptr };
71 const char* eawhpotential_names[eawhpotentialNR + 1] = { "convolved", "umbrella", nullptr };
73 const char* eawhcoordprovider_names[eawhcoordproviderNR + 1] = { "pull", nullptr };
76 * Read parameters of an AWH bias dimension.
78 * \param[in,out] inp Input file entries.
79 * \param[in] prefix Prefix for dimension parameters.
80 * \param[in,out] dimParams AWH dimensional parameters.
81 * \param[in] pull_params Pull parameters.
82 * \param[in,out] wi Struct for bookeeping warnings.
83 * \param[in] bComment True if comments should be printed.
85 static void readDimParams(std::vector<t_inpfile>* inp,
86 const std::string& prefix,
87 AwhDimParams* dimParams,
88 const pull_params_t* pull_params,
96 inp, "The provider of the reaction coordinate, currently only pull is supported");
99 opt = prefix + "-coord-provider";
100 dimParams->eCoordProvider = get_eeenum(inp, opt, eawhcoordprovider_names, wi);
104 printStringNoNewline(inp, "The coordinate index for this dimension");
106 opt = prefix + "-coord-index";
108 coordIndexInput = get_eint(inp, opt, 1, wi);
109 if (coordIndexInput < 1)
112 "Failed to read a valid coordinate index for %s. "
113 "Note that the pull coordinate indexing starts at 1.",
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", EPULLTYPE(epullEXTERNAL));
126 if (dimParams->coordIndex >= pull_params->ncoord)
129 "The given AWH coordinate index (%d) is larger than the number of pull "
131 coordIndexInput, pull_params->ncoord);
133 if (pull_params->coord[dimParams->coordIndex].rate != 0)
135 auto message = formatString(
136 "Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate",
137 coordIndexInput, pull_params->coord[dimParams->coordIndex].rate);
138 warning_error(wi, message);
141 /* Grid params for each axis */
142 int eGeom = pull_params->coord[dimParams->coordIndex].eGeom;
146 printStringNoNewline(inp, "Start and end values for each coordinate dimension");
149 opt = prefix + "-start";
150 dimParams->origin = get_ereal(inp, opt, 0., wi);
152 opt = prefix + "-end";
153 dimParams->end = get_ereal(inp, opt, 0., wi);
155 if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
157 auto message = formatString(
158 "The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
159 "This will result in only one point along this axis in the coordinate value grid.",
160 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end);
161 warning(wi, message);
163 /* Check that the requested interval is in allowed range */
164 if (eGeom == epullgDIST)
166 if (dimParams->origin < 0 || dimParams->end < 0)
169 "%s-start (%g) or %s-end (%g) set to a negative value. With pull geometry "
170 "distance coordinate values are non-negative. "
171 "Perhaps you want to use geometry %s instead?",
172 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
173 EPULLGEOM(epullgDIR));
176 else if (eGeom == epullgANGLE || eGeom == epullgANGLEAXIS)
178 if (dimParams->origin < 0 || dimParams->end > 180)
181 "%s-start (%g) and %s-end (%g) are outside of the allowed range 0 to 180 deg "
182 "for pull geometries %s and %s ",
183 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
184 EPULLGEOM(epullgANGLE), EPULLGEOM(epullgANGLEAXIS));
187 else if (eGeom == epullgDIHEDRAL)
189 if (dimParams->origin < -180 || dimParams->end > 180)
192 "%s-start (%g) and %s-end (%g) are outside of the allowed range -180 to 180 "
193 "deg for pull geometry %s. ",
194 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
195 EPULLGEOM(epullgDIHEDRAL));
201 printStringNoNewline(
202 inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
204 opt = prefix + "-force-constant";
205 dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
206 if (dimParams->forceConstant <= 0)
208 warning_error(wi, "The force AWH bias force constant should be > 0");
213 printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps)");
215 opt = prefix + "-diffusion";
216 dimParams->diffusion = get_ereal(inp, opt, 0, wi);
218 if (dimParams->diffusion <= 0)
220 const double diffusion_default = 1e-5;
221 auto message = formatString(
222 "%s not explicitly set by user. You can choose to use a default "
223 "value (%g nm^2/ps or rad^2/ps) but this may very well be "
224 "non-optimal for your system!",
225 opt.c_str(), diffusion_default);
226 warning(wi, message);
227 dimParams->diffusion = diffusion_default;
232 printStringNoNewline(inp,
233 "Diameter that needs to be sampled around a point before it is "
234 "considered covered.");
236 opt = prefix + "-cover-diameter";
237 dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
239 if (dimParams->coverDiameter < 0)
241 gmx_fatal(FARGS, "%s (%g) cannot be negative.", opt.c_str(), dimParams->coverDiameter);
246 * Check consistency of input at the AWH bias level.
248 * \param[in] awhBiasParams AWH bias parameters.
249 * \param[in,out] wi Struct for bookkeeping warnings.
251 static void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
253 /* Covering diameter and sharing warning. */
254 for (int d = 0; d < awhBiasParams.ndim; d++)
256 double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
257 if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
260 "The covering diameter is only relevant to set for bias sharing simulations.");
266 * Read parameters of an AWH bias.
268 * \param[in,out] inp Input file entries.
269 * \param[in,out] awhBiasParams AWH dimensional parameters.
270 * \param[in] prefix Prefix for bias parameters.
271 * \param[in] ir Input parameter struct.
272 * \param[in,out] wi Struct for bookeeping warnings.
273 * \param[in] bComment True if comments should be printed.
275 static void read_bias_params(std::vector<t_inpfile>* inp,
276 AwhBiasParams* awhBiasParams,
277 const std::string& prefix,
278 const t_inputrec* ir,
284 printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
287 std::string opt = prefix + "-error-init";
288 /* We allow using a default value here without warning (but warn the user if the diffusion constant is not set). */
289 awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi);
290 if (awhBiasParams->errorInitial <= 0)
292 gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
297 printStringNoNewline(inp,
298 "Growth rate of the reference histogram determining the bias update "
299 "size: exp-linear or linear");
301 opt = prefix + "-growth";
302 awhBiasParams->eGrowth = get_eeenum(inp, opt, eawhgrowth_names, wi);
306 printStringNoNewline(inp,
307 "Start the simulation by equilibrating histogram towards the target "
308 "distribution: no or yes");
310 opt = prefix + "-equilibrate-histogram";
311 awhBiasParams->equilibrateHistogram = (get_eeenum(inp, opt, yesno_names, wi) != 0);
312 if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != eawhgrowthEXP_LINEAR)
315 formatString("Option %s will only have an effect for histogram growth type '%s'.",
316 opt.c_str(), EAWHGROWTH(eawhgrowthEXP_LINEAR));
317 warning(wi, message);
322 printStringNoNewline(
323 inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
325 opt = prefix + "-target";
326 awhBiasParams->eTarget = get_eeenum(inp, opt, eawhtarget_names, wi);
328 if ((awhBiasParams->eTarget == eawhtargetLOCALBOLTZMANN)
329 && (awhBiasParams->eGrowth == eawhgrowthEXP_LINEAR))
331 auto message = formatString(
332 "Target type '%s' combined with histogram growth type '%s' is not "
333 "expected to give stable bias updates. You probably want to use growth type "
335 EAWHTARGET(eawhtargetLOCALBOLTZMANN), EAWHGROWTH(eawhgrowthEXP_LINEAR),
336 EAWHGROWTH(eawhgrowthLINEAR));
337 warning(wi, message);
342 printStringNoNewline(inp,
343 "Boltzmann beta scaling factor for target distribution types "
344 "'boltzmann' and 'boltzmann-local'");
346 opt = prefix + "-target-beta-scaling";
347 awhBiasParams->targetBetaScaling = get_ereal(inp, opt, 0, wi);
349 switch (awhBiasParams->eTarget)
351 case eawhtargetBOLTZMANN:
352 case eawhtargetLOCALBOLTZMANN:
353 if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
355 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
356 awhBiasParams->targetBetaScaling, EAWHTARGET(awhBiasParams->eTarget));
360 if (awhBiasParams->targetBetaScaling != 0)
364 "Value for %s (%g) set explicitly but will not be used for target type %s.",
365 opt.c_str(), awhBiasParams->targetBetaScaling,
366 EAWHTARGET(awhBiasParams->eTarget));
373 printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
375 opt = prefix + "-target-cutoff";
376 awhBiasParams->targetCutoff = get_ereal(inp, opt, 0, wi);
378 switch (awhBiasParams->eTarget)
380 case eawhtargetCUTOFF:
381 if (awhBiasParams->targetCutoff <= 0)
383 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
384 awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
388 if (awhBiasParams->targetCutoff != 0)
392 "Value for %s (%g) set explicitly but will not be used for target type %s.",
393 opt.c_str(), awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
400 printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
402 opt = prefix + "-user-data";
403 awhBiasParams->bUserData = get_eeenum(inp, opt, yesno_names, wi);
407 printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
409 opt = prefix + "-share-group";
410 awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
411 if (awhBiasParams->shareGroup < 0)
413 warning_error(wi, "AWH bias share-group should be >= 0");
418 printStringNoNewline(inp, "Dimensionality of the coordinate");
420 opt = prefix + "-ndim";
421 awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
423 if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
425 gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(),
426 awhBiasParams->ndim, c_biasMaxNumDim);
428 if (awhBiasParams->ndim > 2)
431 "For awh-dim > 2 the estimate based on the diffusion and the initial error is "
432 "currently only a rough guideline."
433 " You should verify its usefulness for your system before production runs!");
435 snew(awhBiasParams->dimParams, awhBiasParams->ndim);
436 for (int d = 0; d < awhBiasParams->ndim; d++)
438 bComment = bComment && d == 0;
439 std::string prefixdim = prefix + formatString("-dim%d", d + 1);
440 readDimParams(inp, prefixdim, &awhBiasParams->dimParams[d], ir->pull, wi, bComment);
443 /* Check consistencies here that cannot be checked at read time at a lower level. */
444 checkInputConsistencyAwhBias(*awhBiasParams, wi);
448 * Check consistency of input at the AWH level.
450 * \param[in] awhParams AWH parameters.
451 * \param[in,out] wi Struct for bookkeeping warnings.
453 static void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi)
455 /* Each pull coord can map to at most 1 AWH coord.
456 * Check that we have a shared bias when requesting multisim sharing.
458 bool haveSharedBias = false;
459 for (int k1 = 0; k1 < awhParams.numBias; k1++)
461 const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1];
463 if (awhBiasParams1.shareGroup > 0)
465 haveSharedBias = true;
468 /* k1 is the reference AWH, k2 is the AWH we compare with (can be equal to k1) */
469 for (int k2 = k1; k2 < awhParams.numBias; k2++)
471 for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
473 const AwhBiasParams& awhBiasParams2 = awhParams.awhBiasParams[k2];
475 /* d1 is the reference dimension of the reference AWH. d2 is the dim index of the AWH to compare with. */
476 for (int d2 = 0; d2 < awhBiasParams2.ndim; d2++)
478 /* Give an error if (d1, k1) is different from (d2, k2) but the pull coordinate is the same */
479 if ((d1 != d2 || k1 != k2)
480 && (awhBiasParams1.dimParams[d1].coordIndex == awhBiasParams2.dimParams[d2].coordIndex))
482 char errormsg[STRLEN];
484 "One pull coordinate (%d) cannot be mapped to two separate AWH "
485 "dimensions (awh%d-dim%d and awh%d-dim%d). "
486 "If this is really what you want to do you will have to duplicate "
487 "this pull coordinate.",
488 awhBiasParams1.dimParams[d1].coordIndex + 1, k1 + 1, d1 + 1, k2 + 1,
490 gmx_fatal(FARGS, "%s", errormsg);
497 if (awhParams.shareBiasMultisim && !haveSharedBias)
500 "Sharing of biases over multiple simulations is requested, but no bias is marked "
501 "as shared (share-group > 0)");
504 /* mdrun does not support this (yet), but will check again */
505 if (haveBiasSharingWithinSimulation(awhParams))
508 "You have shared biases within a single simulation, but mdrun does not support "
513 AwhParams* readAndCheckAwhParams(std::vector<t_inpfile>* inp, const t_inputrec* ir, warninp_t wi)
515 AwhParams* awhParams;
519 /* Parameters common for all biases */
521 printStringNoNewline(inp, "The way to apply the biasing potential: convolved or umbrella");
522 opt = "awh-potential";
523 awhParams->ePotential = get_eeenum(inp, opt, eawhpotential_names, wi);
525 printStringNoNewline(inp,
526 "The random seed used for sampling the umbrella center in the case of "
527 "umbrella type potential");
529 awhParams->seed = get_eint(inp, opt, -1, wi);
530 if (awhParams->seed == -1)
532 awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
533 fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
536 printStringNoNewline(inp, "Data output interval in number of steps");
538 awhParams->nstOut = get_eint(inp, opt, 100000, wi);
539 if (awhParams->nstOut <= 0)
541 auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
542 opt.c_str(), awhParams->nstOut);
543 warning_error(wi, message);
545 /* This restriction can be removed by changing a flag of print_ebin() */
546 if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
548 auto message = formatString("%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(),
549 awhParams->nstOut, ir->nstenergy);
550 warning_error(wi, message);
553 printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
554 opt = "awh-nstsample";
555 awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
557 printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
558 opt = "awh-nsamples-update";
559 awhParams->numSamplesUpdateFreeEnergy = get_eint(inp, opt, 10, wi);
560 if (awhParams->numSamplesUpdateFreeEnergy <= 0)
562 warning_error(wi, opt + " needs to be an integer > 0");
565 printStringNoNewline(
566 inp, "When true, biases with share-group>0 are shared between multiple simulations");
567 opt = "awh-share-multisim";
568 awhParams->shareBiasMultisim = (get_eeenum(inp, opt, yesno_names, wi) != 0);
570 printStringNoNewline(inp, "The number of independent AWH biases");
572 awhParams->numBias = get_eint(inp, opt, 1, wi);
573 if (awhParams->numBias <= 0)
575 gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
578 /* Read the parameters specific to each AWH bias */
579 snew(awhParams->awhBiasParams, awhParams->numBias);
581 for (int k = 0; k < awhParams->numBias; k++)
583 bool bComment = (k == 0);
584 std::string prefixawh = formatString("awh%d", k + 1);
585 read_bias_params(inp, &awhParams->awhBiasParams[k], prefixawh, ir, wi, bComment);
588 /* Do a final consistency check before returning */
589 checkInputConsistencyAwh(*awhParams, wi);
591 if (ir->init_step != 0)
593 warning_error(wi, "With AWH init-step should be 0");
600 * Gets the period of a pull coordinate.
602 * \param[in] pullCoordParams The parameters for the pull coordinate.
603 * \param[in] pbc The PBC setup
604 * \param[in] intervalLength The length of the AWH interval for this pull coordinate
605 * \returns the period (or 0 if not periodic).
607 static double get_pull_coord_period(const t_pull_coord& pullCoordParams, const t_pbc& pbc, const real intervalLength)
611 if (pullCoordParams.eGeom == epullgDIR)
613 const real margin = 0.001;
614 // Make dims periodic when the interval covers > 95%
615 const real periodicFraction = 0.95;
617 // Check if the pull direction is along a box vector
618 for (int dim = 0; dim < pbc.ndim_ePBC; dim++)
620 const real boxLength = norm(pbc.box[dim]);
621 const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]);
622 if (innerProduct >= (1 - margin) * boxLength && innerProduct <= (1 + margin) * boxLength)
624 if (intervalLength > (1 + margin) * boxLength)
627 "The AWH interval (%f nm) for a pull coordinate is larger than the "
629 intervalLength, boxLength);
632 if (intervalLength > periodicFraction * boxLength)
639 else if (pullCoordParams.eGeom == epullgDIHEDRAL)
641 /* The dihedral angle is periodic in -180 to 180 deg */
649 * Checks if the given interval is defined in the correct periodic interval.
651 * \param[in] origin Start value of interval.
652 * \param[in] end End value of interval.
653 * \param[in] period Period (or 0 if not periodic).
654 * \returns true if the end point values are in the correct periodic interval.
656 static bool intervalIsInPeriodicInterval(double origin, double end, double period)
658 return (period == 0) || (std::fabs(origin) <= 0.5 * period && std::fabs(end) <= 0.5 * period);
662 * Checks if a value is within an interval.
664 * \param[in] origin Start value of interval.
665 * \param[in] end End value of interval.
666 * \param[in] period Period (or 0 if not periodic).
667 * \param[in] value Value to check.
668 * \returns true if the value is within the interval.
670 static bool valueIsInInterval(double origin, double end, double period, double value)
678 /* The interval closes within the periodic interval */
679 bIn_interval = (value >= origin) && (value <= end);
683 /* The interval wraps around the periodic boundary */
684 bIn_interval = ((value >= origin) && (value <= 0.5 * period))
685 || ((value >= -0.5 * period) && (value <= end));
690 bIn_interval = (value >= origin) && (value <= end);
697 * Check if the starting configuration is consistent with the given interval.
699 * \param[in] awhParams AWH parameters.
700 * \param[in,out] wi Struct for bookeeping warnings.
702 static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t wi)
704 for (int k = 0; k < awhParams->numBias; k++)
706 AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
707 for (int d = 0; d < awhBiasParams->ndim; d++)
709 AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
710 int coordIndex = dimParams->coordIndex;
711 double origin = dimParams->origin, end = dimParams->end, period = dimParams->period;
712 double coordValueInit = dimParams->coordValueInit;
714 if ((period == 0) && (origin > end))
717 "For the non-periodic pull coordinates awh%d-dim%d-start (%f) cannot be "
718 "larger than awh%d-dim%d-end (%f)",
719 k + 1, d + 1, origin, k + 1, d + 1, end);
722 /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
723 Make sure the AWH interval is within this reference interval.
725 Note: we could fairly simply allow using a more general interval (e.g. [x, x + period]) but it complicates
726 things slightly and I don't see that there is a great need for it. It would also mean that the interval would
727 depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180,
728 independent of AWH parameters.
730 if (!intervalIsInPeriodicInterval(origin, end, period))
733 "When using AWH with periodic pull coordinate geometries "
734 "awh%d-dim%d-start (%.8g) and "
735 "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take "
737 "minus half a period and plus half a period, i.e. in the interval [%.8g, "
739 k + 1, d + 1, origin, k + 1, d + 1, end, period, -0.5 * period, 0.5 * period);
742 /* Warn if the pull initial coordinate value is not in the grid */
743 if (!valueIsInInterval(origin, end, period, coordValueInit))
745 auto message = formatString(
746 "The initial coordinate value (%.8g) for pull coordinate index %d falls "
748 "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end "
750 "This can lead to large initial forces pulling the coordinate towards the "
751 "sampling interval.",
752 coordValueInit, coordIndex + 1, k + 1, d + 1, origin, k + 1, d + 1, end);
753 warning(wi, message);
759 void setStateDependentAwhParams(AwhParams* awhParams,
760 const pull_params_t* pull_params,
764 const tensor& compressibility,
765 const t_grpopts* inputrecGroupOptions,
768 /* The temperature is not really state depenendent but is not known
769 * when read_awhParams is called (in get ir).
770 * It is known first after do_index has been called in grompp.cpp.
772 if (inputrecGroupOptions->ref_t == nullptr || inputrecGroupOptions->ref_t[0] <= 0)
774 gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
776 for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
778 if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
781 "AWH biasing is currently only supported for identical temperatures for all "
782 "temperature coupling groups");
787 set_pbc(&pbc, ePBC, box);
789 for (int k = 0; k < awhParams->numBias; k++)
791 AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
792 for (int d = 0; d < awhBiasParams->ndim; d++)
794 AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
795 const t_pull_coord& pullCoordParams = pull_params->coord[dimParams->coordIndex];
797 if (pullCoordParams.eGeom == epullgDIRPBC)
800 "AWH does not support pull geometry '%s'. "
801 "If the maximum distance between the groups is always less than half the "
803 "you can use geometry '%s' instead.",
804 EPULLGEOM(epullgDIRPBC), EPULLGEOM(epullgDIR));
808 get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin);
809 // We would like to check for scaling, but we don't have the full inputrec available here
810 if (dimParams->period > 0
811 && !(pullCoordParams.eGeom == epullgANGLE || pullCoordParams.eGeom == epullgDIHEDRAL))
813 bool coordIsScaled = false;
814 for (int d2 = 0; d2 < DIM; d2++)
816 if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0)
818 coordIsScaled = true;
823 std::string mesg = gmx::formatString(
824 "AWH dimension %d of bias %d is periodic with pull geometry '%s', "
825 "while you should are applying pressure scaling to the corresponding "
826 "box vector, this is not supported.",
827 d + 1, k + 1, EPULLGEOM(pullCoordParams.eGeom));
828 warning(wi, mesg.c_str());
832 /* The initial coordinate value, converted to external user units. */
833 dimParams->coordValueInit = get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
835 dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(&pullCoordParams);
838 checkInputConsistencyInterval(awhParams, wi);
840 /* Register AWH as external potential with pull to check consistency. */
841 Awh::registerAwhWithPull(*awhParams, pull_work);