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 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, 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.
11 *
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.
16 *
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.
21 *
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.
26 *
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.
34 *
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.
37 */
38 #include "gmxpre.h"
39
40 #include "read_params.h"
41
42 #include "gromacs/applied_forces/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/multipletimestepping.h"
52 #include "gromacs/mdtypes/pull_params.h"
53 #include "gromacs/pbcutil/pbc.h"
54 #include "gromacs/pulling/pull.h"
55 #include "gromacs/random/seed.h"
56 #include "gromacs/topology/mtop_util.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/smalloc.h"
60 #include "gromacs/utility/stringutil.h"
61
62 #include "biasparams.h"
63 #include "biassharing.h"
64
65 namespace gmx
66 {
67
68 const char* eawhtarget_names[eawhtargetNR + 1] = { "constant", "cutoff", "boltzmann",
69 "local-boltzmann", nullptr };
70
71 const char* eawhgrowth_names[eawhgrowthNR + 1] = { "exp-linear", "linear", nullptr };
72
73 const char* eawhpotential_names[eawhpotentialNR + 1] = { "convolved", "umbrella", nullptr };
74
75 const char* eawhcoordprovider_names[eawhcoordproviderNR + 1] = { "pull", "fep-lambda", nullptr };
76
77 namespace
78 {
79
80 /*! \brief
81 * Check multiple time-stepping consistency between AWH and pull and/or FEP
82 *
83 * \param[in] inputrec Inputput parameter struct.
84 * \param[in,out] wi Struct for bookeeping warnings.
85 */
checkMtsConsistency(const t_inputrec & inputrec,warninp_t wi)86 void checkMtsConsistency(const t_inputrec& inputrec, warninp_t wi)
87 {
88 if (!inputrec.useMts)
89 {
90 return;
91 }
92
93 GMX_RELEASE_ASSERT(inputrec.mtsLevels.size() == 2, "Only 2 MTS levels supported here");
94
95 bool usesPull = false;
96 bool usesFep = false;
97 for (int b = 0; b < inputrec.awhParams->numBias; b++)
98 {
99 const auto& biasParams = inputrec.awhParams->awhBiasParams[b];
100 for (int d = 0; d < biasParams.ndim; d++)
101 {
102 switch (biasParams.dimParams[d].eCoordProvider)
103 {
104 case eawhcoordproviderPULL: usesPull = true; break;
105 case eawhcoordproviderFREE_ENERGY_LAMBDA: usesFep = true; break;
106 default: GMX_RELEASE_ASSERT(false, "Unsupported coord provider");
107 }
108 }
109 }
110 const int awhMtsLevel = forceGroupMtsLevel(inputrec.mtsLevels, MtsForceGroups::Awh);
111 if (usesPull && forceGroupMtsLevel(inputrec.mtsLevels, MtsForceGroups::Pull) != awhMtsLevel)
112 {
113 warning_error(wi,
114 "When AWH is applied to pull coordinates, pull and AWH should be computed at "
115 "the same MTS level");
116 }
117 if (usesFep && awhMtsLevel != ssize(inputrec.mtsLevels) - 1)
118 {
119 warning_error(wi,
120 "When AWH is applied to the free-energy lambda with MTS, AWH should be "
121 "computed at the slow MTS level");
122 }
123
124 if (inputrec.awhParams->nstSampleCoord % inputrec.mtsLevels[awhMtsLevel].stepFactor != 0)
125 {
126 warning_error(wi,
127 "With MTS applied to AWH, awh-nstsample should be a multiple of mts-factor");
128 }
129 }
130
131 /*! \brief
132 * Check the parameters of an AWH bias pull dimension.
133 *
134 * \param[in] prefix Prefix for dimension parameters.
135 * \param[in,out] dimParams AWH dimensional parameters.
136 * \param[in] pull_params Pull parameters.
137 * \param[in,out] wi Struct for bookeeping warnings.
138 */
checkPullDimParams(const std::string & prefix,AwhDimParams * dimParams,const pull_params_t & pull_params,warninp_t wi)139 void checkPullDimParams(const std::string& prefix,
140 AwhDimParams* dimParams,
141 const pull_params_t& pull_params,
142 warninp_t wi)
143 {
144 if (dimParams->coordIndex < 0)
145 {
146 gmx_fatal(FARGS,
147 "Failed to read a valid coordinate index for %s-coord-index. "
148 "Note that the pull coordinate indexing starts at 1.",
149 prefix.c_str());
150 }
151 if (dimParams->coordIndex >= pull_params.ncoord)
152 {
153 gmx_fatal(FARGS,
154 "The given AWH coordinate index (%d) is larger than the number of pull "
155 "coordinates (%d)",
156 dimParams->coordIndex + 1, pull_params.ncoord);
157 }
158 if (pull_params.coord[dimParams->coordIndex].rate != 0)
159 {
160 auto message = formatString(
161 "Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate",
162 dimParams->coordIndex + 1, pull_params.coord[dimParams->coordIndex].rate);
163 warning_error(wi, message);
164 }
165
166 if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
167 {
168 auto message = formatString(
169 "The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
170 "This will result in only one point along this axis in the coordinate value grid.",
171 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end);
172 warning(wi, message);
173 }
174
175 if (dimParams->forceConstant <= 0)
176 {
177 warning_error(wi, "The force AWH bias force constant should be > 0");
178 }
179
180 /* Grid params for each axis */
181 int eGeom = pull_params.coord[dimParams->coordIndex].eGeom;
182
183 /* Check that the requested interval is in allowed range */
184 if (eGeom == epullgDIST)
185 {
186 if (dimParams->origin < 0 || dimParams->end < 0)
187 {
188 gmx_fatal(FARGS,
189 "%s-start (%g) or %s-end (%g) set to a negative value. With pull "
190 "geometry distance coordinate values are non-negative. "
191 "Perhaps you want to use geometry %s instead?",
192 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
193 EPULLGEOM(epullgDIR));
194 }
195 }
196 else if (eGeom == epullgANGLE || eGeom == epullgANGLEAXIS)
197 {
198 if (dimParams->origin < 0 || dimParams->end > 180)
199 {
200 gmx_fatal(FARGS,
201 "%s-start (%g) and %s-end (%g) are outside of the allowed range "
202 "0 to 180 deg for pull geometries %s and %s ",
203 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
204 EPULLGEOM(epullgANGLE), EPULLGEOM(epullgANGLEAXIS));
205 }
206 }
207 else if (eGeom == epullgDIHEDRAL)
208 {
209 if (dimParams->origin < -180 || dimParams->end > 180)
210 {
211 gmx_fatal(FARGS,
212 "%s-start (%g) and %s-end (%g) are outside of the allowed range "
213 "-180 to 180 deg for pull geometry %s. ",
214 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end,
215 EPULLGEOM(epullgDIHEDRAL));
216 }
217 }
218 }
219
220 /*! \brief
221 * Check parameters of an AWH bias in a free energy lambda state dimension.
222 *
223 * \param[in] prefix Prefix for dimension parameters.
224 * \param[in,out] dimParams AWH dimensional parameters.
225 * \param[in] lambdaParams The free energy lambda related parameters.
226 * \param[in] efep This is the type of FEP calculation (efep enumerator).
227 * \param[in,out] wi Struct for bookeeping warnings.
228 */
checkFepLambdaDimParams(const std::string & prefix,const AwhDimParams * dimParams,const t_lambda * lambdaParams,const int efep,warninp_t wi)229 void checkFepLambdaDimParams(const std::string& prefix,
230 const AwhDimParams* dimParams,
231 const t_lambda* lambdaParams,
232 const int efep,
233 warninp_t wi)
234 {
235 std::string opt;
236
237 if (!lambdaParams)
238 {
239 gmx_fatal(FARGS,
240 "There must be free energy input if using AWH to steer the free energy lambda "
241 "state.");
242 }
243
244 if (lambdaParams->lambda_neighbors != -1)
245 {
246 gmx_fatal(FARGS,
247 "When running AWH coupled to the free energy lambda state all lambda states "
248 "should be used as neighbors in order to get correct probabilities, i.e. "
249 "calc-lambda-neighbors (%d) must be %d.",
250 lambdaParams->lambda_neighbors, -1);
251 }
252
253 if (efep == efepSLOWGROWTH || lambdaParams->delta_lambda != 0)
254 {
255 gmx_fatal(FARGS,
256 "AWH coupled to the free energy lambda state is not compatible with slow-growth "
257 "and delta-lambda must be 0.");
258 }
259
260 if (efep == efepEXPANDED)
261 {
262 gmx_fatal(FARGS,
263 "AWH is not treated like other expanded ensemble methods. Do not use expanded.");
264 }
265
266 if (dimParams->origin < 0)
267 {
268 opt = prefix + "-start";
269 gmx_fatal(FARGS,
270 "When running AWH coupled to the free energy lambda state the lower lambda state "
271 "for AWH, %s (%.0f), must be >= 0.",
272 opt.c_str(), dimParams->origin);
273 }
274 if (dimParams->end >= lambdaParams->n_lambda)
275 {
276 opt = prefix + "-end";
277 gmx_fatal(FARGS,
278 "When running AWH coupled to the free energy lambda state the upper lambda state "
279 "for AWH, %s (%.0f), must be < n_lambda (%d).",
280 opt.c_str(), dimParams->origin, lambdaParams->n_lambda);
281 }
282 if (gmx_within_tol(dimParams->end - dimParams->origin, 0, GMX_REAL_EPS))
283 {
284 auto message = formatString(
285 "The given interval length given by %s-start (%g) and %s-end (%g) is zero. "
286 "This will result in only one lambda point along this free energy lambda state "
287 "axis in the coordinate value grid.",
288 prefix.c_str(), dimParams->origin, prefix.c_str(), dimParams->end);
289 warning(wi, message);
290 }
291
292 if (dimParams->forceConstant != 0)
293 {
294 warning_error(
295 wi,
296 "The force AWH bias force constant is not used with free energy lambda state as "
297 "coordinate provider.");
298 }
299 }
300
301 /*! \brief
302 * Check that AWH FEP is not combined with incompatible decoupling.
303 *
304 * \param[in] mtop System topology.
305 * \param[in,out] wi Struct for bookeeping warnings.
306 */
checkFepLambdaDimDecouplingConsistency(const gmx_mtop_t & mtop,warninp_t wi)307 void checkFepLambdaDimDecouplingConsistency(const gmx_mtop_t& mtop, warninp_t wi)
308 {
309 if (haveFepPerturbedMasses(mtop))
310 {
311 warning(wi,
312 "Masses may not be perturbed when using the free energy lambda state as AWH "
313 "coordinate provider. If you are using fep-lambdas to specify lambda states "
314 "make sure that you also specify mass-lambdas without perturbation.");
315 }
316 if (havePerturbedConstraints(mtop))
317 {
318 warning(wi,
319 "Constraints may not be perturbed when using the free energy lambda state as AWH "
320 "coordinate provider. If you are using fep-lambdas to specify lambda states "
321 "make sure that you also specify mass-lambdas without perturbation.");
322 }
323 }
324
325 /*! \brief
326 * Read parameters of an AWH bias dimension.
327 *
328 * \param[in,out] inp Input file entries.
329 * \param[in] prefix Prefix for dimension parameters.
330 * \param[in,out] dimParams AWH dimensional parameters.
331 * \param[in,out] wi Struct for bookeeping warnings.
332 * \param[in] bComment True if comments should be printed.
333 */
readDimParams(std::vector<t_inpfile> * inp,const std::string & prefix,AwhDimParams * dimParams,warninp_t wi,bool bComment)334 void readDimParams(std::vector<t_inpfile>* inp,
335 const std::string& prefix,
336 AwhDimParams* dimParams,
337 warninp_t wi,
338 bool bComment)
339 {
340 std::string opt;
341 if (bComment)
342 {
343 printStringNoNewline(
344 inp,
345 "The provider of the reaction coordinate, "
346 "currently only 'pull' and 'fep-lambda' (free energy lambda state) is supported");
347 }
348 opt = prefix + "-coord-provider";
349 dimParams->eCoordProvider = get_eeenum(inp, opt, eawhcoordprovider_names, wi);
350
351 if (bComment)
352 {
353 printStringNoNewline(inp, "The coordinate index for this dimension");
354 }
355 opt = prefix + "-coord-index";
356 int coordIndexInput;
357 coordIndexInput = get_eint(inp, opt, 1, wi);
358
359 /* The pull coordinate indices start at 1 in the input file, at 0 internally */
360 dimParams->coordIndex = coordIndexInput - 1;
361
362 if (bComment)
363 {
364 printStringNoNewline(inp, "Start and end values for each coordinate dimension");
365 }
366 opt = prefix + "-start";
367 dimParams->origin = get_ereal(inp, opt, 0., wi);
368 opt = prefix + "-end";
369 dimParams->end = get_ereal(inp, opt, 0., wi);
370
371 if (bComment)
372 {
373 printStringNoNewline(
374 inp, "The force constant for this coordinate (kJ/mol/nm^2 or kJ/mol/rad^2)");
375 }
376 opt = prefix + "-force-constant";
377 dimParams->forceConstant = get_ereal(inp, opt, 0, wi);
378
379 if (bComment)
380 {
381 printStringNoNewline(inp, "Estimated diffusion constant (nm^2/ps, rad^2/ps or ps^-1)");
382 }
383 opt = prefix + "-diffusion";
384 dimParams->diffusion = get_ereal(inp, opt, 0, wi);
385
386 if (bComment)
387 {
388 printStringNoNewline(inp,
389 "Diameter that needs to be sampled around a point before it is "
390 "considered covered. In FEP dimensions the cover diameter is "
391 "specified in lambda states.");
392 }
393 opt = prefix + "-cover-diameter";
394 dimParams->coverDiameter = get_ereal(inp, opt, 0, wi);
395 }
396
397 /*! \brief
398 * Check the parameters of an AWH bias dimension.
399 *
400 * \param[in] prefix Prefix for dimension parameters.
401 * \param[in,out] dimParams AWH dimensional parameters.
402 * \param[in] ir Input parameter struct.
403 * \param[in,out] wi Struct for bookeeping warnings.
404 */
checkDimParams(const std::string & prefix,AwhDimParams * dimParams,const t_inputrec * ir,warninp_t wi)405 void checkDimParams(const std::string& prefix, AwhDimParams* dimParams, const t_inputrec* ir, warninp_t wi)
406 {
407 if (dimParams->eCoordProvider == eawhcoordproviderPULL)
408 {
409 if (!ir->bPull)
410 {
411 gmx_fatal(FARGS,
412 "AWH biasing along a pull dimension is only compatible with COM pulling "
413 "turned on");
414 }
415 checkPullDimParams(prefix, dimParams, *ir->pull, wi);
416 }
417 else if (dimParams->eCoordProvider == eawhcoordproviderFREE_ENERGY_LAMBDA)
418 {
419 if (ir->efep == efepNO)
420 {
421 gmx_fatal(FARGS,
422 "AWH biasing along a free energy lambda state dimension is only compatible "
423 "with free energy turned on");
424 }
425 checkFepLambdaDimParams(prefix, dimParams, ir->fepvals, ir->efep, wi);
426 }
427 else
428 {
429 gmx_fatal(FARGS,
430 "AWH biasing can only be applied to pull and free energy lambda state "
431 "coordinates");
432 }
433 }
434
435 /*! \brief
436 * Check consistency of input at the AWH bias level.
437 *
438 * \param[in] awhBiasParams AWH bias parameters.
439 * \param[in,out] wi Struct for bookkeeping warnings.
440 */
checkInputConsistencyAwhBias(const AwhBiasParams & awhBiasParams,warninp_t wi)441 void checkInputConsistencyAwhBias(const AwhBiasParams& awhBiasParams, warninp_t wi)
442 {
443 /* Covering diameter and sharing warning. */
444 for (int d = 0; d < awhBiasParams.ndim; d++)
445 {
446 double coverDiameter = awhBiasParams.dimParams[d].coverDiameter;
447 if (awhBiasParams.shareGroup <= 0 && coverDiameter > 0)
448 {
449 warning(wi,
450 "The covering diameter is only relevant to set for bias sharing simulations.");
451 }
452 }
453 }
454
455 /*! \brief
456 * Read parameters of an AWH bias.
457 *
458 * \param[in,out] inp Input file entries.
459 * \param[in,out] awhBiasParams AWH dimensional parameters.
460 * \param[in] prefix Prefix for bias parameters.
461 * \param[in,out] wi Struct for bookeeping warnings.
462 * \param[in] bComment True if comments should be printed.
463 */
readBiasParams(std::vector<t_inpfile> * inp,AwhBiasParams * awhBiasParams,const std::string & prefix,warninp_t wi,bool bComment)464 void readBiasParams(std::vector<t_inpfile>* inp,
465 AwhBiasParams* awhBiasParams,
466 const std::string& prefix,
467 warninp_t wi,
468 bool bComment)
469 {
470 if (bComment)
471 {
472 printStringNoNewline(inp, "Estimated initial PMF error (kJ/mol)");
473 }
474
475 std::string opt = prefix + "-error-init";
476 awhBiasParams->errorInitial = get_ereal(inp, opt, 10, wi);
477
478 if (bComment)
479 {
480 printStringNoNewline(inp,
481 "Growth rate of the reference histogram determining the bias update "
482 "size: exp-linear or linear");
483 }
484 opt = prefix + "-growth";
485 awhBiasParams->eGrowth = get_eeenum(inp, opt, eawhgrowth_names, wi);
486
487 if (bComment)
488 {
489 printStringNoNewline(inp,
490 "Start the simulation by equilibrating histogram towards the target "
491 "distribution: no or yes");
492 }
493 opt = prefix + "-equilibrate-histogram";
494 awhBiasParams->equilibrateHistogram = (get_eeenum(inp, opt, yesno_names, wi) != 0);
495
496 if (bComment)
497 {
498 printStringNoNewline(
499 inp, "Target distribution type: constant, cutoff, boltzmann or local-boltzmann");
500 }
501 opt = prefix + "-target";
502 awhBiasParams->eTarget = get_eeenum(inp, opt, eawhtarget_names, wi);
503
504 if (bComment)
505 {
506 printStringNoNewline(inp,
507 "Boltzmann beta scaling factor for target distribution types "
508 "'boltzmann' and 'boltzmann-local'");
509 }
510 opt = prefix + "-target-beta-scaling";
511 awhBiasParams->targetBetaScaling = get_ereal(inp, opt, 0, wi);
512
513 if (bComment)
514 {
515 printStringNoNewline(inp, "Free energy cutoff value for target distribution type 'cutoff'");
516 }
517 opt = prefix + "-target-cutoff";
518 awhBiasParams->targetCutoff = get_ereal(inp, opt, 0, wi);
519
520 if (bComment)
521 {
522 printStringNoNewline(inp, "Initialize PMF and target with user data: no or yes");
523 }
524 opt = prefix + "-user-data";
525 awhBiasParams->bUserData = get_eeenum(inp, opt, yesno_names, wi);
526
527 if (bComment)
528 {
529 printStringNoNewline(inp, "Group index to share the bias with, 0 means not shared");
530 }
531 opt = prefix + "-share-group";
532 awhBiasParams->shareGroup = get_eint(inp, opt, 0, wi);
533
534 if (bComment)
535 {
536 printStringNoNewline(inp, "Dimensionality of the coordinate");
537 }
538 opt = prefix + "-ndim";
539 awhBiasParams->ndim = get_eint(inp, opt, 0, wi);
540
541 /* Check this before starting to read the AWH dimension parameters. */
542 if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
543 {
544 gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(),
545 awhBiasParams->ndim, c_biasMaxNumDim);
546 }
547 snew(awhBiasParams->dimParams, awhBiasParams->ndim);
548 for (int d = 0; d < awhBiasParams->ndim; d++)
549 {
550 bComment = bComment && d == 0;
551 std::string prefixdim = prefix + formatString("-dim%d", d + 1);
552 readDimParams(inp, prefixdim, &awhBiasParams->dimParams[d], wi, bComment);
553 }
554 }
555
556 /*! \brief
557 * Check the parameters of an AWH bias.
558 *
559 * \param[in] awhBiasParams AWH dimensional parameters.
560 * \param[in] prefix Prefix for bias parameters.
561 * \param[in] ir Input parameter struct.
562 * \param[in,out] wi Struct for bookeeping warnings.
563 */
checkBiasParams(const AwhBiasParams * awhBiasParams,const std::string & prefix,const t_inputrec * ir,warninp_t wi)564 void checkBiasParams(const AwhBiasParams* awhBiasParams, const std::string& prefix, const t_inputrec* ir, warninp_t wi)
565 {
566 std::string opt = prefix + "-error-init";
567 if (awhBiasParams->errorInitial <= 0)
568 {
569 gmx_fatal(FARGS, "%s needs to be > 0.", opt.c_str());
570 }
571
572 opt = prefix + "-equilibrate-histogram";
573 if (awhBiasParams->equilibrateHistogram && awhBiasParams->eGrowth != eawhgrowthEXP_LINEAR)
574 {
575 auto message =
576 formatString("Option %s will only have an effect for histogram growth type '%s'.",
577 opt.c_str(), EAWHGROWTH(eawhgrowthEXP_LINEAR));
578 warning(wi, message);
579 }
580
581 if ((awhBiasParams->eTarget == eawhtargetLOCALBOLTZMANN)
582 && (awhBiasParams->eGrowth == eawhgrowthEXP_LINEAR))
583 {
584 auto message = formatString(
585 "Target type '%s' combined with histogram growth type '%s' is not "
586 "expected to give stable bias updates. You probably want to use growth type "
587 "'%s' instead.",
588 EAWHTARGET(eawhtargetLOCALBOLTZMANN), EAWHGROWTH(eawhgrowthEXP_LINEAR),
589 EAWHGROWTH(eawhgrowthLINEAR));
590 warning(wi, message);
591 }
592
593 opt = prefix + "-target-beta-scaling";
594 switch (awhBiasParams->eTarget)
595 {
596 case eawhtargetBOLTZMANN:
597 case eawhtargetLOCALBOLTZMANN:
598 if (awhBiasParams->targetBetaScaling < 0 || awhBiasParams->targetBetaScaling > 1)
599 {
600 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
601 awhBiasParams->targetBetaScaling, EAWHTARGET(awhBiasParams->eTarget));
602 }
603 break;
604 default:
605 if (awhBiasParams->targetBetaScaling != 0)
606 {
607 gmx_fatal(
608 FARGS,
609 "Value for %s (%g) set explicitly but will not be used for target type %s.",
610 opt.c_str(), awhBiasParams->targetBetaScaling,
611 EAWHTARGET(awhBiasParams->eTarget));
612 }
613 break;
614 }
615
616 opt = prefix + "-target-cutoff";
617 switch (awhBiasParams->eTarget)
618 {
619 case eawhtargetCUTOFF:
620 if (awhBiasParams->targetCutoff <= 0)
621 {
622 gmx_fatal(FARGS, "%s = %g is not useful for target type %s.", opt.c_str(),
623 awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
624 }
625 break;
626 default:
627 if (awhBiasParams->targetCutoff != 0)
628 {
629 gmx_fatal(
630 FARGS,
631 "Value for %s (%g) set explicitly but will not be used for target type %s.",
632 opt.c_str(), awhBiasParams->targetCutoff, EAWHTARGET(awhBiasParams->eTarget));
633 }
634 break;
635 }
636
637 opt = prefix + "-share-group";
638 if (awhBiasParams->shareGroup < 0)
639 {
640 warning_error(wi, "AWH bias share-group should be >= 0");
641 }
642
643 opt = prefix + "-ndim";
644 if (awhBiasParams->ndim <= 0 || awhBiasParams->ndim > c_biasMaxNumDim)
645 {
646 gmx_fatal(FARGS, "%s (%d) needs to be > 0 and at most %d\n", opt.c_str(),
647 awhBiasParams->ndim, c_biasMaxNumDim);
648 }
649 if (awhBiasParams->ndim > 2)
650 {
651 warning_note(wi,
652 "For awh-dim > 2 the estimate based on the diffusion and the initial error is "
653 "currently only a rough guideline."
654 " You should verify its usefulness for your system before production runs!");
655 }
656 for (int d = 0; d < awhBiasParams->ndim; d++)
657 {
658 std::string prefixdim = prefix + formatString("-dim%d", d + 1);
659 checkDimParams(prefixdim, &awhBiasParams->dimParams[d], ir, wi);
660 }
661
662 /* Check consistencies here that cannot be checked at read time at a lower level. */
663 checkInputConsistencyAwhBias(*awhBiasParams, wi);
664 }
665
666 /*! \brief
667 * Check consistency of input at the AWH level.
668 *
669 * \param[in] awhParams AWH parameters.
670 * \param[in,out] wi Struct for bookkeeping warnings.
671 */
checkInputConsistencyAwh(const AwhParams & awhParams,warninp_t wi)672 void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi)
673 {
674 /* Each pull coord can map to at most 1 AWH coord.
675 * Check that we have a shared bias when requesting multisim sharing.
676 */
677 bool haveSharedBias = false;
678 for (int k1 = 0; k1 < awhParams.numBias; k1++)
679 {
680 const AwhBiasParams& awhBiasParams1 = awhParams.awhBiasParams[k1];
681
682 if (awhBiasParams1.shareGroup > 0)
683 {
684 haveSharedBias = true;
685 }
686
687 /* k1 is the reference AWH, k2 is the AWH we compare with (can be equal to k1) */
688 for (int k2 = k1; k2 < awhParams.numBias; k2++)
689 {
690 for (int d1 = 0; d1 < awhBiasParams1.ndim; d1++)
691 {
692 if (awhBiasParams1.dimParams[d1].eCoordProvider == eawhcoordproviderFREE_ENERGY_LAMBDA)
693 {
694 continue;
695 }
696 const AwhBiasParams& awhBiasParams2 = awhParams.awhBiasParams[k2];
697
698 /* d1 is the reference dimension of the reference AWH. d2 is the dim index of the AWH to compare with. */
699 for (int d2 = 0; d2 < awhBiasParams2.ndim; d2++)
700 {
701 if (awhBiasParams2.dimParams[d2].eCoordProvider == eawhcoordproviderFREE_ENERGY_LAMBDA)
702 {
703 continue;
704 }
705 /* Give an error if (d1, k1) is different from (d2, k2) but the pull coordinate is the same */
706 if ((d1 != d2 || k1 != k2)
707 && (awhBiasParams1.dimParams[d1].coordIndex == awhBiasParams2.dimParams[d2].coordIndex))
708 {
709 char errormsg[STRLEN];
710 sprintf(errormsg,
711 "One pull coordinate (%d) cannot be mapped to two separate AWH "
712 "dimensions (awh%d-dim%d and awh%d-dim%d). "
713 "If this is really what you want to do you will have to duplicate "
714 "this pull coordinate.",
715 awhBiasParams1.dimParams[d1].coordIndex + 1, k1 + 1, d1 + 1, k2 + 1,
716 d2 + 1);
717 gmx_fatal(FARGS, "%s", errormsg);
718 }
719 }
720 }
721 }
722 }
723
724 if (awhParams.shareBiasMultisim && !haveSharedBias)
725 {
726 warning(wi,
727 "Sharing of biases over multiple simulations is requested, but no bias is marked "
728 "as shared (share-group > 0)");
729 }
730
731 /* mdrun does not support this (yet), but will check again */
732 if (haveBiasSharingWithinSimulation(awhParams))
733 {
734 warning(wi,
735 "You have shared biases within a single simulation, but mdrun does not support "
736 "this (yet)");
737 }
738 }
739 } // namespace
740
readAwhParams(std::vector<t_inpfile> * inp,warninp_t wi)741 AwhParams* readAwhParams(std::vector<t_inpfile>* inp, warninp_t wi)
742 {
743 AwhParams* awhParams;
744 snew(awhParams, 1);
745 std::string opt;
746
747 /* Parameters common for all biases */
748
749 printStringNoNewline(inp, "The way to apply the biasing potential: convolved or umbrella");
750 opt = "awh-potential";
751 awhParams->ePotential = get_eeenum(inp, opt, eawhpotential_names, wi);
752
753 printStringNoNewline(inp,
754 "The random seed used for sampling the umbrella center in the case of "
755 "umbrella type potential");
756 opt = "awh-seed";
757 awhParams->seed = get_eint(inp, opt, -1, wi);
758 if (awhParams->seed == -1)
759 {
760 awhParams->seed = static_cast<int>(gmx::makeRandomSeed());
761 fprintf(stderr, "Setting the AWH bias MC random seed to %" PRId64 "\n", awhParams->seed);
762 }
763
764 printStringNoNewline(inp, "Data output interval in number of steps");
765 opt = "awh-nstout";
766 awhParams->nstOut = get_eint(inp, opt, 100000, wi);
767
768 printStringNoNewline(inp, "Coordinate sampling interval in number of steps");
769 opt = "awh-nstsample";
770 awhParams->nstSampleCoord = get_eint(inp, opt, 10, wi);
771
772 printStringNoNewline(inp, "Free energy and bias update interval in number of samples");
773 opt = "awh-nsamples-update";
774 awhParams->numSamplesUpdateFreeEnergy = get_eint(inp, opt, 10, wi);
775
776 printStringNoNewline(
777 inp, "When true, biases with share-group>0 are shared between multiple simulations");
778 opt = "awh-share-multisim";
779 awhParams->shareBiasMultisim = (get_eeenum(inp, opt, yesno_names, wi) != 0);
780
781 printStringNoNewline(inp, "The number of independent AWH biases");
782 opt = "awh-nbias";
783 awhParams->numBias = get_eint(inp, opt, 1, wi);
784 /* Check this before starting to read the AWH biases. */
785 if (awhParams->numBias <= 0)
786 {
787 gmx_fatal(FARGS, "%s needs to be an integer > 0", opt.c_str());
788 }
789
790 /* Read the parameters specific to each AWH bias */
791 snew(awhParams->awhBiasParams, awhParams->numBias);
792
793 for (int k = 0; k < awhParams->numBias; k++)
794 {
795 bool bComment = (k == 0);
796 std::string prefixawh = formatString("awh%d", k + 1);
797 readBiasParams(inp, &awhParams->awhBiasParams[k], prefixawh, wi, bComment);
798 }
799
800 return awhParams;
801 }
802
checkAwhParams(const AwhParams * awhParams,const t_inputrec * ir,warninp_t wi)803 void checkAwhParams(const AwhParams* awhParams, const t_inputrec* ir, warninp_t wi)
804 {
805 std::string opt;
806
807 checkMtsConsistency(*ir, wi);
808
809 opt = "awh-nstout";
810 if (awhParams->nstOut <= 0)
811 {
812 auto message = formatString("Not writing AWH output with AWH (%s = %d) does not make sense",
813 opt.c_str(), awhParams->nstOut);
814 warning_error(wi, message);
815 }
816 /* This restriction can be removed by changing a flag of print_ebin() */
817 if (ir->nstenergy == 0 || awhParams->nstOut % ir->nstenergy != 0)
818 {
819 auto message = formatString("%s (%d) should be a multiple of nstenergy (%d)", opt.c_str(),
820 awhParams->nstOut, ir->nstenergy);
821 warning_error(wi, message);
822 }
823
824 opt = "awh-nsamples-update";
825 if (awhParams->numSamplesUpdateFreeEnergy <= 0)
826 {
827 warning_error(wi, opt + " needs to be an integer > 0");
828 }
829
830 bool haveFepLambdaDim = false;
831 for (int k = 0; k < awhParams->numBias; k++)
832 {
833 std::string prefixawh = formatString("awh%d", k + 1);
834 checkBiasParams(&awhParams->awhBiasParams[k], prefixawh, ir, wi);
835 /* Check if there is a FEP lambda dimension. */
836 for (int l = 0; l < awhParams->awhBiasParams[k].ndim; l++)
837 {
838 if (awhParams->awhBiasParams[k].dimParams[l].eCoordProvider == eawhcoordproviderFREE_ENERGY_LAMBDA)
839 {
840 haveFepLambdaDim = true;
841 break;
842 }
843 }
844 }
845
846 if (haveFepLambdaDim)
847 {
848 if (awhParams->nstSampleCoord % ir->nstcalcenergy != 0)
849 {
850 opt = "awh-nstsample";
851 auto message = formatString(
852 "%s (%d) should be a multiple of nstcalcenergy (%d) when using AWH for "
853 "sampling an FEP lambda dimension",
854 opt.c_str(), awhParams->nstSampleCoord, ir->nstcalcenergy);
855 warning_error(wi, message);
856 }
857 if (awhParams->ePotential != eawhpotentialUMBRELLA)
858 {
859 opt = "awh-potential";
860 auto message = formatString(
861 "%s (%s) must be set to %s when using AWH for sampling an FEP lambda dimension",
862 opt.c_str(), eawhpotential_names[awhParams->ePotential],
863 eawhpotential_names[eawhpotentialUMBRELLA]);
864 warning_error(wi, message);
865 }
866 }
867
868 /* Do a final consistency check before returning */
869 checkInputConsistencyAwh(*awhParams, wi);
870
871 if (ir->init_step != 0)
872 {
873 warning_error(wi, "With AWH init-step should be 0");
874 }
875 }
876
877 /*! \brief
878 * Gets the period of a pull coordinate.
879 *
880 * \param[in] pullCoordParams The parameters for the pull coordinate.
881 * \param[in] pbc The PBC setup
882 * \param[in] intervalLength The length of the AWH interval for this pull coordinate
883 * \returns the period (or 0 if not periodic).
884 */
get_pull_coord_period(const t_pull_coord & pullCoordParams,const t_pbc & pbc,const real intervalLength)885 static double get_pull_coord_period(const t_pull_coord& pullCoordParams, const t_pbc& pbc, const real intervalLength)
886 {
887 double period = 0;
888
889 if (pullCoordParams.eGeom == epullgDIR)
890 {
891 const real margin = 0.001;
892 // Make dims periodic when the interval covers > 95%
893 const real periodicFraction = 0.95;
894
895 // Check if the pull direction is along a box vector
896 for (int dim = 0; dim < pbc.ndim_ePBC; dim++)
897 {
898 const real boxLength = norm(pbc.box[dim]);
899 const real innerProduct = iprod(pullCoordParams.vec, pbc.box[dim]);
900 if (innerProduct >= (1 - margin) * boxLength && innerProduct <= (1 + margin) * boxLength)
901 {
902 if (intervalLength > (1 + margin) * boxLength)
903 {
904 gmx_fatal(FARGS,
905 "The AWH interval (%f nm) for a pull coordinate is larger than the "
906 "box size (%f nm)",
907 intervalLength, boxLength);
908 }
909
910 if (intervalLength > periodicFraction * boxLength)
911 {
912 period = boxLength;
913 }
914 }
915 }
916 }
917 else if (pullCoordParams.eGeom == epullgDIHEDRAL)
918 {
919 /* The dihedral angle is periodic in -180 to 180 deg */
920 period = 360;
921 }
922
923 return period;
924 }
925
926 /*! \brief
927 * Checks if the given interval is defined in the correct periodic interval.
928 *
929 * \param[in] origin Start value of interval.
930 * \param[in] end End value of interval.
931 * \param[in] period Period (or 0 if not periodic).
932 * \returns true if the end point values are in the correct periodic interval.
933 */
intervalIsInPeriodicInterval(double origin,double end,double period)934 static bool intervalIsInPeriodicInterval(double origin, double end, double period)
935 {
936 return (period == 0) || (std::fabs(origin) <= 0.5 * period && std::fabs(end) <= 0.5 * period);
937 }
938
939 /*! \brief
940 * Checks if a value is within an interval.
941 *
942 * \param[in] origin Start value of interval.
943 * \param[in] end End value of interval.
944 * \param[in] period Period (or 0 if not periodic).
945 * \param[in] value Value to check.
946 * \returns true if the value is within the interval.
947 */
valueIsInInterval(double origin,double end,double period,double value)948 static bool valueIsInInterval(double origin, double end, double period, double value)
949 {
950 bool bIn_interval;
951
952 if (period > 0)
953 {
954 if (origin < end)
955 {
956 /* The interval closes within the periodic interval */
957 bIn_interval = (value >= origin) && (value <= end);
958 }
959 else
960 {
961 /* The interval wraps around the periodic boundary */
962 bIn_interval = ((value >= origin) && (value <= 0.5 * period))
963 || ((value >= -0.5 * period) && (value <= end));
964 }
965 }
966 else
967 {
968 bIn_interval = (value >= origin) && (value <= end);
969 }
970
971 return bIn_interval;
972 }
973
974 /*! \brief
975 * Check if the starting configuration is consistent with the given interval.
976 *
977 * \param[in] awhParams AWH parameters.
978 * \param[in,out] wi Struct for bookeeping warnings.
979 */
checkInputConsistencyInterval(const AwhParams * awhParams,warninp_t wi)980 static void checkInputConsistencyInterval(const AwhParams* awhParams, warninp_t wi)
981 {
982 for (int k = 0; k < awhParams->numBias; k++)
983 {
984 AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
985 for (int d = 0; d < awhBiasParams->ndim; d++)
986 {
987 AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
988 int coordIndex = dimParams->coordIndex;
989 double origin = dimParams->origin, end = dimParams->end, period = dimParams->period;
990 double coordValueInit = dimParams->coordValueInit;
991
992 if ((period == 0) && (origin > end))
993 {
994 gmx_fatal(FARGS,
995 "For the non-periodic pull coordinates awh%d-dim%d-start (%f) cannot be "
996 "larger than awh%d-dim%d-end (%f)",
997 k + 1, d + 1, origin, k + 1, d + 1, end);
998 }
999
1000 /* Currently we assume symmetric periodic intervals, meaning we use [-period/2, period/2] as the reference interval.
1001 Make sure the AWH interval is within this reference interval.
1002
1003 Note: we could fairly simply allow using a more general interval (e.g. [x, x + period]) but it complicates
1004 things slightly and I don't see that there is a great need for it. It would also mean that the interval would
1005 depend on AWH input. Also, for dihedral angles you would always want the reference interval to be -180, +180,
1006 independent of AWH parameters.
1007 */
1008 if (!intervalIsInPeriodicInterval(origin, end, period))
1009 {
1010 gmx_fatal(FARGS,
1011 "When using AWH with periodic pull coordinate geometries "
1012 "awh%d-dim%d-start (%.8g) and "
1013 "awh%d-dim%d-end (%.8g) should cover at most one period (%.8g) and take "
1014 "values in between "
1015 "minus half a period and plus half a period, i.e. in the interval [%.8g, "
1016 "%.8g].",
1017 k + 1, d + 1, origin, k + 1, d + 1, end, period, -0.5 * period, 0.5 * period);
1018 }
1019
1020 /* Warn if the pull initial coordinate value is not in the grid */
1021 if (!valueIsInInterval(origin, end, period, coordValueInit))
1022 {
1023 auto message = formatString(
1024 "The initial coordinate value (%.8g) for pull coordinate index %d falls "
1025 "outside "
1026 "of the sampling nterval awh%d-dim%d-start (%.8g) to awh%d-dim%d-end "
1027 "(%.8g). "
1028 "This can lead to large initial forces pulling the coordinate towards the "
1029 "sampling interval.",
1030 coordValueInit, coordIndex + 1, k + 1, d + 1, origin, k + 1, d + 1, end);
1031 warning(wi, message);
1032 }
1033 }
1034 }
1035 }
1036
1037 /*! \brief
1038 * Sets AWH parameters, for one AWH pull dimension.
1039 *
1040 * \param[in,out] dimParams AWH dimension parameters.
1041 * \param[in] biasIndex The index of the bias containing this AWH pull dimension.
1042 * \param[in] dimIndex The index of this AWH pull dimension.
1043 * \param[in] pull_params Pull parameters.
1044 * \param[in,out] pull_work Pull working struct to register AWH bias in.
1045 * \param[in] pbc A pbc information structure.
1046 * \param[in] compressibility Compressibility matrix for pressure coupling,
1047 * pass all 0 without pressure coupling.
1048 * \param[in,out] wi Struct for bookeeping warnings.
1049 *
1050 */
setStateDependentAwhPullDimParams(AwhDimParams * dimParams,const int biasIndex,const int dimIndex,const pull_params_t * pull_params,pull_t * pull_work,const t_pbc & pbc,const tensor & compressibility,warninp_t wi)1051 static void setStateDependentAwhPullDimParams(AwhDimParams* dimParams,
1052 const int biasIndex,
1053 const int dimIndex,
1054 const pull_params_t* pull_params,
1055 pull_t* pull_work,
1056 const t_pbc& pbc,
1057 const tensor& compressibility,
1058 warninp_t wi)
1059 {
1060 const t_pull_coord& pullCoordParams = pull_params->coord[dimParams->coordIndex];
1061
1062 if (pullCoordParams.eGeom == epullgDIRPBC)
1063 {
1064 gmx_fatal(FARGS,
1065 "AWH does not support pull geometry '%s'. "
1066 "If the maximum distance between the groups is always "
1067 "less than half the box size, "
1068 "you can use geometry '%s' instead.",
1069 EPULLGEOM(epullgDIRPBC), EPULLGEOM(epullgDIR));
1070 }
1071
1072 dimParams->period = get_pull_coord_period(pullCoordParams, pbc, dimParams->end - dimParams->origin);
1073 // We would like to check for scaling, but we don't have the full inputrec available here
1074 if (dimParams->period > 0
1075 && !(pullCoordParams.eGeom == epullgANGLE || pullCoordParams.eGeom == epullgDIHEDRAL))
1076 {
1077 bool coordIsScaled = false;
1078 for (int d2 = 0; d2 < DIM; d2++)
1079 {
1080 if (pullCoordParams.vec[d2] != 0 && norm2(compressibility[d2]) != 0)
1081 {
1082 coordIsScaled = true;
1083 }
1084 }
1085 if (coordIsScaled)
1086 {
1087 std::string mesg = gmx::formatString(
1088 "AWH dimension %d of bias %d is periodic with pull geometry '%s', "
1089 "while you should are applying pressure scaling to the "
1090 "corresponding box vector, this is not supported.",
1091 dimIndex + 1, biasIndex + 1, EPULLGEOM(pullCoordParams.eGeom));
1092 warning(wi, mesg.c_str());
1093 }
1094 }
1095
1096 /* The initial coordinate value, converted to external user units. */
1097 dimParams->coordValueInit = get_pull_coord_value(pull_work, dimParams->coordIndex, &pbc);
1098
1099 dimParams->coordValueInit *= pull_conversion_factor_internal2userinput(&pullCoordParams);
1100 }
1101
setStateDependentAwhParams(AwhParams * awhParams,const pull_params_t & pull_params,pull_t * pull_work,const matrix box,PbcType pbcType,const tensor & compressibility,const t_grpopts * inputrecGroupOptions,const real initLambda,const gmx_mtop_t & mtop,warninp_t wi)1102 void setStateDependentAwhParams(AwhParams* awhParams,
1103 const pull_params_t& pull_params,
1104 pull_t* pull_work,
1105 const matrix box,
1106 PbcType pbcType,
1107 const tensor& compressibility,
1108 const t_grpopts* inputrecGroupOptions,
1109 const real initLambda,
1110 const gmx_mtop_t& mtop,
1111 warninp_t wi)
1112 {
1113 /* The temperature is not really state depenendent but is not known
1114 * when read_awhParams is called (in get ir).
1115 * It is known first after do_index has been called in grompp.cpp.
1116 */
1117 if (inputrecGroupOptions->ref_t == nullptr || inputrecGroupOptions->ref_t[0] <= 0)
1118 {
1119 gmx_fatal(FARGS, "AWH biasing is only supported for temperatures > 0");
1120 }
1121 for (int i = 1; i < inputrecGroupOptions->ngtc; i++)
1122 {
1123 if (inputrecGroupOptions->ref_t[i] != inputrecGroupOptions->ref_t[0])
1124 {
1125 gmx_fatal(FARGS,
1126 "AWH biasing is currently only supported for identical temperatures for all "
1127 "temperature coupling groups");
1128 }
1129 }
1130
1131 t_pbc pbc;
1132 set_pbc(&pbc, pbcType, box);
1133
1134 for (int k = 0; k < awhParams->numBias; k++)
1135 {
1136 AwhBiasParams* awhBiasParams = &awhParams->awhBiasParams[k];
1137 for (int d = 0; d < awhBiasParams->ndim; d++)
1138 {
1139 AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
1140 if (dimParams->eCoordProvider == eawhcoordproviderPULL)
1141 {
1142 setStateDependentAwhPullDimParams(dimParams, k, d, &pull_params, pull_work, pbc,
1143 compressibility, wi);
1144 }
1145 else
1146 {
1147 dimParams->coordValueInit = initLambda;
1148 checkFepLambdaDimDecouplingConsistency(mtop, wi);
1149 }
1150 }
1151 }
1152 checkInputConsistencyInterval(awhParams, wi);
1153
1154 /* Register AWH as external potential with pull (for AWH dimensions that use pull as
1155 * reaction coordinate provider) to check consistency. */
1156 Awh::registerAwhWithPull(*awhParams, pull_work);
1157 }
1158
1159 } // namespace gmx
1160