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,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.
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] = { "constant", "cutoff", "boltzmann",
66 "local-boltzmann", nullptr };
68 const char* eawhgrowth_names[eawhgrowthNR + 1] = { "exp-linear", "linear", nullptr };
70 const char* eawhpotential_names[eawhpotentialNR + 1] = { "convolved", "umbrella", nullptr };
72 const char* eawhcoordprovider_names[eawhcoordproviderNR + 1] = { "pull", nullptr };
75 * Read parameters of an AWH bias dimension.
77 * \param[in,out] inp Input file entries.
78 * \param[in] prefix Prefix for dimension parameters.
79 * \param[in,out] dimParams AWH dimensional parameters.
80 * \param[in] pull_params Pull parameters.
81 * \param[in,out] wi Struct for bookeeping warnings.
82 * \param[in] bComment True if comments should be printed.
84 static void readDimParams(std::vector<t_inpfile>* inp,
85 const std::string& prefix,
86 AwhDimParams* dimParams,
87 const pull_params_t* pull_params,
95 inp, "The provider of the reaction coordinate, currently only pull is supported");
98 opt = prefix + "-coord-provider";
99 dimParams->eCoordProvider = get_eeenum(inp, opt, eawhcoordprovider_names, wi);
103 printStringNoNewline(inp, "The coordinate index for this dimension");
105 opt = prefix + "-coord-index";
107 coordIndexInput = get_eint(inp, opt, 1, wi);
108 if (coordIndexInput < 1)
111 "Failed to read a valid coordinate index for %s. "
112 "Note that the pull coordinate indexing starts at 1.",
116 /* The pull coordinate indices start at 1 in the input file, at 0 internally */
117 dimParams->coordIndex = coordIndexInput - 1;
119 /* The pull settings need to be consistent with the AWH settings */
120 if (!(pull_params->coord[dimParams->coordIndex].eType == epullEXTERNAL))
122 gmx_fatal(FARGS, "AWH biasing can only be applied to pull type %s", EPULLTYPE(epullEXTERNAL));
125 if (dimParams->coordIndex >= pull_params->ncoord)
128 "The given AWH coordinate index (%d) is larger than the number of pull "
130 coordIndexInput, pull_params->ncoord);
132 if (pull_params->coord[dimParams->coordIndex].rate != 0)
134 auto message = formatString(
135 "Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate",
136 coordIndexInput, pull_params->coord[dimParams->coordIndex].rate);
137 warning_error(wi, message);
140 /* Grid params for each axis */
141 int eGeom = pull_params->coord[dimParams->coordIndex].eGeom;
145 printStringNoNewline(inp, "Start and end values for each coordinate dimension");
148 opt = prefix + "-start";
149 dimParams->origin = get_ereal(inp, opt, 0., wi);
151 opt = prefix + "-end";
152 dimParams->end = get_ereal(inp, opt, 0., wi);
154 if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
156 auto message = formatString(
157 "The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
158 "This will result in only one point along this axis in the coordinate value grid.",
159 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end);
160 warning(wi, message);
162 /* Check that the requested interval is in allowed range */
163 if (eGeom == epullgDIST)
165 if (dimParams->origin < 0 || dimParams->end < 0)
168 "%s-start (%g) or %s-end (%g) set to a negative value. With pull geometry "
169 "distance coordinate values are non-negative. "
170 "Perhaps you want to use geometry %s instead?",
171 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
172 EPULLGEOM(epullgDIR));
175 else if (eGeom == epullgANGLE || eGeom == epullgANGLEAXIS)
177 if (dimParams->origin < 0 || dimParams->end > 180)
180 "%s-start (%g) and %s-end (%g) are outside of the allowed range 0 to 180 deg "
181 "for pull geometries %s and %s ",
182 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
183 EPULLGEOM(epullgANGLE), EPULLGEOM(epullgANGLEAXIS));
186 else if (eGeom == epullgDIHEDRAL)
188 if (dimParams->origin < -180 || dimParams->end > 180)
191 "%s-start (%g) and %s-end (%g) are outside of the allowed range -180 to 180 "
192 "deg for pull geometry %s. ",
193 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
194 EPULLGEOM(epullgDIHEDRAL));
200 printStringNoNewline(
201 inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
203 opt = prefix + "-force-constant";
204 dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
205 if (dimParams->forceConstant <= 0)
207 warning_error(wi, "The force AWH bias force constant should be > 0");
212 printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps or rad^2/ps)");
214 opt = prefix + "-diffusion";
215 dimParams->diffusion = get_ereal(inp, opt, 0, wi);
217 if (dimParams->diffusion <= 0)
219 const double diffusion_default = 1e-5;
220 auto message = formatString(
221 "%s not explicitly set by user. You can choose to use a default "
222 "value (%g nm^2/ps or rad^2/ps) but this may very well be "
223 "non-optimal for your system!",
224 opt.c_str(), diffusion_default);
225 warning(wi, message);
226 dimParams->diffusion = diffusion_default;
231 printStringNoNewline(inp,
232 "Diameter that needs to be sampled around a point before it is "
233 "considered covered.");
235 opt = prefix + "-cover-diameter";
236 dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
238 if (dimParams->coverDiameter < 0)
240 gmx_fatal(FARGS, "%s (%g) cannot be negative.", opt.c_str(), dimParams->coverDiameter);
245 * Check consistency of input at the AWH bias level.
247 * \param[in] awhBiasParams AWH bias parameters.
248 * \param[in,out] wi Struct for bookkeeping warnings.
250 static void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
252 /* Covering diameter and sharing warning. */
253 for (int d = 0; d < awhBiasParams.ndim; d++)
255 double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
256 if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
259 "The covering diameter is only relevant to set for bias sharing simulations.");
265 * Read parameters of an AWH bias.
267 * \param[in,out] inp Input file entries.
268 * \param[in,out] awhBiasParams AWH dimensional parameters.
269 * \param[in] prefix Prefix for bias parameters.
270 * \param[in] ir Input parameter struct.
271 * \param[in,out] wi Struct for bookeeping warnings.
272 * \param[in] bComment True if comments should be printed.
274 static void read_bias_params(std::vector<t_inpfile>* inp,
275 AwhBiasParams* awhBiasParams,
276 const std::string& prefix,
277 const t_inputrec* ir,
283 printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
286 std::string opt = prefix + "-error-init";
287 /* We allow using a default value here without warning (but warn the user if the diffusion constant is not set). */
288 awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi);
289 if (awhBiasParams->errorInitial <= 0)
291 gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
296 printStringNoNewline(inp,
297 "Growth rate of the reference histogram determining the bias update "
298 "size: exp-linear or linear");
300 opt = prefix + "-growth";
301 awhBiasParams->eGrowth = get_eeenum(inp, opt, eawhgrowth_names, wi);
305 printStringNoNewline(inp,
306 "Start the simulation by equilibrating histogram towards the target "
307 "distribution: no or yes");
309 opt = prefix + "-equilibrate-histogram";
310 awhBiasParams->equilibrateHistogram = (get_eeenum(inp, opt, yesno_names, wi) != 0);
311 if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != eawhgrowthEXP_LINEAR)
314 formatString("Option %s will only have an effect for histogram growth type '%s'.",
315 opt.c_str(), EAWHGROWTH(eawhgrowthEXP_LINEAR));
316 warning(wi, message);
321 printStringNoNewline(
322 inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
324 opt = prefix + "-target";
325 awhBiasParams->eTarget = get_eeenum(inp, opt, eawhtarget_names, wi);
327 if ((awhBiasParams->eTarget == eawhtargetLOCALBOLTZMANN)
328 && (awhBiasParams->eGrowth == eawhgrowthEXP_LINEAR))
330 auto message = formatString(
331 "Target type '%s' combined with histogram growth type '%s' is not "
332 "expected to give stable bias updates. You probably want to use growth type "
334 EAWHTARGET(eawhtargetLOCALBOLTZMANN), EAWHGROWTH(eawhgrowthEXP_LINEAR),
335 EAWHGROWTH(eawhgrowthLINEAR));
336 warning(wi, message);
341 printStringNoNewline(inp,
342 "Boltzmann beta scaling factor for target distribution types "
343 "'boltzmann' and 'boltzmann-local'");
345 opt = prefix + "-target-beta-scaling";
346 awhBiasParams->targetBetaScaling = get_ereal(inp, opt, 0, wi);
348 switch (awhBiasParams->eTarget)
350 case eawhtargetBOLTZMANN:
351 case eawhtargetLOCALBOLTZMANN:
352 if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
354 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
355 awhBiasParams->targetBetaScaling, EAWHTARGET(awhBiasParams->eTarget));
359 if (awhBiasParams->targetBetaScaling != 0)
363 "Value for %s (%g) set explicitly but will not be used for target type %s.",
364 opt.c_str(), awhBiasParams->targetBetaScaling,
365 EAWHTARGET(awhBiasParams->eTarget));
372 printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
374 opt = prefix + "-target-cutoff";
375 awhBiasParams->targetCutoff = get_ereal(inp, opt, 0, wi);
377 switch (awhBiasParams->eTarget)
379 case eawhtargetCUTOFF:
380 if (awhBiasParams->targetCutoff <= 0)
382 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
383 awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
387 if (awhBiasParams->targetCutoff != 0)
391 "Value for %s (%g) set explicitly but will not be used for target type %s.",
392 opt.c_str(), awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
399 printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
401 opt = prefix + "-user-data";
402 awhBiasParams->bUserData = get_eeenum(inp, opt, yesno_names, wi);
406 printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
408 opt = prefix + "-share-group";
409 awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
410 if (awhBiasParams->shareGroup < 0)
412 warning_error(wi, "AWH bias share-group should be >= 0");
417 printStringNoNewline(inp, "Dimensionality of the coordinate");
419 opt = prefix + "-ndim";
420 awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
422 if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
424 gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(),
425 awhBiasParams->ndim, c_biasMaxNumDim);
427 if (awhBiasParams->ndim > 2)
430 "For awh-dim > 2 the estimate based on the diffusion and the initial error is "
431 "currently only a rough guideline."
432 " You should verify its usefulness for your system before production runs!");
434 snew(awhBiasParams->dimParams, awhBiasParams->ndim);
435 for (int d = 0; d < awhBiasParams->ndim; d++)
437 bComment = bComment && d == 0;
438 std::string prefixdim = prefix + formatString("-dim%d", d + 1);
439 readDimParams(inp, prefixdim, &awhBiasParams->dimParams[d], ir->pull, wi, bComment);
442 /* Check consistencies here that cannot be checked at read time at a lower level. */
443 checkInputConsistencyAwhBias(*awhBiasParams, wi);
447 * Check consistency of input at the AWH level.
449 * \param[in] awhParams AWH parameters.
450 * \param[in,out] wi Struct for bookkeeping warnings.
452 static void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi)
454 /* Each pull coord can map to at most 1 AWH coord.
455 * Check that we have a shared bias when requesting multisim sharing.
457 bool haveSharedBias = false;
458 for (int k1 = 0; k1 < awhParams.numBias; k1++)
460 const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1];
462 if (awhBiasParams1.shareGroup > 0)
464 haveSharedBias = true;
467 /* k1 is the reference AWH, k2 is the AWH we compare with (can be equal to k1) */
468 for (int k2 = k1; k2 < awhParams.numBias; k2++)
470 for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
472 const AwhBiasParams& awhBiasParams2 = awhParams.awhBiasParams[k2];
474 /* d1 is the reference dimension of the reference AWH. d2 is the dim index of the AWH to compare with. */
475 for (int d2 = 0; d2 < awhBiasParams2.ndim; d2++)
477 /* Give an error if (d1, k1) is different from (d2, k2) but the pull coordinate is the same */
478 if ((d1 != d2 || k1 != k2)
479 && (awhBiasParams1.dimParams[d1].coordIndex == awhBiasParams2.dimParams[d2].coordIndex))
481 char errormsg[STRLEN];
483 "One pull coordinate (%d) cannot be mapped to two separate AWH "
484 "dimensions (awh%d-dim%d and awh%d-dim%d). "
485 "If this is really what you want to do you will have to duplicate "
486 "this pull coordinate.",
487 awhBiasParams1.dimParams[d1].coordIndex + 1, k1 + 1, d1 + 1, k2 + 1,
489 gmx_fatal(FARGS, "%s", errormsg);
496 if (awhParams.shareBiasMultisim && !haveSharedBias)
499 "Sharing of biases over multiple simulations is requested, but no bias is marked "
500 "as shared (share-group > 0)");
503 /* mdrun does not support this (yet), but will check again */
504 if (haveBiasSharingWithinSimulation(awhParams))
507 "You have shared biases within a single simulation, but mdrun does not support "
512 AwhParams* readAndCheckAwhParams(std::vector<t_inpfile>* inp, const t_inputrec* ir, warninp_t wi)
514 AwhParams* awhParams;
518 /* Parameters common for all biases */
520 printStringNoNewline(inp, "The way to apply the biasing potential: convolved or umbrella");
521 opt = "awh-potential";
522 awhParams->ePotential = get_eeenum(inp, opt, eawhpotential_names, wi);
524 printStringNoNewline(inp,
525 "The random seed used for sampling the umbrella center in the case of "
526 "umbrella type potential");
528 awhParams->seed = get_eint(inp, opt, -1, wi);
529 if (awhParams->seed == -1)
531 awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
532 fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
535 printStringNoNewline(inp, "Data output interval in number of steps");
537 awhParams->nstOut = get_eint(inp, opt, 100000, wi);
538 if (awhParams->nstOut <= 0)
540 auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
541 opt.c_str(), awhParams->nstOut);
542 warning_error(wi, message);
544 /* This restriction can be removed by changing a flag of print_ebin() */
545 if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
547 auto message = formatString("%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(),
548 awhParams->nstOut, ir->nstenergy);
549 warning_error(wi, message);
552 printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
553 opt = "awh-nstsample";
554 awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
556 printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
557 opt = "awh-nsamples-update";
558 awhParams->numSamplesUpdateFreeEnergy = get_eint(inp, opt, 10, wi);
559 if (awhParams->numSamplesUpdateFreeEnergy <= 0)
561 warning_error(wi, opt + " needs to be an integer > 0");
564 printStringNoNewline(
565 inp, "When true, biases with share-group>0 are shared between multiple simulations");
566 opt = "awh-share-multisim";
567 awhParams->shareBiasMultisim = (get_eeenum(inp, opt, yesno_names, wi) != 0);
569 printStringNoNewline(inp, "The number of independent AWH biases");
571 awhParams->numBias = get_eint(inp, opt, 1, wi);
572 if (awhParams->numBias <= 0)
574 gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
577 /* Read the parameters specific to each AWH bias */
578 snew(awhParams->awhBiasParams, awhParams->numBias);
580 for (int k = 0; k < awhParams->numBias; k++)
582 bool bComment = (k == 0);
583 std::string prefixawh = formatString("awh%d", k + 1);
584 read_bias_params(inp, &awhParams->awhBiasParams[k], prefixawh, ir, wi, bComment);
587 /* Do a final consistency check before returning */
588 checkInputConsistencyAwh(*awhParams, wi);
590 if (ir->init_step != 0)
592 warning_error(wi, "With AWH init-step should be 0");
599 * Gets the period of a pull coordinate.
601 * \param[in] pullCoordParams The parameters for the pull coordinate.
602 * \param[in] pbc The PBC setup
603 * \param[in] intervalLength The length of the AWH interval for this pull coordinate
604 * \returns the period (or 0 if not periodic).
606 static double get_pull_coord_period(const t_pull_coord& pullCoordParams, const t_pbc& pbc, const real intervalLength)
610 if (pullCoordParams.eGeom == epullgDIR)
612 const real margin = 0.001;
613 // Make dims periodic when the interval covers > 95%
614 const real periodicFraction = 0.95;
616 // Check if the pull direction is along a box vector
617 for (int dim = 0; dim < pbc.ndim_ePBC; dim++)
619 const real boxLength = norm(pbc.box[dim]);
620 const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]);
621 if (innerProduct >= (1 - margin) * boxLength && innerProduct <= (1 + margin) * boxLength)
623 GMX_RELEASE_ASSERT(intervalLength < (1 + margin) * boxLength,
624 "We have checked before that interval <= period");
625 if (intervalLength > periodicFraction * boxLength)
632 else if (pullCoordParams.eGeom == epullgDIHEDRAL)
634 /* The dihedral angle is periodic in -180 to 180 deg */
642 * Checks if the given interval is defined in the correct periodic interval.
644 * \param[in] origin Start value of interval.
645 * \param[in] end End value of interval.
646 * \param[in] period Period (or 0 if not periodic).
647 * \returns true if the end point values are in the correct periodic interval.
649 static bool intervalIsInPeriodicInterval(double origin, double end, double period)
651 return (period == 0) || (std::fabs(origin) <= 0.5 * period && std::fabs(end) <= 0.5 * period);
655 * Checks if a value is within an interval.
657 * \param[in] origin Start value of interval.
658 * \param[in] end End value of interval.
659 * \param[in] period Period (or 0 if not periodic).
660 * \param[in] value Value to check.
661 * \returns true if the value is within the interval.
663 static bool valueIsInInterval(double origin, double end, double period, double value)
671 /* The interval closes within the periodic interval */
672 bIn_interval = (value >= origin) && (value <= end);
676 /* The interval wraps around the periodic boundary */
677 bIn_interval = ((value >= origin) && (value <= 0.5 * period))
678 || ((value >= -0.5 * period) && (value <= end));
683 bIn_interval = (value >= origin) && (value <= end);
690 * Check if the starting configuration is consistent with the given interval.
692 * \param[in] awhParams AWH parameters.
693 * \param[in,out] wi Struct for bookeeping warnings.
695 static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t wi)
697 for (int k = 0; k < awhParams->numBias; k++)
699 AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
700 for (int d = 0; d < awhBiasParams->ndim; d++)
702 AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
703 int coordIndex = dimParams->coordIndex;
704 double origin = dimParams->origin, end = dimParams->end, period = dimParams->period;
705 double coordValueInit = dimParams->coordValueInit;
707 if ((period == 0) && (origin > end))
710 "For the non-periodic pull coordinates awh%d-dim%d-start (%f) cannot be "
711 "larger than awh%d-dim%d-end (%f)",
712 k + 1, d + 1, origin, k + 1, d + 1, end);
715 /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
716 Make sure the AWH interval is within this reference interval.
718 Note: we could fairly simply allow using a more general interval (e.g. [x, x + period]) but it complicates
719 things slightly and I don't see that there is a great need for it. It would also mean that the interval would
720 depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180,
721 independent of AWH parameters.
723 if (!intervalIsInPeriodicInterval(origin, end, period))
726 "When using AWH with periodic pull coordinate geometries "
727 "awh%d-dim%d-start (%.8g) and "
728 "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take "
730 "minus half a period and plus half a period, i.e. in the interval [%.8g, "
732 k + 1, d + 1, origin, k + 1, d + 1, end, period, -0.5 * period, 0.5 * period);
735 /* Warn if the pull initial coordinate value is not in the grid */
736 if (!valueIsInInterval(origin, end, period, coordValueInit))
738 auto message = formatString(
739 "The initial coordinate value (%.8g) for pull coordinate index %d falls "
741 "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end "
743 "This can lead to large initial forces pulling the coordinate towards the "
744 "sampling interval.",
745 coordValueInit, coordIndex + 1, k + 1, d + 1, origin, k + 1, d + 1, end);
746 warning(wi, message);
752 void setStateDependentAwhParams(AwhParams* awhParams,
753 const pull_params_t* pull_params,
757 const tensor& compressibility,
758 const t_grpopts* inputrecGroupOptions,
761 /* The temperature is not really state depenendent but is not known
762 * when read_awhParams is called (in get ir).
763 * It is known first after do_index has been called in grompp.cpp.
765 if (inputrecGroupOptions->ref_t == nullptr || inputrecGroupOptions->ref_t[0] <= 0)
767 gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
769 for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
771 if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
774 "AWH biasing is currently only supported for identical temperatures for all "
775 "temperature coupling groups");
780 set_pbc(&pbc, ePBC, box);
782 for (int k = 0; k < awhParams->numBias; k++)
784 AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
785 for (int d = 0; d < awhBiasParams->ndim; d++)
787 AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
788 const t_pull_coord& pullCoordParams = pull_params->coord[dimParams->coordIndex];
790 if (pullCoordParams.eGeom == epullgDIRPBC)
793 "AWH does not support pull geometry '%s'. "
794 "If the maximum distance between the groups is always less than half the "
796 "you can use geometry '%s' instead.",
797 EPULLGEOM(epullgDIRPBC), EPULLGEOM(epullgDIR));
801 get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin);
802 // We would like to check for scaling, but we don't have the full inputrec available here
803 if (dimParams->period > 0
804 && !(pullCoordParams.eGeom == epullgANGLE || pullCoordParams.eGeom == epullgDIHEDRAL))
806 bool coordIsScaled = false;
807 for (int d2 = 0; d2 < DIM; d2++)
809 if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0)
811 coordIsScaled = true;
816 std::string mesg = gmx::formatString(
817 "AWH dimension %d of bias %d is periodic with pull geometry '%s', "
818 "while you should are applying pressure scaling to the corresponding "
819 "box vector, this is not supported.",
820 d + 1, k + 1, EPULLGEOM(pullCoordParams.eGeom));
821 warning(wi, mesg.c_str());
825 /* The initial coordinate value, converted to external user units. */
826 dimParams->coordValueInit = get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
828 dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(&pullCoordParams);
831 checkInputConsistencyInterval(awhParams, wi);
833 /* Register AWH as external potential with pull to check consistency. */
834 Awh::registerAwhWithPull(*awhParams, pull_work);