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, 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.
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.
52 #include "gromacs/applied_forces/awh/read_params.h"
53 #include "gromacs/fileio/readinp.h"
54 #include "gromacs/fileio/warninp.h"
55 #include "gromacs/gmxlib/network.h"
56 #include "gromacs/gmxpreprocess/toputil.h"
57 #include "gromacs/math/functions.h"
58 #include "gromacs/math/units.h"
59 #include "gromacs/math/utilities.h"
60 #include "gromacs/math/vec.h"
61 #include "gromacs/mdlib/calc_verletbuf.h"
62 #include "gromacs/mdlib/vcm.h"
63 #include "gromacs/mdrun/mdmodules.h"
64 #include "gromacs/mdtypes/awh_params.h"
65 #include "gromacs/mdtypes/inputrec.h"
66 #include "gromacs/mdtypes/md_enums.h"
67 #include "gromacs/mdtypes/multipletimestepping.h"
68 #include "gromacs/mdtypes/pull_params.h"
69 #include "gromacs/options/options.h"
70 #include "gromacs/options/treesupport.h"
71 #include "gromacs/pbcutil/pbc.h"
72 #include "gromacs/selection/indexutil.h"
73 #include "gromacs/topology/block.h"
74 #include "gromacs/topology/ifunc.h"
75 #include "gromacs/topology/index.h"
76 #include "gromacs/topology/mtop_util.h"
77 #include "gromacs/topology/symtab.h"
78 #include "gromacs/topology/topology.h"
79 #include "gromacs/utility/arrayref.h"
80 #include "gromacs/utility/cstringutil.h"
81 #include "gromacs/utility/exceptions.h"
82 #include "gromacs/utility/fatalerror.h"
83 #include "gromacs/utility/filestream.h"
84 #include "gromacs/utility/gmxassert.h"
85 #include "gromacs/utility/ikeyvaluetreeerror.h"
86 #include "gromacs/utility/keyvaluetree.h"
87 #include "gromacs/utility/keyvaluetreebuilder.h"
88 #include "gromacs/utility/keyvaluetreemdpwriter.h"
89 #include "gromacs/utility/keyvaluetreetransform.h"
90 #include "gromacs/utility/mdmodulesnotifiers.h"
91 #include "gromacs/utility/smalloc.h"
92 #include "gromacs/utility/strconvert.h"
93 #include "gromacs/utility/stringcompare.h"
94 #include "gromacs/utility/stringutil.h"
95 #include "gromacs/utility/textwriter.h"
100 using gmx::BasicVector;
102 /* Resource parameters
103 * Do not change any of these until you read the instruction
104 * in readinp.h. Some cpp's do not take spaces after the backslash
105 * (like the c-shell), which will give you a very weird compiler
109 struct gmx_inputrec_strings
111 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], freeze[STRLEN], frdim[STRLEN],
112 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
113 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
114 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
115 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::string> fep_lambda;
116 char lambda_weights[STRLEN];
117 std::vector<std::string> pullGroupNames;
118 std::vector<std::string> rotateGroupNames;
119 char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
122 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
123 static gmx_inputrec_strings* inputrecStrings = nullptr;
125 void init_inputrec_strings()
130 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
131 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
133 inputrecStrings = new gmx_inputrec_strings();
136 void done_inputrec_strings()
138 delete inputrecStrings;
139 inputrecStrings = nullptr;
145 egrptpALL, /* All particles have to be a member of a group. */
146 egrptpALL_GENREST, /* A rest group with name is generated for particles *
147 * that are not part of any group. */
148 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
149 * for the rest group. */
150 egrptpONE /* Merge all selected groups into one group, *
151 * make a rest group for the remaining particles. */
154 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
155 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
156 "h-angles", "all-angles", nullptr };
158 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
159 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
161 static void getSimTemps(int ntemps, t_simtemp* simtemp, gmx::ArrayRef<double> temperature_lambdas)
166 for (i = 0; i < ntemps; i++)
168 /* simple linear scaling -- allows more control */
169 if (simtemp->eSimTempScale == SimulatedTempering::Linear)
171 simtemp->temperatures[i] =
173 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
175 else if (simtemp->eSimTempScale
176 == SimulatedTempering::Geometric) /* should give roughly equal acceptance for constant heat capacity . . . */
178 simtemp->temperatures[i] = simtemp->simtemp_low
179 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
180 static_cast<real>((1.0 * i) / (ntemps - 1)));
182 else if (simtemp->eSimTempScale == SimulatedTempering::Exponential)
184 simtemp->temperatures[i] = simtemp->simtemp_low
185 + (simtemp->simtemp_high - simtemp->simtemp_low)
186 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
191 sprintf(errorstr, "eSimTempScale=%s not defined", enumValueToString(simtemp->eSimTempScale));
192 gmx_fatal(FARGS, "%s", errorstr);
198 static void _low_check(bool b, const char* s, warninp_t wi)
202 warning_error(wi, s);
206 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
210 if (*p > 0 && *p % nst != 0)
212 /* Round up to the next multiple of nst */
213 *p = ((*p) / nst + 1) * nst;
214 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
219 //! Convert legacy mdp entries to modern ones.
220 static void process_interaction_modifier(InteractionModifiers* eintmod)
222 if (*eintmod == InteractionModifiers::PotShiftVerletUnsupported)
224 *eintmod = InteractionModifiers::PotShift;
228 void check_ir(const char* mdparin,
229 const gmx::MDModulesNotifiers& mdModulesNotifiers,
233 /* Check internal consistency.
234 * NOTE: index groups are not set here yet, don't check things
235 * like temperature coupling group options here, but in triple_check
238 /* Strange macro: first one fills the err_buf, and then one can check
239 * the condition, which will print the message and increase the error
242 #define CHECK(b) _low_check(b, err_buf, wi)
243 char err_buf[256], warn_buf[STRLEN];
246 t_lambda* fep = ir->fepvals.get();
247 t_expanded* expand = ir->expandedvals.get();
249 set_warning_line(wi, mdparin, -1);
251 /* We cannot check MTS requirements with an invalid MTS setup
252 * and we will already have generated errors with an invalid MTS setup.
254 if (gmx::haveValidMtsSetup(*ir))
256 std::vector<std::string> errorMessages = gmx::checkMtsRequirements(*ir);
258 for (const auto& errorMessage : errorMessages)
260 warning_error(wi, errorMessage.c_str());
264 if (ir->coulombtype == CoulombInteractionType::RFNecUnsupported)
266 std::string message =
267 gmx::formatString("%s electrostatics is no longer supported",
268 enumValueToString(CoulombInteractionType::RFNecUnsupported));
269 warning_error(wi, message);
272 /* BASIC CUT-OFF STUFF */
273 if (ir->rcoulomb < 0)
275 warning_error(wi, "rcoulomb should be >= 0");
279 warning_error(wi, "rvdw should be >= 0");
281 if (ir->rlist < 0 && !(ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0))
283 warning_error(wi, "rlist should be >= 0");
286 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
287 "neighbour-list update scheme for efficient buffering for improved energy "
288 "conservation, please use the Verlet cut-off scheme instead.)");
289 CHECK(ir->nstlist < 0);
291 process_interaction_modifier(&ir->coulomb_modifier);
292 process_interaction_modifier(&ir->vdw_modifier);
294 if (ir->cutoff_scheme == CutoffScheme::Group)
297 "The group cutoff scheme has been removed since GROMACS 2020. "
298 "Please use the Verlet cutoff scheme.");
300 if (ir->cutoff_scheme == CutoffScheme::Verlet)
304 /* Normal Verlet type neighbor-list, currently only limited feature support */
305 if (inputrec2nboundeddim(ir) < 3)
307 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
310 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
311 if (ir->rcoulomb != ir->rvdw)
313 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
314 // for PME load balancing, we can support this exception.
315 bool bUsesPmeTwinRangeKernel =
316 (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == VanDerWaalsType::Cut
317 && ir->rcoulomb > ir->rvdw);
318 if (!bUsesPmeTwinRangeKernel)
321 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
322 "rcoulomb>rvdw with PME electrostatics)");
326 if (ir->vdwtype == VanDerWaalsType::Shift || ir->vdwtype == VanDerWaalsType::Switch)
328 if (ir->vdw_modifier == InteractionModifiers::None
329 || ir->vdw_modifier == InteractionModifiers::PotShift)
332 (ir->vdwtype == VanDerWaalsType::Shift ? InteractionModifiers::ForceSwitch
333 : InteractionModifiers::PotSwitch);
336 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
338 enumValueToString(ir->vdwtype),
339 enumValueToString(VanDerWaalsType::Cut),
340 enumValueToString(ir->vdw_modifier));
341 warning_note(wi, warn_buf);
343 ir->vdwtype = VanDerWaalsType::Cut;
348 "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
349 enumValueToString(ir->vdwtype),
350 enumValueToString(ir->vdw_modifier));
351 warning_error(wi, warn_buf);
355 if (!(ir->vdwtype == VanDerWaalsType::Cut || ir->vdwtype == VanDerWaalsType::Pme))
358 "With Verlet lists only cut-off and PME LJ interactions are supported");
360 if (!(ir->coulombtype == CoulombInteractionType::Cut || EEL_RF(ir->coulombtype)
361 || EEL_PME(ir->coulombtype) || ir->coulombtype == CoulombInteractionType::Ewald))
364 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
365 "electrostatics are supported");
367 if (!(ir->coulomb_modifier == InteractionModifiers::None
368 || ir->coulomb_modifier == InteractionModifiers::PotShift))
370 sprintf(warn_buf, "coulomb_modifier=%s is not supported", enumValueToString(ir->coulomb_modifier));
371 warning_error(wi, warn_buf);
374 if (EEL_USER(ir->coulombtype))
377 "Coulomb type %s is not supported with the verlet scheme",
378 enumValueToString(ir->coulombtype));
379 warning_error(wi, warn_buf);
382 if (ir->nstlist <= 0)
384 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
387 if (ir->nstlist < 10)
390 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
391 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
395 rc_max = std::max(ir->rvdw, ir->rcoulomb);
399 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
400 ir->verletbuf_tol = 0;
403 else if (ir->verletbuf_tol <= 0)
405 if (ir->verletbuf_tol == 0)
407 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
410 if (ir->rlist < rc_max)
413 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
416 if (ir->rlist == rc_max && ir->nstlist > 1)
420 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
421 "buffer. The cluster pair list does have a buffering effect, but choosing "
422 "a larger rlist might be necessary for good energy conservation.");
427 if (ir->rlist > rc_max)
430 "You have set rlist larger than the interaction cut-off, but you also "
431 "have verlet-buffer-tolerance > 0. Will set rlist using "
432 "verlet-buffer-tolerance.");
435 if (ir->nstlist == 1)
437 /* No buffer required */
442 if (EI_DYNAMICS(ir->eI))
444 if (inputrec2nboundeddim(ir) < 3)
447 "The box volume is required for calculating rlist from the "
448 "energy drift with verlet-buffer-tolerance > 0. You are "
449 "using at least one unbounded dimension, so no volume can be "
450 "computed. Either use a finite box, or set rlist yourself "
451 "together with verlet-buffer-tolerance = -1.");
453 /* Set rlist temporarily so we can continue processing */
458 /* Set the buffer to 5% of the cut-off */
459 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
465 /* GENERAL INTEGRATOR STUFF */
468 if (ir->etc != TemperatureCoupling::No)
470 if (EI_RANDOM(ir->eI))
473 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
474 "implicitly. See the documentation for more information on which "
475 "parameters affect temperature for %s.",
476 enumValueToString(ir->etc),
477 enumValueToString(ir->eI),
478 enumValueToString(ir->eI));
483 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
485 enumValueToString(ir->etc),
486 enumValueToString(ir->eI));
488 warning_note(wi, warn_buf);
490 ir->etc = TemperatureCoupling::No;
492 if (ir->eI == IntegrationAlgorithm::VVAK)
495 "Integrator method %s is implemented primarily for validation purposes; for "
496 "molecular dynamics, you should probably be using %s or %s",
497 enumValueToString(IntegrationAlgorithm::VVAK),
498 enumValueToString(IntegrationAlgorithm::MD),
499 enumValueToString(IntegrationAlgorithm::VV));
500 warning_note(wi, warn_buf);
502 if (!EI_DYNAMICS(ir->eI))
504 if (ir->epc != PressureCoupling::No)
507 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
508 enumValueToString(ir->epc),
509 enumValueToString(ir->eI));
510 warning_note(wi, warn_buf);
512 ir->epc = PressureCoupling::No;
514 if (EI_DYNAMICS(ir->eI))
516 if (ir->nstcalcenergy < 0)
518 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
519 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
521 /* nstcalcenergy larger than nstener does not make sense.
522 * We ideally want nstcalcenergy=nstener.
526 ir->nstcalcenergy = std::gcd(ir->nstenergy, ir->nstlist);
530 ir->nstcalcenergy = ir->nstenergy;
534 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
535 || (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl > 0
536 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
539 const char* nsten = "nstenergy";
540 const char* nstdh = "nstdhdl";
541 const char* min_name = nsten;
542 int min_nst = ir->nstenergy;
544 /* find the smallest of ( nstenergy, nstdhdl ) */
545 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl > 0
546 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
548 min_nst = ir->fepvals->nstdhdl;
551 /* If the user sets nstenergy small, we should respect that */
552 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy, min_name, min_nst);
553 warning_note(wi, warn_buf);
554 ir->nstcalcenergy = min_nst;
557 if (ir->epc != PressureCoupling::No)
559 if (ir->nstpcouple < 0)
561 ir->nstpcouple = ir_optimal_nstpcouple(ir);
563 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
566 "With multiple time stepping, nstpcouple should be a mutiple of "
571 if (ir->nstcalcenergy > 0)
573 if (ir->efep != FreeEnergyPerturbationType::No)
575 /* nstdhdl should be a multiple of nstcalcenergy */
576 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
580 /* nstexpanded should be a multiple of nstcalcenergy */
581 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded", &ir->expandedvals->nstexpanded, wi);
583 /* for storing exact averages nstenergy should be
584 * a multiple of nstcalcenergy
586 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
589 // Inquire all MDModules, if their parameters match with the energy
590 // calculation frequency
591 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
592 mdModulesNotifiers.preProcessingNotifier_.notify(&energyCalculationFrequencyErrors);
594 // Emit all errors from the energy calculation frequency checks
595 for (const std::string& energyFrequencyErrorMessage :
596 energyCalculationFrequencyErrors.errorMessages())
598 warning_error(wi, energyFrequencyErrorMessage);
602 if (ir->nsteps == 0 && !ir->bContinuation)
605 "For a correct single-point energy evaluation with nsteps = 0, use "
606 "continuation = yes to avoid constraining the input coordinates.");
610 if ((EI_SD(ir->eI) || ir->eI == IntegrationAlgorithm::BD) && ir->bContinuation && ir->ld_seed != -1)
613 "You are doing a continuation with SD or BD, make sure that ld_seed is "
614 "different from the previous run (using ld_seed=-1 will ensure this)");
620 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
621 CHECK(ir->pbcType != PbcType::Xyz);
622 sprintf(err_buf, "with TPI nstlist should be larger than zero");
623 CHECK(ir->nstlist <= 0);
624 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
625 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
629 if ((opts->nshake > 0) && (opts->bMorse))
631 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
632 warning(wi, warn_buf);
635 if ((EI_SD(ir->eI) || ir->eI == IntegrationAlgorithm::BD) && ir->bContinuation && ir->ld_seed != -1)
638 "You are doing a continuation with SD or BD, make sure that ld_seed is "
639 "different from the previous run (using ld_seed=-1 will ensure this)");
641 /* verify simulated tempering options */
645 bool bAllTempZero = TRUE;
646 for (i = 0; i < fep->n_lambda; i++)
649 "Entry %d for %s must be between 0 and 1, instead is %g",
651 enumValueToString(FreeEnergyPerturbationCouplingType::Temperature),
652 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i]);
653 CHECK((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i] < 0)
654 || (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i]
656 if (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i] > 0)
658 bAllTempZero = FALSE;
661 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
662 CHECK(bAllTempZero == TRUE);
664 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
665 CHECK(ir->eI != IntegrationAlgorithm::VV);
667 /* check compatability of the temperature coupling with simulated tempering */
669 if (ir->etc == TemperatureCoupling::NoseHoover)
672 "Nose-Hoover based temperature control such as [%s] my not be "
673 "entirelyconsistent with simulated tempering",
674 enumValueToString(ir->etc));
675 warning_note(wi, warn_buf);
678 /* check that the temperatures make sense */
681 "Higher simulated tempering temperature (%g) must be >= than the simulated "
682 "tempering lower temperature (%g)",
683 ir->simtempvals->simtemp_high,
684 ir->simtempvals->simtemp_low);
685 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
688 "Higher simulated tempering temperature (%g) must be >= zero",
689 ir->simtempvals->simtemp_high);
690 CHECK(ir->simtempvals->simtemp_high <= 0);
693 "Lower simulated tempering temperature (%g) must be >= zero",
694 ir->simtempvals->simtemp_low);
695 CHECK(ir->simtempvals->simtemp_low <= 0);
698 /* verify free energy options */
700 if (ir->efep != FreeEnergyPerturbationType::No)
702 fep = ir->fepvals.get();
703 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
704 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
707 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
709 static_cast<int>(fep->sc_r_power));
710 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
713 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
716 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
719 "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
721 CHECK(fep->delta_lambda > 0 && (ir->efep == FreeEnergyPerturbationType::Expanded));
723 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
724 CHECK(!(EI_VV(ir->eI)) && (ir->efep == FreeEnergyPerturbationType::Expanded));
726 sprintf(err_buf, "Free-energy not implemented for Ewald");
727 CHECK(ir->coulombtype == CoulombInteractionType::Ewald);
729 /* check validty of lambda inputs */
730 if (fep->n_lambda == 0)
732 /* Clear output in case of no states:*/
733 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
734 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
739 "initial thermodynamic state %d does not exist, only goes to %d",
742 CHECK((fep->init_fep_state >= fep->n_lambda));
746 "Lambda state must be set, either with init-lambda-state or with init-lambda");
747 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
750 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
751 "init-lambda-state or with init-lambda, but not both",
753 fep->init_fep_state);
754 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
757 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
761 for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
763 if (fep->separate_dvdl[i])
768 if (n_lambda_terms > 1)
771 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
772 "use init-lambda to set lambda state (except for slow growth). Use "
773 "init-lambda-state instead.");
774 warning(wi, warn_buf);
777 if (n_lambda_terms < 2 && fep->n_lambda > 0)
780 "init-lambda is deprecated for setting lambda state (except for slow "
781 "growth). Use init-lambda-state instead.");
785 for (j = 0; j < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); j++)
787 for (i = 0; i < fep->n_lambda; i++)
789 auto enumValue = static_cast<FreeEnergyPerturbationCouplingType>(j);
791 "Entry %d for %s must be between 0 and 1, instead is %g",
793 enumValueToString(enumValue),
794 fep->all_lambda[j][i]);
795 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
799 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
801 for (i = 0; i < fep->n_lambda; i++)
804 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
805 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
806 "crashes, and is not supported.",
808 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i],
809 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i]);
810 CHECK((fep->sc_alpha > 0)
811 && (((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i] > 0.0)
812 && (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i] < 1.0))
813 && ((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i] > 0.0)
814 && (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i]
819 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
821 real sigma, lambda, r_sc;
824 /* Maximum estimate for A and B charges equal with lambda power 1 */
826 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
827 1.0 / fep->sc_r_power);
829 "With PME there is a minor soft core effect present at the cut-off, "
830 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
831 "energy conservation, but usually other effects dominate. With a common sigma "
832 "value of %g nm the fraction of the particle-particle potential at the cut-off "
833 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
839 warning_note(wi, warn_buf);
842 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
843 be treated differently, but that's the next step */
845 for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
847 auto enumValue = static_cast<FreeEnergyPerturbationCouplingType>(i);
848 for (j = 0; j < fep->n_lambda; j++)
850 sprintf(err_buf, "%s[%d] must be between 0 and 1", enumValueToString(enumValue), j);
851 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
856 if ((ir->bSimTemp) || (ir->efep == FreeEnergyPerturbationType::Expanded))
858 fep = ir->fepvals.get();
860 /* checking equilibration of weights inputs for validity */
863 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
865 expand->equil_n_at_lam,
866 enumValueToString(LambdaWeightWillReachEquilibrium::NumAtLambda));
867 CHECK((expand->equil_n_at_lam > 0)
868 && (expand->elmceq != LambdaWeightWillReachEquilibrium::NumAtLambda));
871 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
873 expand->equil_samples,
874 enumValueToString(LambdaWeightWillReachEquilibrium::Samples));
875 CHECK((expand->equil_samples > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Samples));
878 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
880 enumValueToString(LambdaWeightWillReachEquilibrium::Steps));
881 CHECK((expand->equil_steps > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Steps));
884 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
885 expand->equil_samples,
886 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta));
887 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::WLDelta));
890 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
892 enumValueToString(LambdaWeightWillReachEquilibrium::Ratio));
893 CHECK((expand->equil_ratio > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Ratio));
896 "weight-equil-number-all-lambda (%d) must be a positive integer if "
897 "lmc-weights-equil=%s",
898 expand->equil_n_at_lam,
899 enumValueToString(LambdaWeightWillReachEquilibrium::NumAtLambda));
900 CHECK((expand->equil_n_at_lam <= 0)
901 && (expand->elmceq == LambdaWeightWillReachEquilibrium::NumAtLambda));
904 "weight-equil-number-samples (%d) must be a positive integer if "
905 "lmc-weights-equil=%s",
906 expand->equil_samples,
907 enumValueToString(LambdaWeightWillReachEquilibrium::Samples));
908 CHECK((expand->equil_samples <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Samples));
911 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
913 enumValueToString(LambdaWeightWillReachEquilibrium::Steps));
914 CHECK((expand->equil_steps <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Steps));
917 "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
918 expand->equil_wl_delta,
919 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta));
920 CHECK((expand->equil_wl_delta <= 0)
921 && (expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta));
924 "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
926 enumValueToString(LambdaWeightWillReachEquilibrium::Ratio));
927 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Ratio));
930 "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
931 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta),
932 enumValueToString(LambdaWeightCalculation::WL),
933 enumValueToString(LambdaWeightCalculation::WWL));
934 CHECK((expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta) && (!EWL(expand->elamstats)));
936 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
937 CHECK((expand->lmc_repeats <= 0));
938 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
939 CHECK((expand->minvarmin <= 0));
940 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
941 CHECK((expand->c_range < 0));
943 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
946 expand->lmc_forced_nstart);
947 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
948 && (expand->elmcmove != LambdaMoveCalculation::No));
949 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
950 CHECK((expand->lmc_forced_nstart < 0));
952 "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
953 fep->init_fep_state);
954 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
956 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
957 CHECK((expand->init_wl_delta < 0));
958 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
959 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
960 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
961 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
963 /* if there is no temperature control, we need to specify an MC temperature */
964 if (!integratorHasReferenceTemperature(ir)
965 && (expand->elmcmove != LambdaMoveCalculation::No) && (expand->mc_temp <= 0.0))
968 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
969 "to a positive number");
970 warning_error(wi, err_buf);
972 if (expand->nstTij > 0)
974 sprintf(err_buf, "nstlog must be non-zero");
975 CHECK(ir->nstlog == 0);
976 // Avoid modulus by zero in the case that already triggered an error exit.
980 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
983 CHECK((expand->nstTij % ir->nstlog) != 0);
989 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
990 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
993 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
995 if (ir->pbcType == PbcType::No)
997 if (ir->epc != PressureCoupling::No)
999 warning(wi, "Turning off pressure coupling for vacuum system");
1000 ir->epc = PressureCoupling::No;
1006 "Can not have pressure coupling with pbc=%s",
1007 c_pbcTypeNames[ir->pbcType].c_str());
1008 CHECK(ir->epc != PressureCoupling::No);
1010 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
1011 CHECK(EEL_FULL(ir->coulombtype));
1014 "Can not have dispersion correction with pbc=%s",
1015 c_pbcTypeNames[ir->pbcType].c_str());
1016 CHECK(ir->eDispCorr != DispersionCorrectionType::No);
1019 if (ir->rlist == 0.0)
1022 "can only have neighborlist cut-off zero (=infinite)\n"
1023 "with coulombtype = %s or coulombtype = %s\n"
1024 "without periodic boundary conditions (pbc = %s) and\n"
1025 "rcoulomb and rvdw set to zero",
1026 enumValueToString(CoulombInteractionType::Cut),
1027 enumValueToString(CoulombInteractionType::User),
1028 c_pbcTypeNames[PbcType::No].c_str());
1029 CHECK(((ir->coulombtype != CoulombInteractionType::Cut)
1030 && (ir->coulombtype != CoulombInteractionType::User))
1031 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1033 if (ir->nstlist > 0)
1036 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1037 "nstype=simple and only one MPI rank");
1042 if (ir->nstcomm == 0)
1044 // TODO Change this behaviour. There should be exactly one way
1045 // to turn off an algorithm.
1046 ir->comm_mode = ComRemovalAlgorithm::No;
1048 if (ir->comm_mode != ComRemovalAlgorithm::No)
1050 if (ir->nstcomm < 0)
1052 // TODO Such input was once valid. Now that we've been
1053 // helpful for a few years, we should reject such input,
1054 // lest we have to support every historical decision
1057 "If you want to remove the rotation around the center of mass, you should set "
1058 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1059 "its absolute value");
1060 ir->nstcomm = abs(ir->nstcomm);
1063 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
1066 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
1067 "nstcomm to nstcalcenergy");
1068 ir->nstcomm = ir->nstcalcenergy;
1071 if (ir->comm_mode == ComRemovalAlgorithm::Angular)
1074 "Can not remove the rotation around the center of mass with periodic "
1076 CHECK(ir->bPeriodicMols);
1077 if (ir->pbcType != PbcType::No)
1080 "Removing the rotation around the center of mass in a periodic system, "
1081 "this can lead to artifacts. Only use this on a single (cluster of) "
1082 "molecules. This cluster should not cross periodic boundaries.");
1087 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No
1088 && ir->comm_mode != ComRemovalAlgorithm::Angular)
1091 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1092 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1094 enumValueToString(IntegrationAlgorithm::SD1));
1095 warning_note(wi, warn_buf);
1098 /* TEMPERATURE COUPLING */
1099 if (ir->etc == TemperatureCoupling::Yes)
1101 ir->etc = TemperatureCoupling::Berendsen;
1103 "Old option for temperature coupling given: "
1104 "changing \"yes\" to \"Berendsen\"\n");
1107 if ((ir->etc == TemperatureCoupling::NoseHoover) || (ir->epc == PressureCoupling::Mttk))
1109 if (ir->opts.nhchainlength < 1)
1112 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1114 ir->opts.nhchainlength);
1115 ir->opts.nhchainlength = 1;
1116 warning(wi, warn_buf);
1119 if (ir->etc == TemperatureCoupling::NoseHoover && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1123 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1124 ir->opts.nhchainlength = 1;
1129 ir->opts.nhchainlength = 0;
1132 if (ir->eI == IntegrationAlgorithm::VVAK)
1135 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1137 enumValueToString(IntegrationAlgorithm::VVAK));
1138 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1141 if (ETC_ANDERSEN(ir->etc))
1144 "%s temperature control not supported for integrator %s.",
1145 enumValueToString(ir->etc),
1146 enumValueToString(ir->eI));
1147 CHECK(!(EI_VV(ir->eI)));
1149 if (ir->nstcomm > 0 && (ir->etc == TemperatureCoupling::Andersen))
1152 "Center of mass removal not necessary for %s. All velocities of coupled "
1153 "groups are rerandomized periodically, so flying ice cube errors will not "
1155 enumValueToString(ir->etc));
1156 warning_note(wi, warn_buf);
1160 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1161 "randomized every time step",
1163 enumValueToString(ir->etc));
1164 CHECK(ir->nstcomm > 1 && (ir->etc == TemperatureCoupling::Andersen));
1167 if (ir->etc == TemperatureCoupling::Berendsen)
1170 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1171 "might want to consider using the %s thermostat.",
1172 enumValueToString(ir->etc),
1173 enumValueToString(TemperatureCoupling::VRescale));
1174 warning_note(wi, warn_buf);
1177 if ((ir->etc == TemperatureCoupling::NoseHoover || ETC_ANDERSEN(ir->etc))
1178 && ir->epc == PressureCoupling::Berendsen)
1181 "Using Berendsen pressure coupling invalidates the "
1182 "true ensemble for the thermostat");
1183 warning(wi, warn_buf);
1186 /* PRESSURE COUPLING */
1187 if (ir->epc == PressureCoupling::Isotropic)
1189 ir->epc = PressureCoupling::Berendsen;
1191 "Old option for pressure coupling given: "
1192 "changing \"Isotropic\" to \"Berendsen\"\n");
1195 if (ir->epc != PressureCoupling::No)
1197 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1199 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1200 CHECK(ir->tau_p <= 0);
1202 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1205 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1206 "times larger than nstpcouple*dt (%g)",
1207 enumValueToString(ir->epc),
1209 pcouple_min_integration_steps(ir->epc),
1211 warning(wi, warn_buf);
1215 "compressibility must be > 0 when using pressure"
1217 enumValueToString(ir->epc));
1218 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1219 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1220 && ir->compress[ZZ][YY] <= 0));
1222 if (PressureCoupling::ParrinelloRahman == ir->epc && opts->bGenVel)
1225 "You are generating velocities so I am assuming you "
1226 "are equilibrating a system. You are using "
1227 "%s pressure coupling, but this can be "
1228 "unstable for equilibration. If your system crashes, try "
1229 "equilibrating first with Berendsen pressure coupling. If "
1230 "you are not equilibrating the system, you can probably "
1231 "ignore this warning.",
1232 enumValueToString(ir->epc));
1233 warning(wi, warn_buf);
1239 if (ir->epc == PressureCoupling::Mttk)
1241 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1245 /* ELECTROSTATICS */
1246 /* More checks are in triple check (grompp.c) */
1248 if (ir->coulombtype == CoulombInteractionType::Switch)
1251 "coulombtype = %s is only for testing purposes and can lead to serious "
1252 "artifacts, advice: use coulombtype = %s",
1253 enumValueToString(ir->coulombtype),
1254 enumValueToString(CoulombInteractionType::RFZero));
1255 warning(wi, warn_buf);
1258 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1261 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1262 "format and exchanging epsilon-r and epsilon-rf",
1264 warning(wi, warn_buf);
1265 ir->epsilon_rf = ir->epsilon_r;
1266 ir->epsilon_r = 1.0;
1269 if (ir->epsilon_r == 0)
1272 "It is pointless to use long-range electrostatics with infinite relative "
1274 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1276 CHECK(EEL_FULL(ir->coulombtype));
1279 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1281 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1282 CHECK(ir->epsilon_r < 0);
1285 if (EEL_RF(ir->coulombtype))
1287 /* reaction field (at the cut-off) */
1289 if (ir->coulombtype == CoulombInteractionType::RFZero && ir->epsilon_rf != 0)
1292 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1293 enumValueToString(ir->coulombtype));
1294 warning(wi, warn_buf);
1295 ir->epsilon_rf = 0.0;
1298 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1299 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1300 if (ir->epsilon_rf == ir->epsilon_r)
1303 "Using epsilon-rf = epsilon-r with %s does not make sense",
1304 enumValueToString(ir->coulombtype));
1305 warning(wi, warn_buf);
1308 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1309 * means the interaction is zero outside rcoulomb, but it helps to
1310 * provide accurate energy conservation.
1312 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1314 if (ir_coulomb_switched(ir))
1317 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1318 "potential modifier options!",
1319 enumValueToString(ir->coulombtype));
1320 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1324 if (ir->coulombtype == CoulombInteractionType::Switch || ir->coulombtype == CoulombInteractionType::Shift)
1327 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1328 "secondary coulomb-modifier.");
1329 CHECK(ir->coulomb_modifier != InteractionModifiers::None);
1331 if (ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift)
1334 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1335 "secondary vdw-modifier.");
1336 CHECK(ir->vdw_modifier != InteractionModifiers::None);
1339 if (ir->coulombtype == CoulombInteractionType::Switch || ir->coulombtype == CoulombInteractionType::Shift
1340 || ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift)
1343 "The switch/shift interaction settings are just for compatibility; you will get "
1345 "performance from applying potential modifiers to your interactions!\n");
1346 warning_note(wi, warn_buf);
1349 if (ir->coulombtype == CoulombInteractionType::PmeSwitch
1350 || ir->coulomb_modifier == InteractionModifiers::PotSwitch)
1352 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1354 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1356 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1357 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1358 "will be good regardless, since ewald_rtol = %g.",
1360 ir->rcoulomb_switch,
1363 warning(wi, warn_buf);
1367 if (ir->vdwtype == VanDerWaalsType::Switch || ir->vdw_modifier == InteractionModifiers::PotSwitch)
1369 if (ir->rvdw_switch == 0)
1372 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1373 "potential. This suggests it was not set in the mdp, which can lead to large "
1374 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1375 "switching range.");
1376 warning(wi, warn_buf);
1380 if (EEL_FULL(ir->coulombtype))
1382 if (ir->coulombtype == CoulombInteractionType::PmeSwitch
1383 || ir->coulombtype == CoulombInteractionType::PmeUser
1384 || ir->coulombtype == CoulombInteractionType::PmeUserSwitch)
1387 "With coulombtype = %s, rcoulomb must be <= rlist",
1388 enumValueToString(ir->coulombtype));
1389 CHECK(ir->rcoulomb > ir->rlist);
1393 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1395 // TODO: Move these checks into the ewald module with the options class
1397 int orderMax = (ir->coulombtype == CoulombInteractionType::P3mAD ? 8 : 12);
1399 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1402 "With coulombtype = %s, you should have %d <= pme-order <= %d",
1403 enumValueToString(ir->coulombtype),
1406 warning_error(wi, warn_buf);
1410 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1412 if (ir->ewald_geometry == EwaldGeometry::ThreeD)
1415 "With pbc=%s you should use ewald-geometry=%s",
1416 c_pbcTypeNames[ir->pbcType].c_str(),
1417 enumValueToString(EwaldGeometry::ThreeDC));
1418 warning(wi, warn_buf);
1420 /* This check avoids extra pbc coding for exclusion corrections */
1421 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1422 CHECK(ir->wall_ewald_zfac < 2);
1424 if ((ir->ewald_geometry == EwaldGeometry::ThreeDC) && (ir->pbcType != PbcType::XY)
1425 && EEL_FULL(ir->coulombtype))
1428 "With %s and ewald_geometry = %s you should use pbc = %s",
1429 enumValueToString(ir->coulombtype),
1430 enumValueToString(EwaldGeometry::ThreeDC),
1431 c_pbcTypeNames[PbcType::XY].c_str());
1432 warning(wi, warn_buf);
1434 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1436 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1437 CHECK(ir->bPeriodicMols);
1438 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1439 warning_note(wi, warn_buf);
1441 "With epsilon_surface > 0 you can only use domain decomposition "
1442 "when there are only small molecules with all bonds constrained (mdrun will check "
1444 warning_note(wi, warn_buf);
1447 if (ir_vdw_switched(ir))
1449 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1450 CHECK(ir->rvdw_switch >= ir->rvdw);
1452 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1455 "You are applying a switch function to vdw forces or potentials from %g to %g "
1456 "nm, which is more than half the interaction range, whereas switch functions "
1457 "are intended to act only close to the cut-off.",
1460 warning_note(wi, warn_buf);
1464 if (ir->vdwtype == VanDerWaalsType::Pme)
1466 if (!(ir->vdw_modifier == InteractionModifiers::None
1467 || ir->vdw_modifier == InteractionModifiers::PotShift))
1470 "With vdwtype = %s, the only supported modifiers are %s and %s",
1471 enumValueToString(ir->vdwtype),
1472 enumValueToString(InteractionModifiers::PotShift),
1473 enumValueToString(InteractionModifiers::None));
1474 warning_error(wi, err_buf);
1478 if (ir->vdwtype == VanDerWaalsType::User && ir->eDispCorr != DispersionCorrectionType::No)
1481 "You have selected user tables with dispersion correction, the dispersion "
1482 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1483 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1484 "really want dispersion correction to -C6/r^6.");
1487 if (ir->eI == IntegrationAlgorithm::LBFGS
1488 && (ir->coulombtype == CoulombInteractionType::Cut || ir->vdwtype == VanDerWaalsType::Cut)
1491 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1494 if (ir->eI == IntegrationAlgorithm::LBFGS && ir->nbfgscorr <= 0)
1496 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1499 /* IMPLICIT SOLVENT */
1500 if (ir->coulombtype == CoulombInteractionType::GBNotused)
1502 sprintf(warn_buf, "Invalid option %s for coulombtype", enumValueToString(ir->coulombtype));
1503 warning_error(wi, warn_buf);
1508 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1513 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1516 // cosine acceleration is only supported in leap-frog
1517 if (ir->cos_accel != 0.0 && ir->eI != IntegrationAlgorithm::MD)
1519 warning_error(wi, "cos-acceleration is only supported by integrator = md");
1523 /* interpret a number of doubles from a string and put them in an array,
1524 after allocating space for them.
1525 str = the input string
1526 n = the (pre-allocated) number of doubles read
1527 r = the output array of doubles. */
1528 static std::vector<real> parse_n_real(const std::string& str, int* n, warninp_t wi)
1530 auto values = gmx::splitString(str);
1533 std::vector<real> r;
1534 for (int i = 0; i < *n; i++)
1538 r.emplace_back(gmx::fromString<real>(values[i]));
1540 catch (gmx::GromacsException&)
1543 "Invalid value " + values[i]
1544 + " in string in mdp file. Expected a real number.");
1551 static void do_fep_params(t_inputrec* ir, gmx::ArrayRef<std::string> fep_lambda, char weights[STRLEN], warninp_t wi)
1554 int i, j, max_n_lambda, nweights;
1555 t_lambda* fep = ir->fepvals.get();
1556 t_expanded* expand = ir->expandedvals.get();
1557 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::vector<real>> count_fep_lambdas;
1558 bool bOneLambda = TRUE;
1559 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, int> nfep;
1561 /* FEP input processing */
1562 /* first, identify the number of lambda values for each type.
1563 All that are nonzero must have the same number */
1565 for (auto i : keysOf(nfep))
1567 count_fep_lambdas[i] = parse_n_real(fep_lambda[static_cast<int>(i)], &(nfep[i]), wi);
1570 /* now, determine the number of components. All must be either zero, or equal. */
1573 for (auto i : keysOf(nfep))
1575 if (nfep[i] > max_n_lambda)
1577 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1578 must have the same number if its not zero.*/
1583 for (auto i : keysOf(nfep))
1587 ir->fepvals->separate_dvdl[i] = FALSE;
1589 else if (nfep[i] == max_n_lambda)
1591 if (i != FreeEnergyPerturbationCouplingType::Temperature) /* we treat this differently -- not really a reason to compute
1592 the derivative with respect to the temperature currently */
1594 ir->fepvals->separate_dvdl[i] = TRUE;
1600 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1603 enumValueToString(i),
1607 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1608 ir->fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Temperature] = FALSE;
1610 /* the number of lambdas is the number we've read in, which is either zero
1611 or the same for all */
1612 fep->n_lambda = max_n_lambda;
1614 /* if init_lambda is defined, we need to set lambda */
1615 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1617 ir->fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep] = TRUE;
1619 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1620 for (auto i : keysOf(nfep))
1622 fep->all_lambda[i].resize(fep->n_lambda);
1623 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1626 for (j = 0; j < fep->n_lambda; j++)
1628 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1633 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1634 internal bookkeeping -- for now, init_lambda */
1636 if ((nfep[FreeEnergyPerturbationCouplingType::Fep] == 0) && (fep->init_lambda >= 0))
1638 for (i = 0; i < fep->n_lambda; i++)
1640 fep->all_lambda[FreeEnergyPerturbationCouplingType::Fep][i] = fep->init_lambda;
1644 /* check to see if only a single component lambda is defined, and soft core is defined.
1645 In this case, turn on coulomb soft core */
1647 if (max_n_lambda == 0)
1653 for (auto i : keysOf(nfep))
1655 if ((nfep[i] != 0) && (i != FreeEnergyPerturbationCouplingType::Fep))
1661 if ((bOneLambda) && (fep->sc_alpha > 0))
1663 fep->bScCoul = TRUE;
1666 /* Fill in the others with the efptFEP if they are not explicitly
1667 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1668 they are all zero. */
1670 for (auto i : keysOf(nfep))
1672 if ((nfep[i] == 0) && (i != FreeEnergyPerturbationCouplingType::Fep))
1674 for (j = 0; j < fep->n_lambda; j++)
1676 fep->all_lambda[i][j] = fep->all_lambda[FreeEnergyPerturbationCouplingType::Fep][j];
1682 /* now read in the weights */
1683 expand->init_lambda_weights = parse_n_real(weights, &nweights, wi);
1686 expand->init_lambda_weights.resize(fep->n_lambda); /* initialize to zero */
1688 else if (nweights != fep->n_lambda)
1691 "Number of weights (%d) is not equal to number of lambda values (%d)",
1695 if ((expand->nstexpanded < 0) && (ir->efep != FreeEnergyPerturbationType::No))
1697 expand->nstexpanded = fep->nstdhdl;
1698 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1703 static void do_simtemp_params(t_inputrec* ir)
1705 ir->simtempvals->temperatures.resize(ir->fepvals->n_lambda);
1706 getSimTemps(ir->fepvals->n_lambda,
1707 ir->simtempvals.get(),
1708 ir->fepvals->all_lambda[FreeEnergyPerturbationCouplingType::Temperature]);
1711 template<typename T>
1712 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1715 for (const auto& input : inputs)
1719 outputs[i] = gmx::fromStdString<T>(input);
1721 catch (gmx::GromacsException&)
1723 auto message = gmx::formatString(
1724 "Invalid value for mdp option %s. %s should only consist of integers separated "
1728 warning_error(wi, message);
1734 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1737 for (const auto& input : inputs)
1741 outputs[i] = gmx::fromString<real>(input);
1743 catch (gmx::GromacsException&)
1745 auto message = gmx::formatString(
1746 "Invalid value for mdp option %s. %s should only consist of real numbers "
1747 "separated by spaces.",
1750 warning_error(wi, message);
1756 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1758 opts->wall_atomtype[0] = nullptr;
1759 opts->wall_atomtype[1] = nullptr;
1761 ir->wall_atomtype[0] = -1;
1762 ir->wall_atomtype[1] = -1;
1763 ir->wall_density[0] = 0;
1764 ir->wall_density[1] = 0;
1768 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1769 if (wallAtomTypes.size() != size_t(ir->nwall))
1772 "Expected %d elements for wall_atomtype, found %zu",
1774 wallAtomTypes.size());
1776 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1777 for (int i = 0; i < ir->nwall; i++)
1779 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1782 if (ir->wall_type == WallType::NineThree || ir->wall_type == WallType::TenFour)
1784 auto wallDensity = gmx::splitString(wall_density);
1785 if (wallDensity.size() != size_t(ir->nwall))
1788 "Expected %d elements for wall-density, found %zu",
1790 wallDensity.size());
1792 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1793 for (int i = 0; i < ir->nwall; i++)
1795 if (ir->wall_density[i] <= 0)
1797 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1804 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1808 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1809 for (int i = 0; i < nwall; i++)
1811 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1812 grps->emplace_back(groups->groupNames.size() - 1);
1817 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1819 /* read expanded ensemble parameters */
1820 printStringNewline(inp, "expanded ensemble variables");
1821 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1822 expand->elamstats = getEnum<LambdaWeightCalculation>(inp, "lmc-stats", wi);
1823 expand->elmcmove = getEnum<LambdaMoveCalculation>(inp, "lmc-move", wi);
1824 expand->elmceq = getEnum<LambdaWeightWillReachEquilibrium>(inp, "lmc-weights-equil", wi);
1825 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1826 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1827 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1828 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1829 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1830 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1831 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1832 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1833 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1834 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1835 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1836 expand->bSymmetrizedTMatrix =
1837 (getEnum<Boolean>(inp, "symmetrized-transition-matrix", wi) != Boolean::No);
1838 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1839 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1840 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1841 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1842 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1843 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1844 expand->bWLoneovert = (getEnum<Boolean>(inp, "wl-oneovert", wi) != Boolean::No);
1847 /*! \brief Return whether an end state with the given coupling-lambda
1848 * value describes fully-interacting VDW.
1850 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1851 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1853 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1855 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1861 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1864 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1866 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1868 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1871 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1872 std::string message = gmx::formatExceptionMessageToString(*ex);
1873 warning_error(wi_, message.c_str());
1878 std::string getOptionName(const gmx::KeyValueTreePath& context)
1880 if (mapping_ != nullptr)
1882 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1883 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1886 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1891 const gmx::IKeyValueTreeBackMapping* mapping_;
1896 void get_ir(const char* mdparin,
1897 const char* mdparout,
1898 gmx::MDModules* mdModules,
1901 WriteMdpHeader writeMdpHeader,
1905 double dumdub[2][6];
1907 char warn_buf[STRLEN];
1908 t_lambda* fep = ir->fepvals.get();
1909 t_expanded* expand = ir->expandedvals.get();
1911 const char* no_names[] = { "no", nullptr };
1913 init_inputrec_strings();
1914 gmx::TextInputFile stream(mdparin);
1915 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1917 snew(dumstr[0], STRLEN);
1918 snew(dumstr[1], STRLEN);
1920 /* ignore the following deprecated commands */
1921 replace_inp_entry(inp, "title", nullptr);
1922 replace_inp_entry(inp, "cpp", nullptr);
1923 replace_inp_entry(inp, "domain-decomposition", nullptr);
1924 replace_inp_entry(inp, "andersen-seed", nullptr);
1925 replace_inp_entry(inp, "dihre", nullptr);
1926 replace_inp_entry(inp, "dihre-fc", nullptr);
1927 replace_inp_entry(inp, "dihre-tau", nullptr);
1928 replace_inp_entry(inp, "nstdihreout", nullptr);
1929 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1930 replace_inp_entry(inp, "optimize-fft", nullptr);
1931 replace_inp_entry(inp, "adress_type", nullptr);
1932 replace_inp_entry(inp, "adress_const_wf", nullptr);
1933 replace_inp_entry(inp, "adress_ex_width", nullptr);
1934 replace_inp_entry(inp, "adress_hy_width", nullptr);
1935 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1936 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1937 replace_inp_entry(inp, "adress_site", nullptr);
1938 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1939 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1940 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1941 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1942 replace_inp_entry(inp, "rlistlong", nullptr);
1943 replace_inp_entry(inp, "nstcalclr", nullptr);
1944 replace_inp_entry(inp, "pull-print-com2", nullptr);
1945 replace_inp_entry(inp, "gb-algorithm", nullptr);
1946 replace_inp_entry(inp, "nstgbradii", nullptr);
1947 replace_inp_entry(inp, "rgbradii", nullptr);
1948 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1949 replace_inp_entry(inp, "gb-saltconc", nullptr);
1950 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1951 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1952 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1953 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1954 replace_inp_entry(inp, "sa-algorithm", nullptr);
1955 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1956 replace_inp_entry(inp, "ns-type", nullptr);
1958 /* replace the following commands with the clearer new versions*/
1959 replace_inp_entry(inp, "unconstrained-start", "continuation");
1960 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1961 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1962 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1963 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1964 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1965 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1967 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1968 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1969 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1970 setStringEntry(&inp, "include", opts->include, nullptr);
1971 printStringNoNewline(
1972 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1973 setStringEntry(&inp, "define", opts->define, nullptr);
1975 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1976 ir->eI = getEnum<IntegrationAlgorithm>(&inp, "integrator", wi);
1977 printStringNoNewline(&inp, "Start time and timestep in ps");
1978 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1979 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1980 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1981 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1982 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1983 printStringNoNewline(
1984 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1985 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1986 printStringNoNewline(&inp, "Multiple time-stepping");
1987 ir->useMts = (getEnum<Boolean>(&inp, "mts", wi) != Boolean::No);
1990 gmx::GromppMtsOpts& mtsOpts = opts->mtsOpts;
1991 mtsOpts.numLevels = get_eint(&inp, "mts-levels", 2, wi);
1992 mtsOpts.level2Forces = setStringEntry(&inp, "mts-level2-forces", "longrange-nonbonded");
1993 mtsOpts.level2Factor = get_eint(&inp, "mts-level2-factor", 2, wi);
1995 // We clear after reading without dynamics to not force the user to remove MTS mdp options
1996 if (!EI_DYNAMICS(ir->eI))
2001 printStringNoNewline(&inp, "mode for center of mass motion removal");
2002 ir->comm_mode = getEnum<ComRemovalAlgorithm>(&inp, "comm-mode", wi);
2003 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
2004 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
2005 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
2006 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
2008 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2009 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2010 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2011 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2014 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2015 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2016 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2017 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2018 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2019 ir->niter = get_eint(&inp, "niter", 20, wi);
2020 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2021 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2022 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2023 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2024 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2026 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2027 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2029 /* Output options */
2030 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2031 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2032 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2033 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2034 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2035 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2036 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2037 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2038 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2039 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2040 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2041 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2042 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2043 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2044 printStringNoNewline(&inp, "default, all atoms will be written.");
2045 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2046 printStringNoNewline(&inp, "Selection of energy groups");
2047 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2049 /* Neighbor searching */
2050 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2051 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2052 ir->cutoff_scheme = getEnum<CutoffScheme>(&inp, "cutoff-scheme", wi);
2053 printStringNoNewline(&inp, "nblist update frequency");
2054 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2055 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2056 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2057 std::vector<const char*> pbcTypesNamesChar;
2058 for (const auto& pbcTypeName : c_pbcTypeNames)
2060 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2062 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2063 ir->bPeriodicMols = getEnum<Boolean>(&inp, "periodic-molecules", wi) != Boolean::No;
2064 printStringNoNewline(&inp,
2065 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2066 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2067 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2068 printStringNoNewline(&inp, "nblist cut-off");
2069 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2070 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2072 /* Electrostatics */
2073 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2074 printStringNoNewline(&inp, "Method for doing electrostatics");
2075 ir->coulombtype = getEnum<CoulombInteractionType>(&inp, "coulombtype", wi);
2076 ir->coulomb_modifier = getEnum<InteractionModifiers>(&inp, "coulomb-modifier", wi);
2077 printStringNoNewline(&inp, "cut-off lengths");
2078 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2079 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2080 printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
2081 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2082 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2083 printStringNoNewline(&inp, "Method for doing Van der Waals");
2084 ir->vdwtype = getEnum<VanDerWaalsType>(&inp, "vdw-type", wi);
2085 ir->vdw_modifier = getEnum<InteractionModifiers>(&inp, "vdw-modifier", wi);
2086 printStringNoNewline(&inp, "cut-off lengths");
2087 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2088 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2089 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2090 ir->eDispCorr = getEnum<DispersionCorrectionType>(&inp, "DispCorr", wi);
2091 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2092 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2093 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2094 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2095 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2096 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2097 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2098 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2099 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2100 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2101 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2102 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2103 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2104 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2105 ir->ljpme_combination_rule = getEnum<LongRangeVdW>(&inp, "lj-pme-comb-rule", wi);
2106 ir->ewald_geometry = getEnum<EwaldGeometry>(&inp, "ewald-geometry", wi);
2107 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2109 /* Implicit solvation is no longer supported, but we need grompp
2110 to be able to refuse old .mdp files that would have built a tpr
2111 to run it. Thus, only "no" is accepted. */
2112 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2114 /* Coupling stuff */
2115 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2116 printStringNoNewline(&inp, "Temperature coupling");
2117 ir->etc = getEnum<TemperatureCoupling>(&inp, "tcoupl", wi);
2118 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2119 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2120 ir->bPrintNHChains = (getEnum<Boolean>(&inp, "print-nose-hoover-chain-variables", wi) != Boolean::No);
2121 printStringNoNewline(&inp, "Groups to couple separately");
2122 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2123 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2124 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2125 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2126 printStringNoNewline(&inp, "pressure coupling");
2127 ir->epc = getEnum<PressureCoupling>(&inp, "pcoupl", wi);
2128 ir->epct = getEnum<PressureCouplingType>(&inp, "pcoupltype", wi);
2129 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2130 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2131 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2132 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2133 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2134 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2135 ir->refcoord_scaling = getEnum<RefCoordScaling>(&inp, "refcoord-scaling", wi);
2138 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2139 ir->bQMMM = (getEnum<Boolean>(&inp, "QMMM", wi) != Boolean::No);
2140 printStringNoNewline(&inp, "Groups treated with MiMiC");
2141 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2143 /* Simulated annealing */
2144 printStringNewline(&inp, "SIMULATED ANNEALING");
2145 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2146 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2147 printStringNoNewline(&inp,
2148 "Number of time points to use for specifying annealing in each group");
2149 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2150 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2151 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2152 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2153 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2156 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2157 opts->bGenVel = (getEnum<Boolean>(&inp, "gen-vel", wi) != Boolean::No);
2158 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2159 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2162 printStringNewline(&inp, "OPTIONS FOR BONDS");
2163 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2164 printStringNoNewline(&inp, "Type of constraint algorithm");
2165 ir->eConstrAlg = getEnum<ConstraintAlgorithm>(&inp, "constraint-algorithm", wi);
2166 printStringNoNewline(&inp, "Do not constrain the start configuration");
2167 ir->bContinuation = (getEnum<Boolean>(&inp, "continuation", wi) != Boolean::No);
2168 printStringNoNewline(&inp,
2169 "Use successive overrelaxation to reduce the number of shake iterations");
2170 ir->bShakeSOR = (getEnum<Boolean>(&inp, "Shake-SOR", wi) != Boolean::No);
2171 printStringNoNewline(&inp, "Relative tolerance of shake");
2172 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2173 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2174 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2175 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2176 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2177 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2178 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2179 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2180 printStringNoNewline(&inp, "rotates over more degrees than");
2181 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2182 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2183 opts->bMorse = (getEnum<Boolean>(&inp, "morse", wi) != Boolean::No);
2185 /* Energy group exclusions */
2186 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2187 printStringNoNewline(
2188 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2189 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2192 printStringNewline(&inp, "WALLS");
2193 printStringNoNewline(
2194 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2195 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2196 ir->wall_type = getEnum<WallType>(&inp, "wall-type", wi);
2197 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2198 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2199 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2200 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2203 printStringNewline(&inp, "COM PULLING");
2204 ir->bPull = (getEnum<Boolean>(&inp, "pull", wi) != Boolean::No);
2207 ir->pull = std::make_unique<pull_params_t>();
2208 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
2212 for (int c = 0; c < ir->pull->ncoord; c++)
2214 if (ir->pull->coord[c].eType == PullingAlgorithm::Constraint)
2217 "Constraint COM pulling is not supported in combination with "
2218 "multiple time stepping");
2226 NOTE: needs COM pulling or free energy input */
2227 printStringNewline(&inp, "AWH biasing");
2228 ir->bDoAwh = (getEnum<Boolean>(&inp, "awh", wi) != Boolean::No);
2231 ir->awhParams = std::make_unique<gmx::AwhParams>(&inp, wi);
2234 /* Enforced rotation */
2235 printStringNewline(&inp, "ENFORCED ROTATION");
2236 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2237 ir->bRot = (getEnum<Boolean>(&inp, "rotation", wi) != Boolean::No);
2241 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2244 /* Interactive MD */
2246 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2247 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2248 if (inputrecStrings->imd_grp[0] != '\0')
2255 printStringNewline(&inp, "NMR refinement stuff");
2256 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2257 ir->eDisre = getEnum<DistanceRestraintRefinement>(&inp, "disre", wi);
2258 printStringNoNewline(
2259 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2260 ir->eDisreWeighting = getEnum<DistanceRestraintWeighting>(&inp, "disre-weighting", wi);
2261 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2262 ir->bDisreMixed = (getEnum<Boolean>(&inp, "disre-mixed", wi) != Boolean::No);
2263 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2264 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2265 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2266 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2267 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2268 opts->bOrire = (getEnum<Boolean>(&inp, "orire", wi) != Boolean::No);
2269 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2270 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2271 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2272 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2273 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2274 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2276 /* free energy variables */
2277 printStringNewline(&inp, "Free energy variables");
2278 ir->efep = getEnum<FreeEnergyPerturbationType>(&inp, "free-energy", wi);
2279 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2280 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2281 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2282 opts->bCoupleIntra = (getEnum<Boolean>(&inp, "couple-intramol", wi) != Boolean::No);
2284 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2286 it was not entered */
2287 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2288 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2289 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2290 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Fep] =
2291 setStringEntry(&inp, "fep-lambdas", "");
2292 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Mass] =
2293 setStringEntry(&inp, "mass-lambdas", "");
2294 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Coul] =
2295 setStringEntry(&inp, "coul-lambdas", "");
2296 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Vdw] =
2297 setStringEntry(&inp, "vdw-lambdas", "");
2298 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Bonded] =
2299 setStringEntry(&inp, "bonded-lambdas", "");
2300 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Restraint] =
2301 setStringEntry(&inp, "restraint-lambdas", "");
2302 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Temperature] =
2303 setStringEntry(&inp, "temperature-lambdas", "");
2304 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2305 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2306 fep->edHdLPrintEnergy = getEnum<FreeEnergyPrintEnergy>(&inp, "dhdl-print-energy", wi);
2307 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2308 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2309 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2310 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2311 fep->bScCoul = (getEnum<Boolean>(&inp, "sc-coul", wi) != Boolean::No);
2312 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2313 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2314 fep->separate_dhdl_file = getEnum<SeparateDhdlFile>(&inp, "separate-dhdl-file", wi);
2315 fep->dhdl_derivatives = getEnum<DhDlDerivativeCalculation>(&inp, "dhdl-derivatives", wi);
2316 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2317 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2319 /* Non-equilibrium MD stuff */
2320 printStringNewline(&inp, "Non-equilibrium MD stuff");
2321 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2322 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2323 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2324 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2326 /* simulated tempering variables */
2327 printStringNewline(&inp, "simulated tempering variables");
2328 ir->bSimTemp = (getEnum<Boolean>(&inp, "simulated-tempering", wi) != Boolean::No);
2329 ir->simtempvals->eSimTempScale = getEnum<SimulatedTempering>(&inp, "simulated-tempering-scaling", wi);
2330 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2331 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2333 /* expanded ensemble variables */
2334 if (ir->efep == FreeEnergyPerturbationType::Expanded || ir->bSimTemp)
2336 read_expandedparams(&inp, expand, wi);
2339 /* Electric fields */
2341 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2342 gmx::KeyValueTreeTransformer transform;
2343 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2344 mdModules->initMdpTransform(transform.rules());
2345 for (const auto& path : transform.mappedPaths())
2347 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2348 mark_einp_set(inp, path[0].c_str());
2350 MdpErrorHandler errorHandler(wi);
2351 auto result = transform.transform(convertedValues, &errorHandler);
2352 ir->params = new gmx::KeyValueTreeObject(result.object());
2353 mdModules->adjustInputrecBasedOnModules(ir);
2354 errorHandler.setBackMapping(result.backMapping());
2355 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2358 /* Ion/water position swapping ("computational electrophysiology") */
2359 printStringNewline(&inp,
2360 "Ion/water position swapping for computational electrophysiology setups");
2361 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2362 ir->eSwapCoords = getEnum<SwapType>(&inp, "swapcoords", wi);
2363 if (ir->eSwapCoords != SwapType::No)
2370 printStringNoNewline(&inp, "Swap attempt frequency");
2371 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2372 printStringNoNewline(&inp, "Number of ion types to be controlled");
2373 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2376 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2378 ir->swap->ngrp = nIonTypes + static_cast<int>(SwapGroupSplittingType::Count);
2379 snew(ir->swap->grp, ir->swap->ngrp);
2380 for (i = 0; i < ir->swap->ngrp; i++)
2382 snew(ir->swap->grp[i].molname, STRLEN);
2384 printStringNoNewline(&inp,
2385 "Two index groups that contain the compartment-partitioning atoms");
2386 setStringEntry(&inp,
2388 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
2390 setStringEntry(&inp,
2392 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname,
2394 printStringNoNewline(&inp,
2395 "Use center of mass of split groups (yes/no), otherwise center of "
2396 "geometry is used");
2397 ir->swap->massw_split[0] = (getEnum<Boolean>(&inp, "massw-split0", wi) != Boolean::No);
2398 ir->swap->massw_split[1] = (getEnum<Boolean>(&inp, "massw-split1", wi) != Boolean::No);
2400 printStringNoNewline(&inp, "Name of solvent molecules");
2401 setStringEntry(&inp,
2403 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].molname,
2406 printStringNoNewline(&inp,
2407 "Split cylinder: radius, upper and lower extension (nm) (this will "
2408 "define the channels)");
2409 printStringNoNewline(&inp,
2410 "Note that the split cylinder settings do not have an influence on "
2411 "the swapping protocol,");
2412 printStringNoNewline(
2414 "however, if correctly defined, the permeation events are recorded per channel");
2415 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2416 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2417 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2418 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2419 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2420 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2422 printStringNoNewline(
2424 "Average the number of ions per compartment over these many swap attempt steps");
2425 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2427 printStringNoNewline(
2428 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2429 printStringNoNewline(
2430 &inp, "and the requested number of ions of this type in compartments A and B");
2431 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2432 for (i = 0; i < nIonTypes; i++)
2434 int ig = static_cast<int>(SwapGroupSplittingType::Count) + i;
2436 sprintf(buf, "iontype%d-name", i);
2437 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2438 sprintf(buf, "iontype%d-in-A", i);
2439 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2440 sprintf(buf, "iontype%d-in-B", i);
2441 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2444 printStringNoNewline(
2446 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2447 printStringNoNewline(
2449 "at maximum distance (= bulk concentration) to the split group layers. However,");
2450 printStringNoNewline(&inp,
2451 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2452 "layer from the middle at 0.0");
2453 printStringNoNewline(&inp,
2454 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2455 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2456 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2457 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2458 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2460 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2463 printStringNoNewline(
2464 &inp, "Start to swap ions if threshold difference to requested count is reached");
2465 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2468 /* AdResS is no longer supported, but we need grompp to be able to
2469 refuse to process old .mdp files that used it. */
2470 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2472 /* User defined thingies */
2473 printStringNewline(&inp, "User defined thingies");
2474 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2475 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2476 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2477 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2478 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2479 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2480 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2481 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2482 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2483 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2487 gmx::TextOutputFile stream(mdparout);
2488 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2490 // Transform module data into a flat key-value tree for output.
2491 gmx::KeyValueTreeBuilder builder;
2492 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2493 mdModules->buildMdpOutput(&builderObject);
2495 gmx::TextWriter writer(&stream);
2496 writeKeyValueTreeAsMdp(&writer, builder.build());
2501 /* Process options if necessary */
2502 for (m = 0; m < 2; m++)
2504 for (i = 0; i < 2 * DIM; i++)
2508 if (ir->epc != PressureCoupling::No)
2512 case PressureCouplingType::Isotropic:
2513 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2517 "Pressure coupling incorrect number of values (I need exactly 1)");
2519 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2521 case PressureCouplingType::SemiIsotropic:
2522 case PressureCouplingType::SurfaceTension:
2523 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2527 "Pressure coupling incorrect number of values (I need exactly 2)");
2529 dumdub[m][YY] = dumdub[m][XX];
2531 case PressureCouplingType::Anisotropic:
2532 if (sscanf(dumstr[m],
2533 "%lf%lf%lf%lf%lf%lf",
2544 "Pressure coupling incorrect number of values (I need exactly 6)");
2549 "Pressure coupling type %s not implemented yet",
2550 enumValueToString(ir->epct));
2554 clear_mat(ir->ref_p);
2555 clear_mat(ir->compress);
2556 for (i = 0; i < DIM; i++)
2558 ir->ref_p[i][i] = dumdub[1][i];
2559 ir->compress[i][i] = dumdub[0][i];
2561 if (ir->epct == PressureCouplingType::Anisotropic)
2563 ir->ref_p[XX][YY] = dumdub[1][3];
2564 ir->ref_p[XX][ZZ] = dumdub[1][4];
2565 ir->ref_p[YY][ZZ] = dumdub[1][5];
2566 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2569 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2570 "apply a threefold shear stress?\n");
2572 ir->compress[XX][YY] = dumdub[0][3];
2573 ir->compress[XX][ZZ] = dumdub[0][4];
2574 ir->compress[YY][ZZ] = dumdub[0][5];
2575 for (i = 0; i < DIM; i++)
2577 for (m = 0; m < i; m++)
2579 ir->ref_p[i][m] = ir->ref_p[m][i];
2580 ir->compress[i][m] = ir->compress[m][i];
2585 if (ir->comm_mode == ComRemovalAlgorithm::No)
2590 opts->couple_moltype = nullptr;
2591 if (strlen(inputrecStrings->couple_moltype) > 0)
2593 if (ir->efep != FreeEnergyPerturbationType::No)
2595 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2596 if (opts->couple_lam0 == opts->couple_lam1)
2598 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2600 if (ir->eI == IntegrationAlgorithm::MD
2601 && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2605 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2612 "Free energy is turned off, so we will not decouple the molecule listed "
2616 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2617 if (ir->efep != FreeEnergyPerturbationType::No)
2619 if (fep->delta_lambda != 0)
2621 ir->efep = FreeEnergyPerturbationType::SlowGrowth;
2625 if (fep->edHdLPrintEnergy == FreeEnergyPrintEnergy::Yes)
2627 fep->edHdLPrintEnergy = FreeEnergyPrintEnergy::Total;
2629 "Old option for dhdl-print-energy given: "
2630 "changing \"yes\" to \"total\"\n");
2633 if (ir->bSimTemp && (fep->edHdLPrintEnergy == FreeEnergyPrintEnergy::No))
2635 /* always print out the energy to dhdl if we are doing
2636 expanded ensemble, since we need the total energy for
2637 analysis if the temperature is changing. In some
2638 conditions one may only want the potential energy, so
2639 we will allow that if the appropriate mdp setting has
2640 been enabled. Otherwise, total it is:
2642 fep->edHdLPrintEnergy = FreeEnergyPrintEnergy::Total;
2645 if ((ir->efep != FreeEnergyPerturbationType::No) || ir->bSimTemp)
2647 ir->bExpanded = FALSE;
2648 if ((ir->efep == FreeEnergyPerturbationType::Expanded) || ir->bSimTemp)
2650 ir->bExpanded = TRUE;
2652 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2653 if (ir->bSimTemp) /* done after fep params */
2655 do_simtemp_params(ir);
2658 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2659 * setup and not on the old way of specifying the free-energy setup,
2660 * we should check for using soft-core when not needed, since that
2661 * can complicate the sampling significantly.
2662 * Note that we only check for the automated coupling setup.
2663 * If the (advanced) user does FEP through manual topology changes,
2664 * this check will not be triggered.
2666 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->n_lambda == 0
2667 && ir->fepvals->sc_alpha != 0
2668 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2671 "You are using soft-core interactions while the Van der Waals interactions are "
2672 "not decoupled (note that the sc-coul option is only active when using lambda "
2673 "states). Although this will not lead to errors, you will need much more "
2674 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2679 ir->fepvals->n_lambda = 0;
2682 /* WALL PARAMETERS */
2684 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2686 /* ORIENTATION RESTRAINT PARAMETERS */
2688 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2690 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2693 /* DEFORMATION PARAMETERS */
2695 clear_mat(ir->deform);
2696 for (i = 0; i < 6; i++)
2701 double gmx_unused canary;
2702 int ndeform = sscanf(inputrecStrings->deform,
2703 "%lf %lf %lf %lf %lf %lf %lf",
2712 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2716 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2717 inputrecStrings->deform)
2720 for (i = 0; i < 3; i++)
2722 ir->deform[i][i] = dumdub[0][i];
2724 ir->deform[YY][XX] = dumdub[0][3];
2725 ir->deform[ZZ][XX] = dumdub[0][4];
2726 ir->deform[ZZ][YY] = dumdub[0][5];
2727 if (ir->epc != PressureCoupling::No)
2729 for (i = 0; i < 3; i++)
2731 for (j = 0; j <= i; j++)
2733 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2735 warning_error(wi, "A box element has deform set and compressibility > 0");
2739 for (i = 0; i < 3; i++)
2741 for (j = 0; j < i; j++)
2743 if (ir->deform[i][j] != 0)
2745 for (m = j; m < DIM; m++)
2747 if (ir->compress[m][j] != 0)
2750 "An off-diagonal box element has deform set while "
2751 "compressibility > 0 for the same component of another box "
2752 "vector, this might lead to spurious periodicity effects.");
2753 warning(wi, warn_buf);
2761 /* Ion/water position swapping checks */
2762 if (ir->eSwapCoords != SwapType::No)
2764 if (ir->swap->nstswap < 1)
2766 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2768 if (ir->swap->nAverage < 1)
2770 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2772 if (ir->swap->threshold < 1.0)
2774 warning_error(wi, "Ion count threshold must be at least 1.\n");
2778 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2781 std::vector<std::string> errorMessages;
2782 ir->mtsLevels = gmx::setupMtsLevels(opts->mtsOpts, &errorMessages);
2784 for (const auto& errorMessage : errorMessages)
2786 warning_error(wi, errorMessage.c_str());
2792 gmx::checkAwhParams(*ir->awhParams, *ir, wi);
2799 /* We would like gn to be const as well, but C doesn't allow this */
2800 /* TODO this is utility functionality (search for the index of a
2801 string in a collection), so should be refactored and located more
2803 int search_string(const char* s, int ng, char* gn[])
2807 for (i = 0; (i < ng); i++)
2809 if (gmx_strcasecmp(s, gn[i]) == 0)
2816 "Group %s referenced in the .mdp file was not found in the index file.\n"
2817 "Group names must match either [moleculetype] names or custom index group\n"
2818 "names, in which case you must supply an index file to the '-n' option\n"
2823 static void atomGroupRangeValidation(int natoms, int groupIndex, const t_blocka& block)
2825 /* Now go over the atoms in the group */
2826 for (int j = block.index[groupIndex]; (j < block.index[groupIndex + 1]); j++)
2828 int aj = block.a[j];
2830 /* Range checking */
2831 if ((aj < 0) || (aj >= natoms))
2833 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2838 static void do_numbering(int natoms,
2839 SimulationGroups* groups,
2840 gmx::ArrayRef<std::string> groupsFromMdpFile,
2843 SimulationAtomGroupType gtype,
2849 unsigned short* cbuf;
2850 AtomGroupIndices* grps = &(groups->groups[gtype]);
2853 char warn_buf[STRLEN];
2855 title = shortName(gtype);
2858 /* Mark all id's as not set */
2859 for (int i = 0; (i < natoms); i++)
2864 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2866 /* Lookup the group name in the block structure */
2867 const int gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2868 if ((grptp != egrptpONE) || (i == 0))
2870 grps->emplace_back(gid);
2872 GMX_ASSERT(block, "Can't have a nullptr block");
2873 atomGroupRangeValidation(natoms, gid, *block);
2874 /* Now go over the atoms in the group */
2875 for (int j = block->index[gid]; (j < block->index[gid + 1]); j++)
2877 const int aj = block->a[j];
2878 /* Lookup up the old group number */
2879 const int ognr = cbuf[aj];
2882 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title, ognr + 1, i + 1);
2886 /* Store the group number in buffer */
2887 if (grptp == egrptpONE)
2900 /* Now check whether we have done all atoms */
2903 if (grptp == egrptpALL)
2905 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2907 else if (grptp == egrptpPART)
2909 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2910 warning_note(wi, warn_buf);
2912 /* Assign all atoms currently unassigned to a rest group */
2913 for (int j = 0; (j < natoms); j++)
2915 if (cbuf[j] == NOGID)
2917 cbuf[j] = grps->size();
2920 if (grptp != egrptpPART)
2924 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title, natoms - ntot);
2926 /* Add group name "rest" */
2927 grps->emplace_back(restnm);
2929 /* Assign the rest name to all atoms not currently assigned to a group */
2930 for (int j = 0; (j < natoms); j++)
2932 if (cbuf[j] == NOGID)
2934 // group size was not updated before this here, so need to use -1.
2935 cbuf[j] = grps->size() - 1;
2941 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2943 /* All atoms are part of one (or no) group, no index required */
2944 groups->groupNumbers[gtype].clear();
2948 for (int j = 0; (j < natoms); j++)
2950 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2957 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2960 pull_params_t* pull;
2961 int natoms, imin, jmin;
2962 int * nrdf2, *na_vcm, na_tot;
2963 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2968 * First calc 3xnr-atoms for each group
2969 * then subtract half a degree of freedom for each constraint
2971 * Only atoms and nuclei contribute to the degrees of freedom...
2976 const SimulationGroups& groups = mtop->groups;
2977 natoms = mtop->natoms;
2979 /* Allocate one more for a possible rest group */
2980 /* We need to sum degrees of freedom into doubles,
2981 * since floats give too low nrdf's above 3 million atoms.
2983 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
2984 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2985 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2986 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2987 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2989 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2993 for (gmx::index i = 0;
2994 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
2998 clear_ivec(dof_vcm[i]);
3000 nrdf_vcm_sub[i] = 0;
3002 snew(nrdf2, natoms);
3003 for (const AtomProxy atomP : AtomRange(*mtop))
3005 const t_atom& local = atomP.atom();
3006 int i = atomP.globalAtomNumber();
3008 if (local.ptype == ParticleType::Atom || local.ptype == ParticleType::Nucleus)
3010 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
3011 for (int d = 0; d < DIM; d++)
3013 if (opts->nFreeze[g][d] == 0)
3015 /* Add one DOF for particle i (counted as 2*1) */
3017 /* VCM group i has dim d as a DOF */
3018 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
3022 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
3024 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
3030 for (const gmx_molblock_t& molb : mtop->molblock)
3032 const gmx_moltype_t& molt = mtop->moltype[molb.type];
3033 const t_atom* atom = molt.atoms.atom;
3034 for (int mol = 0; mol < molb.nmol; mol++)
3036 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
3038 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
3039 for (int i = 0; i < molt.ilist[ftype].size();)
3041 /* Subtract degrees of freedom for the constraints,
3042 * if the particles still have degrees of freedom left.
3043 * If one of the particles is a vsite or a shell, then all
3044 * constraint motion will go there, but since they do not
3045 * contribute to the constraints the degrees of freedom do not
3048 int ai = as + ia[i + 1];
3049 int aj = as + ia[i + 2];
3050 if (((atom[ia[i + 1]].ptype == ParticleType::Nucleus)
3051 || (atom[ia[i + 1]].ptype == ParticleType::Atom))
3052 && ((atom[ia[i + 2]].ptype == ParticleType::Nucleus)
3053 || (atom[ia[i + 2]].ptype == ParticleType::Atom)))
3071 imin = std::min(imin, nrdf2[ai]);
3072 jmin = std::min(jmin, nrdf2[aj]);
3075 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3077 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3079 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3081 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3084 i += interaction_function[ftype].nratoms + 1;
3087 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3088 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3090 /* Subtract 1 dof from every atom in the SETTLE */
3091 for (int j = 0; j < 3; j++)
3093 int ai = as + ia[i + 1 + j];
3094 imin = std::min(2, nrdf2[ai]);
3096 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3098 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3103 as += molt.atoms.nr;
3109 /* Correct nrdf for the COM constraints.
3110 * We correct using the TC and VCM group of the first atom
3111 * in the reference and pull group. If atoms in one pull group
3112 * belong to different TC or VCM groups it is anyhow difficult
3113 * to determine the optimal nrdf assignment.
3115 pull = ir->pull.get();
3117 for (int i = 0; i < pull->ncoord; i++)
3119 if (pull->coord[i].eType != PullingAlgorithm::Constraint)
3126 for (int j = 0; j < 2; j++)
3128 const t_pull_group* pgrp;
3130 pgrp = &pull->group[pull->coord[i].group[j]];
3132 if (!pgrp->ind.empty())
3134 /* Subtract 1/2 dof from each group */
3135 int ai = pgrp->ind[0];
3136 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3138 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3140 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3143 "Center of mass pulling constraints caused the number of degrees "
3144 "of freedom for temperature coupling group %s to be negative",
3145 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3146 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3151 /* We need to subtract the whole DOF from group j=1 */
3158 if (ir->nstcomm != 0)
3160 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3161 "Expect at least one group when removing COM motion");
3163 /* We remove COM motion up to dim ndof_com() */
3164 const int ndim_rm_vcm = ndof_com(ir);
3166 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3167 * the number of degrees of freedom in each vcm group when COM
3168 * translation is removed and 6 when rotation is removed as well.
3169 * Note that we do not and should not include the rest group here.
3171 for (gmx::index j = 0;
3172 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]);
3175 switch (ir->comm_mode)
3177 case ComRemovalAlgorithm::Linear:
3178 case ComRemovalAlgorithm::LinearAccelerationCorrection:
3179 nrdf_vcm_sub[j] = 0;
3180 for (int d = 0; d < ndim_rm_vcm; d++)
3188 case ComRemovalAlgorithm::Angular: nrdf_vcm_sub[j] = 6; break;
3189 default: gmx_incons("Checking comm_mode");
3193 for (gmx::index i = 0;
3194 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]);
3197 /* Count the number of atoms of TC group i for every VCM group */
3198 for (gmx::index j = 0;
3199 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3205 for (int ai = 0; ai < natoms; ai++)
3207 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3209 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3213 /* Correct for VCM removal according to the fraction of each VCM
3214 * group present in this TC group.
3216 nrdf_uc = nrdf_tc[i];
3218 for (gmx::index j = 0;
3219 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3222 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3224 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3225 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3230 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3232 opts->nrdf[i] = nrdf_tc[i];
3233 if (opts->nrdf[i] < 0)
3238 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3239 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]],
3248 sfree(nrdf_vcm_sub);
3251 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3253 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3254 * But since this is much larger than STRLEN, such a line can not be parsed.
3255 * The real maximum is the number of names that fit in a string: STRLEN/2.
3257 #define EGP_MAX (STRLEN / 2)
3261 auto names = gmx::splitString(val);
3262 if (names.size() % 2 != 0)
3264 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3266 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3268 for (size_t i = 0; i < names.size() / 2; i++)
3270 // TODO this needs to be replaced by a solution using std::find_if
3274 names[2 * i].c_str(),
3275 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3281 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3286 names[2 * i + 1].c_str(),
3287 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3293 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3295 if ((j < nr) && (k < nr))
3297 ir->opts.egp_flags[nr * j + k] |= flag;
3298 ir->opts.egp_flags[nr * k + j] |= flag;
3307 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3309 int ig = -1, i = 0, gind;
3313 /* Just a quick check here, more thorough checks are in mdrun */
3314 if (strcmp(swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
3315 swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname)
3319 "The split groups can not both be '%s'.",
3320 swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname);
3323 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3324 for (ig = 0; ig < swap->ngrp; ig++)
3326 swapg = &swap->grp[ig];
3327 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3328 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3333 "%s group '%s' contains %d atoms.\n",
3334 ig < 3 ? enumValueToString(static_cast<SwapGroupSplittingType>(ig)) : "Swap",
3335 swap->grp[ig].molname,
3337 snew(swapg->ind, swapg->nat);
3338 for (i = 0; i < swapg->nat; i++)
3340 swapg->ind[i] = grps->a[grps->index[gind] + i];
3345 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3351 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3356 ig = search_string(IMDgname, grps->nr, gnames);
3357 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3359 if (IMDgroup->nat > 0)
3362 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3366 snew(IMDgroup->ind, IMDgroup->nat);
3367 for (i = 0; i < IMDgroup->nat; i++)
3369 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3374 /* Checks whether atoms are both part of a COM removal group and frozen.
3375 * If a fully frozen atom is part of a COM removal group, it is removed
3376 * from the COM removal group. A note is issued if such atoms are present.
3377 * A warning is issued for atom with one or two dimensions frozen that
3378 * are part of a COM removal group (mdrun would need to compute COM mass
3379 * per dimension to handle this correctly).
3380 * Also issues a warning when non-frozen atoms are not part of a COM
3381 * removal group while COM removal is active.
3383 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3385 const t_grpopts& opts,
3388 const int vcmRestGroup =
3389 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3391 int numFullyFrozenVcmAtoms = 0;
3392 int numPartiallyFrozenVcmAtoms = 0;
3393 int numNonVcmAtoms = 0;
3394 for (int a = 0; a < numAtoms; a++)
3396 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3397 int numFrozenDims = 0;
3398 for (int d = 0; d < DIM; d++)
3400 numFrozenDims += opts.nFreeze[freezeGroup][d];
3403 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3404 if (vcmGroup < vcmRestGroup)
3406 if (numFrozenDims == DIM)
3408 /* Do not remove COM motion for this fully frozen atom */
3409 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3411 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3414 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3415 numFullyFrozenVcmAtoms++;
3417 else if (numFrozenDims > 0)
3419 numPartiallyFrozenVcmAtoms++;
3422 else if (numFrozenDims < DIM)
3428 if (numFullyFrozenVcmAtoms > 0)
3430 std::string warningText = gmx::formatString(
3431 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3432 "removing these atoms from the COMM removal group(s)",
3433 numFullyFrozenVcmAtoms);
3434 warning_note(wi, warningText.c_str());
3436 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3438 std::string warningText = gmx::formatString(
3439 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3440 "removal group(s), due to limitations in the code these still contribute to the "
3441 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3443 numPartiallyFrozenVcmAtoms,
3445 warning(wi, warningText.c_str());
3447 if (numNonVcmAtoms > 0)
3449 std::string warningText = gmx::formatString(
3450 "%d atoms are not part of any center of mass motion removal group.\n"
3451 "This may lead to artifacts.\n"
3452 "In most cases one should use one group for the whole system.",
3454 warning(wi, warningText.c_str());
3458 void do_index(const char* mdparin,
3462 const gmx::MDModulesNotifiers& mdModulesNotifiers,
3466 t_blocka* defaultIndexGroups;
3474 int i, j, k, restnm;
3475 bool bExcl, bTable, bAnneal;
3476 char warn_buf[STRLEN];
3480 fprintf(stderr, "processing index file...\n");
3484 snew(defaultIndexGroups, 1);
3485 snew(defaultIndexGroups->index, 1);
3487 atoms_all = gmx_mtop_global_atoms(*mtop);
3488 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3489 done_atom(&atoms_all);
3493 defaultIndexGroups = init_index(ndx, &gnames);
3496 SimulationGroups* groups = &mtop->groups;
3497 natoms = mtop->natoms;
3498 symtab = &mtop->symtab;
3500 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3502 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3504 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3505 restnm = groups->groupNames.size() - 1;
3506 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3507 srenew(gnames, defaultIndexGroups->nr + 1);
3508 gnames[restnm] = *(groups->groupNames.back());
3510 set_warning_line(wi, mdparin, -1);
3512 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3513 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3514 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3515 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3516 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3519 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3521 temperatureCouplingGroupNames.size(),
3522 temperatureCouplingReferenceValues.size(),
3523 temperatureCouplingTauValues.size());
3526 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3527 do_numbering(natoms,
3529 temperatureCouplingGroupNames,
3532 SimulationAtomGroupType::TemperatureCoupling,
3534 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST,
3537 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3539 snew(ir->opts.nrdf, nr);
3540 snew(ir->opts.tau_t, nr);
3541 snew(ir->opts.ref_t, nr);
3542 if (ir->eI == IntegrationAlgorithm::BD && ir->bd_fric == 0)
3544 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3547 if (useReferenceTemperature)
3549 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3551 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3555 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3556 for (i = 0; (i < nr); i++)
3558 if ((ir->eI == IntegrationAlgorithm::BD) && ir->opts.tau_t[i] <= 0)
3561 "With integrator %s tau-t should be larger than 0",
3562 enumValueToString(ir->eI));
3563 warning_error(wi, warn_buf);
3566 if (ir->etc != TemperatureCoupling::VRescale && ir->opts.tau_t[i] == 0)
3570 "tau-t = -1 is the value to signal that a group should not have "
3571 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3574 if (ir->opts.tau_t[i] >= 0)
3576 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3579 if (ir->etc != TemperatureCoupling::No && ir->nsttcouple == -1)
3581 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3586 if ((ir->etc == TemperatureCoupling::NoseHoover) && (ir->epc == PressureCoupling::Berendsen))
3589 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3590 "md-vv; use either vrescale temperature with berendsen pressure or "
3591 "Nose-Hoover temperature with MTTK pressure");
3593 if (ir->epc == PressureCoupling::Mttk)
3595 if (ir->etc != TemperatureCoupling::NoseHoover)
3598 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3603 if (ir->nstpcouple != ir->nsttcouple)
3605 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3606 ir->nstpcouple = ir->nsttcouple = mincouple;
3608 "for current Trotter decomposition methods with vv, nsttcouple and "
3609 "nstpcouple must be equal. Both have been reset to "
3610 "min(nsttcouple,nstpcouple) = %d",
3612 warning_note(wi, warn_buf);
3617 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3618 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3620 if (ETC_ANDERSEN(ir->etc))
3622 if (ir->nsttcouple != 1)
3626 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3627 "need for larger nsttcouple > 1, since no global parameters are computed. "
3628 "nsttcouple has been reset to 1");
3629 warning_note(wi, warn_buf);
3632 nstcmin = tcouple_min_integration_steps(ir->etc);
3635 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3638 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3639 "least %d times larger than nsttcouple*dt (%g)",
3640 enumValueToString(ir->etc),
3643 ir->nsttcouple * ir->delta_t);
3644 warning(wi, warn_buf);
3647 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3648 for (i = 0; (i < nr); i++)
3650 if (ir->opts.ref_t[i] < 0)
3652 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3655 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3656 if we are in this conditional) if mc_temp is negative */
3657 if (ir->expandedvals->mc_temp < 0)
3659 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3663 /* Simulated annealing for each group. There are nr groups */
3664 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3665 if (simulatedAnnealingGroupNames.size() == 1
3666 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3668 simulatedAnnealingGroupNames.resize(0);
3670 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3673 "Wrong number of annealing values: %zu (for %d groups)\n",
3674 simulatedAnnealingGroupNames.size(),
3679 snew(ir->opts.annealing, nr);
3680 snew(ir->opts.anneal_npoints, nr);
3681 snew(ir->opts.anneal_time, nr);
3682 snew(ir->opts.anneal_temp, nr);
3683 for (i = 0; i < nr; i++)
3685 ir->opts.annealing[i] = SimulatedAnnealing::No;
3686 ir->opts.anneal_npoints[i] = 0;
3687 ir->opts.anneal_time[i] = nullptr;
3688 ir->opts.anneal_temp[i] = nullptr;
3690 if (!simulatedAnnealingGroupNames.empty())
3693 for (i = 0; i < nr; i++)
3695 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3697 ir->opts.annealing[i] = SimulatedAnnealing::No;
3699 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3701 ir->opts.annealing[i] = SimulatedAnnealing::Single;
3704 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3706 ir->opts.annealing[i] = SimulatedAnnealing::Periodic;
3712 /* Read the other fields too */
3713 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3714 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3717 "Found %zu annealing-npoints values for %zu groups\n",
3718 simulatedAnnealingPoints.size(),
3719 simulatedAnnealingGroupNames.size());
3721 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3722 size_t numSimulatedAnnealingFields = 0;
3723 for (i = 0; i < nr; i++)
3725 if (ir->opts.anneal_npoints[i] == 1)
3729 "Please specify at least a start and an end point for annealing\n");
3731 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3732 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3733 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3736 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3738 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3741 "Found %zu annealing-time values, wanted %zu\n",
3742 simulatedAnnealingTimes.size(),
3743 numSimulatedAnnealingFields);
3745 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3746 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3749 "Found %zu annealing-temp values, wanted %zu\n",
3750 simulatedAnnealingTemperatures.size(),
3751 numSimulatedAnnealingFields);
3754 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3755 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3756 convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
3758 simulatedAnnealingTemperatures,
3760 allSimulatedAnnealingTemperatures.data());
3761 for (i = 0, k = 0; i < nr; i++)
3763 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3765 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3766 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3769 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3771 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3777 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3780 "Annealing timepoints out of order: t=%f comes after "
3782 ir->opts.anneal_time[i][j],
3783 ir->opts.anneal_time[i][j - 1]);
3786 if (ir->opts.anneal_temp[i][j] < 0)
3789 "Found negative temperature in annealing: %f\n",
3790 ir->opts.anneal_temp[i][j]);
3795 /* Print out some summary information, to make sure we got it right */
3796 for (i = 0; i < nr; i++)
3798 if (ir->opts.annealing[i] != SimulatedAnnealing::No)
3800 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3802 "Simulated annealing for group %s: %s, %d timepoints\n",
3803 *(groups->groupNames[j]),
3804 enumValueToString(ir->opts.annealing[i]),
3805 ir->opts.anneal_npoints[i]);
3806 fprintf(stderr, "Time (ps) Temperature (K)\n");
3807 /* All terms except the last one */
3808 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3812 ir->opts.anneal_time[i][j],
3813 ir->opts.anneal_temp[i][j]);
3816 /* Finally the last one */
3817 j = ir->opts.anneal_npoints[i] - 1;
3818 if (ir->opts.annealing[i] == SimulatedAnnealing::Single)
3822 ir->opts.anneal_time[i][j],
3823 ir->opts.anneal_temp[i][j]);
3829 ir->opts.anneal_time[i][j],
3830 ir->opts.anneal_temp[i][j]);
3831 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3834 "There is a temperature jump when your annealing "
3846 for (int i = 1; i < ir->pull->ngroup; i++)
3848 const int gid = search_string(
3849 inputrecStrings->pullGroupNames[i].c_str(), defaultIndexGroups->nr, gnames);
3850 GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
3851 atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
3854 process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3856 checkPullCoords(ir->pull->group, ir->pull->coord);
3861 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3864 if (ir->eSwapCoords != SwapType::No)
3866 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3869 /* Make indices for IMD session */
3872 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3875 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3876 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3877 mdModulesNotifiers.preProcessingNotifier_.notify(defaultIndexGroupsAndNames);
3879 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3880 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3881 if (freezeDims.size() != DIM * freezeGroupNames.size())
3884 "Invalid Freezing input: %zu groups and %zu freeze values",
3885 freezeGroupNames.size(),
3888 do_numbering(natoms,
3893 SimulationAtomGroupType::Freeze,
3898 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3899 ir->opts.ngfrz = nr;
3900 snew(ir->opts.nFreeze, nr);
3901 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3903 for (j = 0; (j < DIM); j++, k++)
3905 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3906 if (!ir->opts.nFreeze[i][j])
3908 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3911 "Please use Y(ES) or N(O) for freezedim only "
3913 freezeDims[k].c_str());
3914 warning(wi, warn_buf);
3919 for (; (i < nr); i++)
3921 for (j = 0; (j < DIM); j++)
3923 ir->opts.nFreeze[i][j] = 0;
3927 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3928 do_numbering(natoms,
3933 SimulationAtomGroupType::EnergyOutput,
3938 add_wall_energrps(groups, ir->nwall, symtab);
3939 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3940 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3941 do_numbering(natoms,
3946 SimulationAtomGroupType::MassCenterVelocityRemoval,
3948 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART,
3952 if (ir->comm_mode != ComRemovalAlgorithm::No)
3954 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3957 /* Now we have filled the freeze struct, so we can calculate NRDF */
3958 calc_nrdf(mtop, ir, gnames);
3960 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3961 do_numbering(natoms,
3966 SimulationAtomGroupType::User1,
3971 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3972 do_numbering(natoms,
3977 SimulationAtomGroupType::User2,
3982 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3983 do_numbering(natoms,
3985 compressedXGroupNames,
3988 SimulationAtomGroupType::CompressedPositionOutput,
3993 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3994 do_numbering(natoms,
3996 orirefFitGroupNames,
3999 SimulationAtomGroupType::OrientationRestraintsFit,
4005 /* MiMiC QMMM input processing */
4006 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
4007 if (qmGroupNames.size() > 1)
4009 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
4011 /* group rest, if any, is always MM! */
4012 do_numbering(natoms,
4017 SimulationAtomGroupType::QuantumMechanics,
4022 ir->opts.ngQM = qmGroupNames.size();
4024 /* end of MiMiC QMMM input */
4028 for (auto group : gmx::keysOf(groups->groups))
4030 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
4031 for (const auto& entry : groups->groups[group])
4033 fprintf(stderr, " %s", *(groups->groupNames[entry]));
4035 fprintf(stderr, "\n");
4039 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
4040 snew(ir->opts.egp_flags, nr * nr);
4042 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
4043 if (bExcl && ir->cutoff_scheme == CutoffScheme::Verlet)
4045 warning_error(wi, "Energy group exclusions are currently not supported");
4047 if (bExcl && EEL_FULL(ir->coulombtype))
4049 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
4052 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
4053 if (bTable && !(ir->vdwtype == VanDerWaalsType::User)
4054 && !(ir->coulombtype == CoulombInteractionType::User)
4055 && !(ir->coulombtype == CoulombInteractionType::PmeUser)
4056 && !(ir->coulombtype == CoulombInteractionType::PmeUserSwitch))
4059 "Can only have energy group pair tables in combination with user tables for VdW "
4063 /* final check before going out of scope if simulated tempering variables
4064 * need to be set to default values.
4066 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
4068 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
4071 "the value for nstexpanded was not specified for "
4072 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
4073 "by default, but it is recommended to set it to an explicit value!",
4074 ir->expandedvals->nstexpanded));
4076 for (i = 0; (i < defaultIndexGroups->nr); i++)
4081 done_blocka(defaultIndexGroups);
4082 sfree(defaultIndexGroups);
4086 static void check_disre(const gmx_mtop_t& mtop)
4088 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
4090 const gmx_ffparams_t& ffparams = mtop.ffparams;
4093 for (int i = 0; i < ffparams.numTypes(); i++)
4095 int ftype = ffparams.functype[i];
4096 if (ftype == F_DISRES)
4098 int label = ffparams.iparams[i].disres.label;
4099 if (label == old_label)
4101 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
4110 "Found %d double distance restraint indices,\n"
4111 "probably the parameters for multiple pairs in one restraint "
4112 "are not identical\n",
4118 //! Returns whether dimensions have an absolute reference due to walls, pbc or freezing
4119 static BasicVector<bool> haveAbsoluteReference(const t_inputrec& ir)
4121 BasicVector<bool> absRef = { false, false, false };
4123 /* Check the degrees of freedom of the COM (not taking COMM removal into account) */
4124 for (int d = 0; d < DIM; d++)
4126 absRef[d] = (d >= ndof_com(&ir));
4128 /* Check for freeze groups */
4129 for (int g = 0; g < ir.opts.ngfrz; g++)
4131 for (int d = 0; d < DIM; d++)
4133 if (ir.opts.nFreeze[g][d] != 0)
4143 //! Returns whether position restraints are used for dimensions
4144 static BasicVector<bool> havePositionRestraints(const gmx_mtop_t& sys)
4146 BasicVector<bool> havePosres = { false, false, false };
4148 for (const auto ilists : IListRange(sys))
4150 const auto& posResList = ilists.list()[F_POSRES];
4151 const auto& fbPosResList = ilists.list()[F_FBPOSRES];
4152 if (ilists.nmol() > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
4154 for (int i = 0; i < posResList.size(); i += 2)
4156 const t_iparams& pr = sys.ffparams.iparams[posResList.iatoms[i]];
4157 for (int d = 0; d < DIM; d++)
4159 if (pr.posres.fcA[d] != 0)
4161 havePosres[d] = true;
4165 for (int i = 0; i < fbPosResList.size(); i += 2)
4167 /* Check for flat-bottom posres */
4168 const t_iparams& pr = sys.ffparams.iparams[fbPosResList.iatoms[i]];
4169 if (pr.fbposres.k != 0)
4171 switch (pr.fbposres.geom)
4173 case efbposresSPHERE: havePosres = { true, true, true }; break;
4174 case efbposresCYLINDERX: havePosres[YY] = havePosres[ZZ] = true; break;
4175 case efbposresCYLINDERY: havePosres[XX] = havePosres[ZZ] = true; break;
4176 case efbposresCYLINDER:
4177 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4178 case efbposresCYLINDERZ: havePosres[XX] = havePosres[YY] = true; break;
4179 case efbposresX: /* d=XX */
4180 case efbposresY: /* d=YY */
4181 case efbposresZ: /* d=ZZ */
4182 havePosres[pr.fbposres.geom - efbposresX] = true;
4186 "Invalid geometry for flat-bottom position restraint.\n"
4187 "Expected nr between 1 and %d. Found %d\n",
4199 static void check_combination_rule_differences(const gmx_mtop_t& mtop,
4201 bool* bC6ParametersWorkWithGeometricRules,
4202 bool* bC6ParametersWorkWithLBRules,
4203 bool* bLBRulesPossible)
4205 int ntypes, tpi, tpj;
4208 double c6i, c6j, c12i, c12j;
4209 double c6, c6_geometric, c6_LB;
4210 double sigmai, sigmaj, epsi, epsj;
4211 bool bCanDoLBRules, bCanDoGeometricRules;
4214 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4215 * force-field floating point parameters.
4218 ptr = getenv("GMX_LJCOMB_TOL");
4222 double gmx_unused canary;
4224 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4227 FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4232 *bC6ParametersWorkWithLBRules = TRUE;
4233 *bC6ParametersWorkWithGeometricRules = TRUE;
4234 bCanDoLBRules = TRUE;
4235 ntypes = mtop.ffparams.atnr;
4236 snew(typecount, ntypes);
4237 gmx_mtop_count_atomtypes(mtop, state, typecount);
4238 *bLBRulesPossible = TRUE;
4239 for (tpi = 0; tpi < ntypes; ++tpi)
4241 c6i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4242 c12i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4243 for (tpj = tpi; tpj < ntypes; ++tpj)
4245 c6j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4246 c12j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4247 c6 = mtop.ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4248 c6_geometric = std::sqrt(c6i * c6j);
4249 if (!gmx_numzero(c6_geometric))
4251 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4253 sigmai = gmx::sixthroot(c12i / c6i);
4254 sigmaj = gmx::sixthroot(c12j / c6j);
4255 epsi = c6i * c6i / (4.0 * c12i);
4256 epsj = c6j * c6j / (4.0 * c12j);
4257 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4261 *bLBRulesPossible = FALSE;
4262 c6_LB = c6_geometric;
4264 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4269 *bC6ParametersWorkWithLBRules = FALSE;
4272 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4274 if (!bCanDoGeometricRules)
4276 *bC6ParametersWorkWithGeometricRules = FALSE;
4283 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t& mtop, warninp_t wi)
4285 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4287 check_combination_rule_differences(
4288 mtop, 0, &bC6ParametersWorkWithGeometricRules, &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4289 if (ir->ljpme_combination_rule == LongRangeVdW::LB)
4291 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4294 "You are using arithmetic-geometric combination rules "
4295 "in LJ-PME, but your non-bonded C6 parameters do not "
4296 "follow these rules.");
4301 if (!bC6ParametersWorkWithGeometricRules)
4303 if (ir->eDispCorr != DispersionCorrectionType::No)
4306 "You are using geometric combination rules in "
4307 "LJ-PME, but your non-bonded C6 parameters do "
4308 "not follow these rules. "
4309 "This will introduce very small errors in the forces and energies in "
4310 "your simulations. Dispersion correction will correct total energy "
4311 "and/or pressure for isotropic systems, but not forces or surface "
4317 "You are using geometric combination rules in "
4318 "LJ-PME, but your non-bonded C6 parameters do "
4319 "not follow these rules. "
4320 "This will introduce very small errors in the forces and energies in "
4321 "your simulations. If your system is homogeneous, consider using "
4322 "dispersion correction "
4323 "for the total energy and pressure.");
4329 static bool allTrue(const BasicVector<bool>& boolVector)
4331 return boolVector[0] && boolVector[1] && boolVector[2];
4334 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4336 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4337 GMX_ASSERT(gmx::checkMtsRequirements(*ir).empty(), "All MTS requirements should be met here");
4339 char err_buf[STRLEN];
4342 gmx_mtop_atomloop_block_t aloopb;
4343 char warn_buf[STRLEN];
4345 set_warning_line(wi, mdparin, -1);
4347 if (allTrue(haveAbsoluteReference(*ir)) && allTrue(havePositionRestraints(*sys)))
4350 "Removing center of mass motion in the presence of position restraints might "
4351 "cause artifacts. When you are using position restraints to equilibrate a "
4352 "macro-molecule, the artifacts are usually negligible.");
4355 if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0 && ir->nstlist > 1
4356 && ((EI_MD(ir->eI) || EI_SD(ir->eI))
4357 && (ir->etc == TemperatureCoupling::VRescale || ir->etc == TemperatureCoupling::Berendsen)))
4359 /* Check if a too small Verlet buffer might potentially
4360 * cause more drift than the thermostat can couple off.
4362 /* Temperature error fraction for warning and suggestion */
4363 const real T_error_warn = 0.002;
4364 const real T_error_suggest = 0.001;
4365 /* For safety: 2 DOF per atom (typical with constraints) */
4366 const real nrdf_at = 2;
4367 real T, tau, max_T_error;
4372 for (i = 0; i < ir->opts.ngtc; i++)
4374 T = std::max(T, ir->opts.ref_t[i]);
4375 tau = std::max(tau, ir->opts.tau_t[i]);
4379 /* This is a worst case estimate of the temperature error,
4380 * assuming perfect buffer estimation and no cancelation
4381 * of errors. The factor 0.5 is because energy distributes
4382 * equally over Ekin and Epot.
4384 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * gmx::c_boltz * T);
4385 if (max_T_error > T_error_warn)
4388 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4389 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4390 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4391 "%.0e or decrease tau_t.",
4396 100 * T_error_suggest,
4397 ir->verletbuf_tol * T_error_suggest / max_T_error);
4398 warning(wi, warn_buf);
4403 if (ETC_ANDERSEN(ir->etc))
4407 for (i = 0; i < ir->opts.ngtc; i++)
4410 "all tau_t must currently be equal using Andersen temperature control, "
4411 "violated for group %d",
4413 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4415 "all tau_t must be positive using Andersen temperature control, "
4419 CHECK(ir->opts.tau_t[i] < 0);
4422 if (ir->etc == TemperatureCoupling::AndersenMassive && ir->comm_mode != ComRemovalAlgorithm::No)
4424 for (i = 0; i < ir->opts.ngtc; i++)
4426 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4428 "tau_t/delta_t for group %d for temperature control method %s must be a "
4429 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4430 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4433 enumValueToString(ir->etc),
4437 CHECK(nsteps % ir->nstcomm != 0);
4442 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != IntegrationAlgorithm::BD
4443 && ir->comm_mode == ComRemovalAlgorithm::No
4444 && !(allTrue(haveAbsoluteReference(*ir)) || allTrue(havePositionRestraints(*sys)) || ir->nsteps <= 10)
4445 && !ETC_ANDERSEN(ir->etc))
4448 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4449 "rounding errors can lead to build up of kinetic energy of the center of mass");
4452 if (ir->epc == PressureCoupling::ParrinelloRahman && ir->etc == TemperatureCoupling::NoseHoover)
4455 for (int g = 0; g < ir->opts.ngtc; g++)
4457 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4459 if (ir->tau_p < 1.9 * tau_t_max)
4461 std::string message = gmx::formatString(
4462 "With %s T-coupling and %s p-coupling, "
4463 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4464 enumValueToString(ir->etc),
4465 enumValueToString(ir->epc),
4470 warning(wi, message.c_str());
4474 /* Check for pressure coupling with absolute position restraints */
4475 if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == RefCoordScaling::No)
4477 const BasicVector<bool> havePosres = havePositionRestraints(*sys);
4479 for (m = 0; m < DIM; m++)
4481 if (havePosres[m] && norm2(ir->compress[m]) > 0)
4484 "You are using pressure coupling with absolute position restraints, "
4485 "this will give artifacts. Use the refcoord_scaling option.");
4493 aloopb = gmx_mtop_atomloop_block_init(*sys);
4495 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4497 if (atom->q != 0 || atom->qB != 0)
4505 if (EEL_FULL(ir->coulombtype))
4508 "You are using full electrostatics treatment %s for a system without charges.\n"
4509 "This costs a lot of performance for just processing zeros, consider using %s "
4511 enumValueToString(ir->coulombtype),
4512 enumValueToString(CoulombInteractionType::Cut));
4513 warning(wi, err_buf);
4518 if (ir->coulombtype == CoulombInteractionType::Cut && ir->rcoulomb > 0)
4521 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4522 "You might want to consider using %s electrostatics.\n",
4523 enumValueToString(CoulombInteractionType::Pme));
4524 warning_note(wi, err_buf);
4528 /* Check if combination rules used in LJ-PME are the same as in the force field */
4529 if (EVDW_PME(ir->vdwtype))
4531 check_combination_rules(ir, *sys, wi);
4534 /* Generalized reaction field */
4535 if (ir->coulombtype == CoulombInteractionType::GRFNotused)
4538 "Generalized reaction-field electrostatics is no longer supported. "
4539 "You can use normal reaction-field instead and compute the reaction-field "
4540 "constant by hand.");
4543 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->sc_alpha != 0
4544 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4546 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4554 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4556 if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
4558 const auto absRef = haveAbsoluteReference(*ir);
4559 const auto havePosres = havePositionRestraints(*sys);
4560 for (m = 0; m < DIM; m++)
4562 if (ir->pull->coord[i].dim[m] && !(absRef[m] || havePosres[m]))
4565 "You are using an absolute reference for pulling, but the rest of "
4566 "the system does not have an absolute reference. This will lead to "
4575 for (i = 0; i < 3; i++)
4577 for (m = 0; m <= i; m++)
4579 if ((ir->epc != PressureCoupling::No && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4581 for (c = 0; c < ir->pull->ncoord; c++)
4583 if (ir->pull->coord[c].eGeom == PullGroupGeometry::DirectionPBC
4584 && ir->pull->coord[c].vec[m] != 0)
4587 "Can not have dynamic box while using pull geometry '%s' "
4589 enumValueToString(ir->pull->coord[c].eGeom),
4601 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4603 char warn_buf[STRLEN];
4606 ptr = check_box(ir->pbcType, box);
4609 warning_error(wi, ptr);
4612 if (bHasNormalConstraints && ir->eConstrAlg == ConstraintAlgorithm::Shake)
4614 if (ir->shake_tol <= 0.0)
4616 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4617 warning_error(wi, warn_buf);
4621 if ((ir->eConstrAlg == ConstraintAlgorithm::Lincs) && bHasNormalConstraints)
4623 /* If we have Lincs constraints: */
4624 if (ir->eI == IntegrationAlgorithm::MD && ir->etc == TemperatureCoupling::No
4625 && ir->eConstrAlg == ConstraintAlgorithm::Lincs && ir->nLincsIter == 1)
4628 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4629 warning_note(wi, warn_buf);
4632 if ((ir->eI == IntegrationAlgorithm::CG || ir->eI == IntegrationAlgorithm::LBFGS)
4633 && (ir->nProjOrder < 8))
4636 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4637 enumValueToString(ir->eI));
4638 warning_note(wi, warn_buf);
4640 if (ir->epc == PressureCoupling::Mttk)
4642 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4646 if (bHasAnyConstraints && ir->epc == PressureCoupling::Mttk)
4648 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4651 if (ir->LincsWarnAngle > 90.0)
4653 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4654 warning(wi, warn_buf);
4655 ir->LincsWarnAngle = 90.0;
4658 if (ir->pbcType != PbcType::No)
4660 if (ir->nstlist == 0)
4663 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4664 "atoms might cause the simulation to crash.");
4666 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4669 "ERROR: The cut-off length is longer than half the shortest box vector or "
4670 "longer than the smallest box diagonal element. Increase the box size or "
4671 "decrease rlist.\n");
4672 warning_error(wi, warn_buf);