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"
99 using gmx::BasicVector;
101 /* Resource parameters
102 * Do not change any of these until you read the instruction
103 * in readinp.h. Some cpp's do not take spaces after the backslash
104 * (like the c-shell), which will give you a very weird compiler
108 struct gmx_inputrec_strings
110 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], freeze[STRLEN], frdim[STRLEN],
111 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
112 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
113 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
114 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::string> fep_lambda;
115 char lambda_weights[STRLEN];
116 std::vector<std::string> pullGroupNames;
117 std::vector<std::string> rotateGroupNames;
118 char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
121 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
122 static gmx_inputrec_strings* inputrecStrings = nullptr;
124 void init_inputrec_strings()
129 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
130 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
132 inputrecStrings = new gmx_inputrec_strings();
135 void done_inputrec_strings()
137 delete inputrecStrings;
138 inputrecStrings = nullptr;
144 egrptpALL, /* All particles have to be a member of a group. */
145 egrptpALL_GENREST, /* A rest group with name is generated for particles *
146 * that are not part of any group. */
147 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
148 * for the rest group. */
149 egrptpONE /* Merge all selected groups into one group, *
150 * make a rest group for the remaining particles. */
153 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
154 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
155 "h-angles", "all-angles", nullptr };
157 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
158 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
160 static void getSimTemps(int ntemps, t_simtemp* simtemp, gmx::ArrayRef<double> temperature_lambdas)
165 for (i = 0; i < ntemps; i++)
167 /* simple linear scaling -- allows more control */
168 if (simtemp->eSimTempScale == SimulatedTempering::Linear)
170 simtemp->temperatures[i] =
172 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
174 else if (simtemp->eSimTempScale
175 == SimulatedTempering::Geometric) /* should give roughly equal acceptance for constant heat capacity . . . */
177 simtemp->temperatures[i] = simtemp->simtemp_low
178 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
179 static_cast<real>((1.0 * i) / (ntemps - 1)));
181 else if (simtemp->eSimTempScale == SimulatedTempering::Exponential)
183 simtemp->temperatures[i] = simtemp->simtemp_low
184 + (simtemp->simtemp_high - simtemp->simtemp_low)
185 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
190 sprintf(errorstr, "eSimTempScale=%s not defined", enumValueToString(simtemp->eSimTempScale));
191 gmx_fatal(FARGS, "%s", errorstr);
197 static void _low_check(bool b, const char* s, warninp_t wi)
201 warning_error(wi, s);
205 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
209 if (*p > 0 && *p % nst != 0)
211 /* Round up to the next multiple of nst */
212 *p = ((*p) / nst + 1) * nst;
213 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
218 //! Convert legacy mdp entries to modern ones.
219 static void process_interaction_modifier(InteractionModifiers* eintmod)
221 if (*eintmod == InteractionModifiers::PotShiftVerletUnsupported)
223 *eintmod = InteractionModifiers::PotShift;
227 void check_ir(const char* mdparin,
228 const gmx::MDModulesNotifiers& mdModulesNotifiers,
232 /* Check internal consistency.
233 * NOTE: index groups are not set here yet, don't check things
234 * like temperature coupling group options here, but in triple_check
237 /* Strange macro: first one fills the err_buf, and then one can check
238 * the condition, which will print the message and increase the error
241 #define CHECK(b) _low_check(b, err_buf, wi)
242 char err_buf[256], warn_buf[STRLEN];
245 t_lambda* fep = ir->fepvals.get();
246 t_expanded* expand = ir->expandedvals.get();
248 set_warning_line(wi, mdparin, -1);
250 /* We cannot check MTS requirements with an invalid MTS setup
251 * and we will already have generated errors with an invalid MTS setup.
253 if (gmx::haveValidMtsSetup(*ir))
255 std::vector<std::string> errorMessages = gmx::checkMtsRequirements(*ir);
257 for (const auto& errorMessage : errorMessages)
259 warning_error(wi, errorMessage.c_str());
263 if (ir->coulombtype == CoulombInteractionType::RFNecUnsupported)
265 std::string message =
266 gmx::formatString("%s electrostatics is no longer supported",
267 enumValueToString(CoulombInteractionType::RFNecUnsupported));
268 warning_error(wi, message);
271 /* BASIC CUT-OFF STUFF */
272 if (ir->rcoulomb < 0)
274 warning_error(wi, "rcoulomb should be >= 0");
278 warning_error(wi, "rvdw should be >= 0");
280 if (ir->rlist < 0 && !(ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0))
282 warning_error(wi, "rlist should be >= 0");
285 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
286 "neighbour-list update scheme for efficient buffering for improved energy "
287 "conservation, please use the Verlet cut-off scheme instead.)");
288 CHECK(ir->nstlist < 0);
290 process_interaction_modifier(&ir->coulomb_modifier);
291 process_interaction_modifier(&ir->vdw_modifier);
293 if (ir->cutoff_scheme == CutoffScheme::Group)
296 "The group cutoff scheme has been removed since GROMACS 2020. "
297 "Please use the Verlet cutoff scheme.");
299 if (ir->cutoff_scheme == CutoffScheme::Verlet)
303 /* Normal Verlet type neighbor-list, currently only limited feature support */
304 if (inputrec2nboundeddim(ir) < 3)
306 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
309 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
310 if (ir->rcoulomb != ir->rvdw)
312 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
313 // for PME load balancing, we can support this exception.
314 bool bUsesPmeTwinRangeKernel =
315 (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == VanDerWaalsType::Cut
316 && ir->rcoulomb > ir->rvdw);
317 if (!bUsesPmeTwinRangeKernel)
320 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
321 "rcoulomb>rvdw with PME electrostatics)");
325 if (ir->vdwtype == VanDerWaalsType::Shift || ir->vdwtype == VanDerWaalsType::Switch)
327 if (ir->vdw_modifier == InteractionModifiers::None
328 || ir->vdw_modifier == InteractionModifiers::PotShift)
331 (ir->vdwtype == VanDerWaalsType::Shift ? InteractionModifiers::ForceSwitch
332 : InteractionModifiers::PotSwitch);
335 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
337 enumValueToString(ir->vdwtype),
338 enumValueToString(VanDerWaalsType::Cut),
339 enumValueToString(ir->vdw_modifier));
340 warning_note(wi, warn_buf);
342 ir->vdwtype = VanDerWaalsType::Cut;
347 "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
348 enumValueToString(ir->vdwtype),
349 enumValueToString(ir->vdw_modifier));
350 warning_error(wi, warn_buf);
354 if (!(ir->vdwtype == VanDerWaalsType::Cut || ir->vdwtype == VanDerWaalsType::Pme))
357 "With Verlet lists only cut-off and PME LJ interactions are supported");
359 if (!(ir->coulombtype == CoulombInteractionType::Cut || EEL_RF(ir->coulombtype)
360 || EEL_PME(ir->coulombtype) || ir->coulombtype == CoulombInteractionType::Ewald))
363 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
364 "electrostatics are supported");
366 if (!(ir->coulomb_modifier == InteractionModifiers::None
367 || ir->coulomb_modifier == InteractionModifiers::PotShift))
369 sprintf(warn_buf, "coulomb_modifier=%s is not supported", enumValueToString(ir->coulomb_modifier));
370 warning_error(wi, warn_buf);
373 if (EEL_USER(ir->coulombtype))
376 "Coulomb type %s is not supported with the verlet scheme",
377 enumValueToString(ir->coulombtype));
378 warning_error(wi, warn_buf);
381 if (ir->nstlist <= 0)
383 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
386 if (ir->nstlist < 10)
389 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
390 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
394 rc_max = std::max(ir->rvdw, ir->rcoulomb);
398 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
399 ir->verletbuf_tol = 0;
402 else if (ir->verletbuf_tol <= 0)
404 if (ir->verletbuf_tol == 0)
406 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
409 if (ir->rlist < rc_max)
412 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
415 if (ir->rlist == rc_max && ir->nstlist > 1)
419 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
420 "buffer. The cluster pair list does have a buffering effect, but choosing "
421 "a larger rlist might be necessary for good energy conservation.");
426 if (ir->rlist > rc_max)
429 "You have set rlist larger than the interaction cut-off, but you also "
430 "have verlet-buffer-tolerance > 0. Will set rlist using "
431 "verlet-buffer-tolerance.");
434 if (ir->nstlist == 1)
436 /* No buffer required */
441 if (EI_DYNAMICS(ir->eI))
443 if (inputrec2nboundeddim(ir) < 3)
446 "The box volume is required for calculating rlist from the "
447 "energy drift with verlet-buffer-tolerance > 0. You are "
448 "using at least one unbounded dimension, so no volume can be "
449 "computed. Either use a finite box, or set rlist yourself "
450 "together with verlet-buffer-tolerance = -1.");
452 /* Set rlist temporarily so we can continue processing */
457 /* Set the buffer to 5% of the cut-off */
458 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
464 /* GENERAL INTEGRATOR STUFF */
467 if (ir->etc != TemperatureCoupling::No)
469 if (EI_RANDOM(ir->eI))
472 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
473 "implicitly. See the documentation for more information on which "
474 "parameters affect temperature for %s.",
475 enumValueToString(ir->etc),
476 enumValueToString(ir->eI),
477 enumValueToString(ir->eI));
482 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
484 enumValueToString(ir->etc),
485 enumValueToString(ir->eI));
487 warning_note(wi, warn_buf);
489 ir->etc = TemperatureCoupling::No;
491 if (ir->eI == IntegrationAlgorithm::VVAK)
494 "Integrator method %s is implemented primarily for validation purposes; for "
495 "molecular dynamics, you should probably be using %s or %s",
496 enumValueToString(IntegrationAlgorithm::VVAK),
497 enumValueToString(IntegrationAlgorithm::MD),
498 enumValueToString(IntegrationAlgorithm::VV));
499 warning_note(wi, warn_buf);
501 if (!EI_DYNAMICS(ir->eI))
503 if (ir->epc != PressureCoupling::No)
506 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
507 enumValueToString(ir->epc),
508 enumValueToString(ir->eI));
509 warning_note(wi, warn_buf);
511 ir->epc = PressureCoupling::No;
513 if (EI_DYNAMICS(ir->eI))
515 if (ir->nstcalcenergy < 0)
517 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
518 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
520 /* nstcalcenergy larger than nstener does not make sense.
521 * We ideally want nstcalcenergy=nstener.
525 ir->nstcalcenergy = std::gcd(ir->nstenergy, ir->nstlist);
529 ir->nstcalcenergy = ir->nstenergy;
533 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
534 || (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl > 0
535 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
538 const char* nsten = "nstenergy";
539 const char* nstdh = "nstdhdl";
540 const char* min_name = nsten;
541 int min_nst = ir->nstenergy;
543 /* find the smallest of ( nstenergy, nstdhdl ) */
544 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl > 0
545 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
547 min_nst = ir->fepvals->nstdhdl;
550 /* If the user sets nstenergy small, we should respect that */
551 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy, min_name, min_nst);
552 warning_note(wi, warn_buf);
553 ir->nstcalcenergy = min_nst;
556 if (ir->epc != PressureCoupling::No)
558 if (ir->nstpcouple < 0)
560 ir->nstpcouple = ir_optimal_nstpcouple(ir);
562 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
565 "With multiple time stepping, nstpcouple should be a mutiple of "
570 if (ir->nstcalcenergy > 0)
572 if (ir->efep != FreeEnergyPerturbationType::No)
574 /* nstdhdl should be a multiple of nstcalcenergy */
575 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
579 /* nstexpanded should be a multiple of nstcalcenergy */
580 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded", &ir->expandedvals->nstexpanded, wi);
582 /* for storing exact averages nstenergy should be
583 * a multiple of nstcalcenergy
585 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
588 // Inquire all MDModules, if their parameters match with the energy
589 // calculation frequency
590 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
591 mdModulesNotifiers.preProcessingNotifier_.notify(&energyCalculationFrequencyErrors);
593 // Emit all errors from the energy calculation frequency checks
594 for (const std::string& energyFrequencyErrorMessage :
595 energyCalculationFrequencyErrors.errorMessages())
597 warning_error(wi, energyFrequencyErrorMessage);
601 if (ir->nsteps == 0 && !ir->bContinuation)
604 "For a correct single-point energy evaluation with nsteps = 0, use "
605 "continuation = yes to avoid constraining the input coordinates.");
609 if ((EI_SD(ir->eI) || ir->eI == IntegrationAlgorithm::BD) && ir->bContinuation && ir->ld_seed != -1)
612 "You are doing a continuation with SD or BD, make sure that ld_seed is "
613 "different from the previous run (using ld_seed=-1 will ensure this)");
619 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
620 CHECK(ir->pbcType != PbcType::Xyz);
621 sprintf(err_buf, "with TPI nstlist should be larger than zero");
622 CHECK(ir->nstlist <= 0);
623 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
624 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
628 if ((opts->nshake > 0) && (opts->bMorse))
630 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
631 warning(wi, warn_buf);
634 if ((EI_SD(ir->eI) || ir->eI == IntegrationAlgorithm::BD) && ir->bContinuation && ir->ld_seed != -1)
637 "You are doing a continuation with SD or BD, make sure that ld_seed is "
638 "different from the previous run (using ld_seed=-1 will ensure this)");
640 /* verify simulated tempering options */
644 bool bAllTempZero = TRUE;
645 for (i = 0; i < fep->n_lambda; i++)
648 "Entry %d for %s must be between 0 and 1, instead is %g",
650 enumValueToString(FreeEnergyPerturbationCouplingType::Temperature),
651 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i]);
652 CHECK((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i] < 0)
653 || (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i]
655 if (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i] > 0)
657 bAllTempZero = FALSE;
660 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
661 CHECK(bAllTempZero == TRUE);
663 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
664 CHECK(ir->eI != IntegrationAlgorithm::VV);
666 /* check compatability of the temperature coupling with simulated tempering */
668 if (ir->etc == TemperatureCoupling::NoseHoover)
671 "Nose-Hoover based temperature control such as [%s] my not be "
672 "entirelyconsistent with simulated tempering",
673 enumValueToString(ir->etc));
674 warning_note(wi, warn_buf);
677 /* check that the temperatures make sense */
680 "Higher simulated tempering temperature (%g) must be >= than the simulated "
681 "tempering lower temperature (%g)",
682 ir->simtempvals->simtemp_high,
683 ir->simtempvals->simtemp_low);
684 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
687 "Higher simulated tempering temperature (%g) must be >= zero",
688 ir->simtempvals->simtemp_high);
689 CHECK(ir->simtempvals->simtemp_high <= 0);
692 "Lower simulated tempering temperature (%g) must be >= zero",
693 ir->simtempvals->simtemp_low);
694 CHECK(ir->simtempvals->simtemp_low <= 0);
697 /* verify free energy options */
699 if (ir->efep != FreeEnergyPerturbationType::No)
701 fep = ir->fepvals.get();
702 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
703 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
706 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
708 static_cast<int>(fep->sc_r_power));
709 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
712 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
715 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
718 "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
720 CHECK(fep->delta_lambda > 0 && (ir->efep == FreeEnergyPerturbationType::Expanded));
722 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
723 CHECK(!(EI_VV(ir->eI)) && (ir->efep == FreeEnergyPerturbationType::Expanded));
725 sprintf(err_buf, "Free-energy not implemented for Ewald");
726 CHECK(ir->coulombtype == CoulombInteractionType::Ewald);
728 /* check validty of lambda inputs */
729 if (fep->n_lambda == 0)
731 /* Clear output in case of no states:*/
732 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
733 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
738 "initial thermodynamic state %d does not exist, only goes to %d",
741 CHECK((fep->init_fep_state >= fep->n_lambda));
745 "Lambda state must be set, either with init-lambda-state or with init-lambda");
746 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
749 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
750 "init-lambda-state or with init-lambda, but not both",
752 fep->init_fep_state);
753 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
756 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
760 for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
762 if (fep->separate_dvdl[i])
767 if (n_lambda_terms > 1)
770 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
771 "use init-lambda to set lambda state (except for slow growth). Use "
772 "init-lambda-state instead.");
773 warning(wi, warn_buf);
776 if (n_lambda_terms < 2 && fep->n_lambda > 0)
779 "init-lambda is deprecated for setting lambda state (except for slow "
780 "growth). Use init-lambda-state instead.");
784 for (j = 0; j < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); j++)
786 for (i = 0; i < fep->n_lambda; i++)
788 auto enumValue = static_cast<FreeEnergyPerturbationCouplingType>(j);
790 "Entry %d for %s must be between 0 and 1, instead is %g",
792 enumValueToString(enumValue),
793 fep->all_lambda[j][i]);
794 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
798 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
800 for (i = 0; i < fep->n_lambda; i++)
803 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
804 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
805 "crashes, and is not supported.",
807 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i],
808 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i]);
809 CHECK((fep->sc_alpha > 0)
810 && (((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i] > 0.0)
811 && (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i] < 1.0))
812 && ((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i] > 0.0)
813 && (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i]
818 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
820 real sigma, lambda, r_sc;
823 /* Maximum estimate for A and B charges equal with lambda power 1 */
825 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
826 1.0 / fep->sc_r_power);
828 "With PME there is a minor soft core effect present at the cut-off, "
829 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
830 "energy conservation, but usually other effects dominate. With a common sigma "
831 "value of %g nm the fraction of the particle-particle potential at the cut-off "
832 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
838 warning_note(wi, warn_buf);
841 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
842 be treated differently, but that's the next step */
844 for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
846 auto enumValue = static_cast<FreeEnergyPerturbationCouplingType>(i);
847 for (j = 0; j < fep->n_lambda; j++)
849 sprintf(err_buf, "%s[%d] must be between 0 and 1", enumValueToString(enumValue), j);
850 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
855 if ((ir->bSimTemp) || (ir->efep == FreeEnergyPerturbationType::Expanded))
857 fep = ir->fepvals.get();
859 /* checking equilibration of weights inputs for validity */
862 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
864 expand->equil_n_at_lam,
865 enumValueToString(LambdaWeightWillReachEquilibrium::NumAtLambda));
866 CHECK((expand->equil_n_at_lam > 0)
867 && (expand->elmceq != LambdaWeightWillReachEquilibrium::NumAtLambda));
870 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
872 expand->equil_samples,
873 enumValueToString(LambdaWeightWillReachEquilibrium::Samples));
874 CHECK((expand->equil_samples > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Samples));
877 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
879 enumValueToString(LambdaWeightWillReachEquilibrium::Steps));
880 CHECK((expand->equil_steps > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Steps));
883 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
884 expand->equil_samples,
885 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta));
886 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::WLDelta));
889 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
891 enumValueToString(LambdaWeightWillReachEquilibrium::Ratio));
892 CHECK((expand->equil_ratio > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Ratio));
895 "weight-equil-number-all-lambda (%d) must be a positive integer if "
896 "lmc-weights-equil=%s",
897 expand->equil_n_at_lam,
898 enumValueToString(LambdaWeightWillReachEquilibrium::NumAtLambda));
899 CHECK((expand->equil_n_at_lam <= 0)
900 && (expand->elmceq == LambdaWeightWillReachEquilibrium::NumAtLambda));
903 "weight-equil-number-samples (%d) must be a positive integer if "
904 "lmc-weights-equil=%s",
905 expand->equil_samples,
906 enumValueToString(LambdaWeightWillReachEquilibrium::Samples));
907 CHECK((expand->equil_samples <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Samples));
910 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
912 enumValueToString(LambdaWeightWillReachEquilibrium::Steps));
913 CHECK((expand->equil_steps <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Steps));
916 "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
917 expand->equil_wl_delta,
918 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta));
919 CHECK((expand->equil_wl_delta <= 0)
920 && (expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta));
923 "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
925 enumValueToString(LambdaWeightWillReachEquilibrium::Ratio));
926 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Ratio));
929 "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
930 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta),
931 enumValueToString(LambdaWeightCalculation::WL),
932 enumValueToString(LambdaWeightCalculation::WWL));
933 CHECK((expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta) && (!EWL(expand->elamstats)));
935 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
936 CHECK((expand->lmc_repeats <= 0));
937 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
938 CHECK((expand->minvarmin <= 0));
939 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
940 CHECK((expand->c_range < 0));
942 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
945 expand->lmc_forced_nstart);
946 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
947 && (expand->elmcmove != LambdaMoveCalculation::No));
948 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
949 CHECK((expand->lmc_forced_nstart < 0));
951 "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
952 fep->init_fep_state);
953 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
955 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
956 CHECK((expand->init_wl_delta < 0));
957 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
958 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
959 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
960 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
962 /* if there is no temperature control, we need to specify an MC temperature */
963 if (!integratorHasReferenceTemperature(ir)
964 && (expand->elmcmove != LambdaMoveCalculation::No) && (expand->mc_temp <= 0.0))
967 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
968 "to a positive number");
969 warning_error(wi, err_buf);
971 if (expand->nstTij > 0)
973 sprintf(err_buf, "nstlog must be non-zero");
974 CHECK(ir->nstlog == 0);
975 // Avoid modulus by zero in the case that already triggered an error exit.
979 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
982 CHECK((expand->nstTij % ir->nstlog) != 0);
988 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
989 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
992 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
994 if (ir->pbcType == PbcType::No)
996 if (ir->epc != PressureCoupling::No)
998 warning(wi, "Turning off pressure coupling for vacuum system");
999 ir->epc = PressureCoupling::No;
1005 "Can not have pressure coupling with pbc=%s",
1006 c_pbcTypeNames[ir->pbcType].c_str());
1007 CHECK(ir->epc != PressureCoupling::No);
1009 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
1010 CHECK(EEL_FULL(ir->coulombtype));
1013 "Can not have dispersion correction with pbc=%s",
1014 c_pbcTypeNames[ir->pbcType].c_str());
1015 CHECK(ir->eDispCorr != DispersionCorrectionType::No);
1018 if (ir->rlist == 0.0)
1021 "can only have neighborlist cut-off zero (=infinite)\n"
1022 "with coulombtype = %s or coulombtype = %s\n"
1023 "without periodic boundary conditions (pbc = %s) and\n"
1024 "rcoulomb and rvdw set to zero",
1025 enumValueToString(CoulombInteractionType::Cut),
1026 enumValueToString(CoulombInteractionType::User),
1027 c_pbcTypeNames[PbcType::No].c_str());
1028 CHECK(((ir->coulombtype != CoulombInteractionType::Cut)
1029 && (ir->coulombtype != CoulombInteractionType::User))
1030 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1032 if (ir->nstlist > 0)
1035 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1036 "nstype=simple and only one MPI rank");
1041 if (ir->nstcomm == 0)
1043 // TODO Change this behaviour. There should be exactly one way
1044 // to turn off an algorithm.
1045 ir->comm_mode = ComRemovalAlgorithm::No;
1047 if (ir->comm_mode != ComRemovalAlgorithm::No)
1049 if (ir->nstcomm < 0)
1051 // TODO Such input was once valid. Now that we've been
1052 // helpful for a few years, we should reject such input,
1053 // lest we have to support every historical decision
1056 "If you want to remove the rotation around the center of mass, you should set "
1057 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1058 "its absolute value");
1059 ir->nstcomm = abs(ir->nstcomm);
1062 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy
1063 && ir->comm_mode != ComRemovalAlgorithm::LinearAccelerationCorrection)
1066 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, consider "
1067 "setting nstcomm equal to nstcalcenergy for less overhead");
1070 if (ir->comm_mode == ComRemovalAlgorithm::Angular)
1073 "Can not remove the rotation around the center of mass with periodic "
1075 CHECK(ir->bPeriodicMols);
1076 if (ir->pbcType != PbcType::No)
1079 "Removing the rotation around the center of mass in a periodic system, "
1080 "this can lead to artifacts. Only use this on a single (cluster of) "
1081 "molecules. This cluster should not cross periodic boundaries.");
1086 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No
1087 && ir->comm_mode != ComRemovalAlgorithm::Angular)
1090 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1091 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1093 enumValueToString(IntegrationAlgorithm::SD1));
1094 warning_note(wi, warn_buf);
1097 /* TEMPERATURE COUPLING */
1098 if (ir->etc == TemperatureCoupling::Yes)
1100 ir->etc = TemperatureCoupling::Berendsen;
1102 "Old option for temperature coupling given: "
1103 "changing \"yes\" to \"Berendsen\"\n");
1106 if ((ir->etc == TemperatureCoupling::NoseHoover) || (ir->epc == PressureCoupling::Mttk))
1108 if (ir->opts.nhchainlength < 1)
1111 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1113 ir->opts.nhchainlength);
1114 ir->opts.nhchainlength = 1;
1115 warning(wi, warn_buf);
1118 if (ir->etc == TemperatureCoupling::NoseHoover && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1122 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1123 ir->opts.nhchainlength = 1;
1128 ir->opts.nhchainlength = 0;
1131 if (ir->eI == IntegrationAlgorithm::VVAK)
1134 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1136 enumValueToString(IntegrationAlgorithm::VVAK));
1137 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1140 if (ETC_ANDERSEN(ir->etc))
1143 "%s temperature control not supported for integrator %s.",
1144 enumValueToString(ir->etc),
1145 enumValueToString(ir->eI));
1146 CHECK(!(EI_VV(ir->eI)));
1148 if (ir->nstcomm > 0 && (ir->etc == TemperatureCoupling::Andersen))
1151 "Center of mass removal not necessary for %s. All velocities of coupled "
1152 "groups are rerandomized periodically, so flying ice cube errors will not "
1154 enumValueToString(ir->etc));
1155 warning_note(wi, warn_buf);
1159 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1160 "randomized every time step",
1162 enumValueToString(ir->etc));
1163 CHECK(ir->nstcomm > 1 && (ir->etc == TemperatureCoupling::Andersen));
1166 if (ir->etc == TemperatureCoupling::Berendsen)
1169 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1170 "might want to consider using the %s thermostat.",
1171 enumValueToString(ir->etc),
1172 enumValueToString(TemperatureCoupling::VRescale));
1173 warning_note(wi, warn_buf);
1176 if ((ir->etc == TemperatureCoupling::NoseHoover || ETC_ANDERSEN(ir->etc))
1177 && ir->epc == PressureCoupling::Berendsen)
1180 "Using Berendsen pressure coupling invalidates the "
1181 "true ensemble for the thermostat");
1182 warning(wi, warn_buf);
1185 /* PRESSURE COUPLING */
1186 if (ir->epc == PressureCoupling::Isotropic)
1188 ir->epc = PressureCoupling::Berendsen;
1190 "Old option for pressure coupling given: "
1191 "changing \"Isotropic\" to \"Berendsen\"\n");
1194 if (ir->epc != PressureCoupling::No)
1196 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1198 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1199 CHECK(ir->tau_p <= 0);
1201 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1204 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1205 "times larger than nstpcouple*dt (%g)",
1206 enumValueToString(ir->epc),
1208 pcouple_min_integration_steps(ir->epc),
1210 warning(wi, warn_buf);
1214 "compressibility must be > 0 when using pressure"
1216 enumValueToString(ir->epc));
1217 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1218 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1219 && ir->compress[ZZ][YY] <= 0));
1221 if (PressureCoupling::ParrinelloRahman == ir->epc && opts->bGenVel)
1224 "You are generating velocities so I am assuming you "
1225 "are equilibrating a system. You are using "
1226 "%s pressure coupling, but this can be "
1227 "unstable for equilibration. If your system crashes, try "
1228 "equilibrating first with Berendsen pressure coupling. If "
1229 "you are not equilibrating the system, you can probably "
1230 "ignore this warning.",
1231 enumValueToString(ir->epc));
1232 warning(wi, warn_buf);
1238 if (ir->epc == PressureCoupling::Mttk)
1240 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1244 /* ELECTROSTATICS */
1245 /* More checks are in triple check (grompp.c) */
1247 if (ir->coulombtype == CoulombInteractionType::Switch)
1250 "coulombtype = %s is only for testing purposes and can lead to serious "
1251 "artifacts, advice: use coulombtype = %s",
1252 enumValueToString(ir->coulombtype),
1253 enumValueToString(CoulombInteractionType::RFZero));
1254 warning(wi, warn_buf);
1257 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1260 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1261 "format and exchanging epsilon-r and epsilon-rf",
1263 warning(wi, warn_buf);
1264 ir->epsilon_rf = ir->epsilon_r;
1265 ir->epsilon_r = 1.0;
1268 if (ir->epsilon_r == 0)
1271 "It is pointless to use long-range electrostatics with infinite relative "
1273 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1275 CHECK(EEL_FULL(ir->coulombtype));
1278 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1280 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1281 CHECK(ir->epsilon_r < 0);
1284 if (EEL_RF(ir->coulombtype))
1286 /* reaction field (at the cut-off) */
1288 if (ir->coulombtype == CoulombInteractionType::RFZero && ir->epsilon_rf != 0)
1291 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1292 enumValueToString(ir->coulombtype));
1293 warning(wi, warn_buf);
1294 ir->epsilon_rf = 0.0;
1297 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1298 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1299 if (ir->epsilon_rf == ir->epsilon_r)
1302 "Using epsilon-rf = epsilon-r with %s does not make sense",
1303 enumValueToString(ir->coulombtype));
1304 warning(wi, warn_buf);
1307 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1308 * means the interaction is zero outside rcoulomb, but it helps to
1309 * provide accurate energy conservation.
1311 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1313 if (ir_coulomb_switched(ir))
1316 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1317 "potential modifier options!",
1318 enumValueToString(ir->coulombtype));
1319 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1323 if (ir->coulombtype == CoulombInteractionType::Switch || ir->coulombtype == CoulombInteractionType::Shift)
1326 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1327 "secondary coulomb-modifier.");
1328 CHECK(ir->coulomb_modifier != InteractionModifiers::None);
1330 if (ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift)
1333 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1334 "secondary vdw-modifier.");
1335 CHECK(ir->vdw_modifier != InteractionModifiers::None);
1338 if (ir->coulombtype == CoulombInteractionType::Switch || ir->coulombtype == CoulombInteractionType::Shift
1339 || ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift)
1342 "The switch/shift interaction settings are just for compatibility; you will get "
1344 "performance from applying potential modifiers to your interactions!\n");
1345 warning_note(wi, warn_buf);
1348 if (ir->coulombtype == CoulombInteractionType::PmeSwitch
1349 || ir->coulomb_modifier == InteractionModifiers::PotSwitch)
1351 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1353 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1355 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1356 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1357 "will be good regardless, since ewald_rtol = %g.",
1359 ir->rcoulomb_switch,
1362 warning(wi, warn_buf);
1366 if (ir->vdwtype == VanDerWaalsType::Switch || ir->vdw_modifier == InteractionModifiers::PotSwitch)
1368 if (ir->rvdw_switch == 0)
1371 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1372 "potential. This suggests it was not set in the mdp, which can lead to large "
1373 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1374 "switching range.");
1375 warning(wi, warn_buf);
1379 if (EEL_FULL(ir->coulombtype))
1381 if (ir->coulombtype == CoulombInteractionType::PmeSwitch
1382 || ir->coulombtype == CoulombInteractionType::PmeUser
1383 || ir->coulombtype == CoulombInteractionType::PmeUserSwitch)
1386 "With coulombtype = %s, rcoulomb must be <= rlist",
1387 enumValueToString(ir->coulombtype));
1388 CHECK(ir->rcoulomb > ir->rlist);
1392 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1394 // TODO: Move these checks into the ewald module with the options class
1396 int orderMax = (ir->coulombtype == CoulombInteractionType::P3mAD ? 8 : 12);
1398 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1401 "With coulombtype = %s, you should have %d <= pme-order <= %d",
1402 enumValueToString(ir->coulombtype),
1405 warning_error(wi, warn_buf);
1409 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1411 if (ir->ewald_geometry == EwaldGeometry::ThreeD)
1414 "With pbc=%s you should use ewald-geometry=%s",
1415 c_pbcTypeNames[ir->pbcType].c_str(),
1416 enumValueToString(EwaldGeometry::ThreeDC));
1417 warning(wi, warn_buf);
1419 /* This check avoids extra pbc coding for exclusion corrections */
1420 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1421 CHECK(ir->wall_ewald_zfac < 2);
1423 if ((ir->ewald_geometry == EwaldGeometry::ThreeDC) && (ir->pbcType != PbcType::XY)
1424 && EEL_FULL(ir->coulombtype))
1427 "With %s and ewald_geometry = %s you should use pbc = %s",
1428 enumValueToString(ir->coulombtype),
1429 enumValueToString(EwaldGeometry::ThreeDC),
1430 c_pbcTypeNames[PbcType::XY].c_str());
1431 warning(wi, warn_buf);
1433 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1435 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1436 CHECK(ir->bPeriodicMols);
1437 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1438 warning_note(wi, warn_buf);
1440 "With epsilon_surface > 0 you can only use domain decomposition "
1441 "when there are only small molecules with all bonds constrained (mdrun will check "
1443 warning_note(wi, warn_buf);
1446 if (ir_vdw_switched(ir))
1448 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1449 CHECK(ir->rvdw_switch >= ir->rvdw);
1451 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1454 "You are applying a switch function to vdw forces or potentials from %g to %g "
1455 "nm, which is more than half the interaction range, whereas switch functions "
1456 "are intended to act only close to the cut-off.",
1459 warning_note(wi, warn_buf);
1463 if (ir->vdwtype == VanDerWaalsType::Pme)
1465 if (!(ir->vdw_modifier == InteractionModifiers::None
1466 || ir->vdw_modifier == InteractionModifiers::PotShift))
1469 "With vdwtype = %s, the only supported modifiers are %s and %s",
1470 enumValueToString(ir->vdwtype),
1471 enumValueToString(InteractionModifiers::PotShift),
1472 enumValueToString(InteractionModifiers::None));
1473 warning_error(wi, err_buf);
1477 if (ir->vdwtype == VanDerWaalsType::User && ir->eDispCorr != DispersionCorrectionType::No)
1480 "You have selected user tables with dispersion correction, the dispersion "
1481 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1482 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1483 "really want dispersion correction to -C6/r^6.");
1486 if (ir->eI == IntegrationAlgorithm::LBFGS
1487 && (ir->coulombtype == CoulombInteractionType::Cut || ir->vdwtype == VanDerWaalsType::Cut)
1490 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1493 if (ir->eI == IntegrationAlgorithm::LBFGS && ir->nbfgscorr <= 0)
1495 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1498 /* IMPLICIT SOLVENT */
1499 if (ir->coulombtype == CoulombInteractionType::GBNotused)
1501 sprintf(warn_buf, "Invalid option %s for coulombtype", enumValueToString(ir->coulombtype));
1502 warning_error(wi, warn_buf);
1507 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1512 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1515 // cosine acceleration is only supported in leap-frog
1516 if (ir->cos_accel != 0.0 && ir->eI != IntegrationAlgorithm::MD)
1518 warning_error(wi, "cos-acceleration is only supported by integrator = md");
1522 /* interpret a number of doubles from a string and put them in an array,
1523 after allocating space for them.
1524 str = the input string
1525 n = the (pre-allocated) number of doubles read
1526 r = the output array of doubles. */
1527 static std::vector<real> parse_n_real(const std::string& str, int* n, warninp_t wi)
1529 auto values = gmx::splitString(str);
1532 std::vector<real> r;
1533 for (int i = 0; i < *n; i++)
1537 r.emplace_back(gmx::fromString<real>(values[i]));
1539 catch (gmx::GromacsException&)
1542 "Invalid value " + values[i]
1543 + " in string in mdp file. Expected a real number.");
1550 static void do_fep_params(t_inputrec* ir, gmx::ArrayRef<std::string> fep_lambda, char weights[STRLEN], warninp_t wi)
1553 int i, j, max_n_lambda, nweights;
1554 t_lambda* fep = ir->fepvals.get();
1555 t_expanded* expand = ir->expandedvals.get();
1556 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::vector<real>> count_fep_lambdas;
1557 bool bOneLambda = TRUE;
1558 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, int> nfep;
1560 /* FEP input processing */
1561 /* first, identify the number of lambda values for each type.
1562 All that are nonzero must have the same number */
1564 for (auto i : keysOf(nfep))
1566 count_fep_lambdas[i] = parse_n_real(fep_lambda[static_cast<int>(i)], &(nfep[i]), wi);
1569 /* now, determine the number of components. All must be either zero, or equal. */
1572 for (auto i : keysOf(nfep))
1574 if (nfep[i] > max_n_lambda)
1576 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1577 must have the same number if its not zero.*/
1582 for (auto i : keysOf(nfep))
1586 ir->fepvals->separate_dvdl[i] = FALSE;
1588 else if (nfep[i] == max_n_lambda)
1590 if (i != FreeEnergyPerturbationCouplingType::Temperature) /* we treat this differently -- not really a reason to compute
1591 the derivative with respect to the temperature currently */
1593 ir->fepvals->separate_dvdl[i] = TRUE;
1599 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1602 enumValueToString(i),
1606 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1607 ir->fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Temperature] = FALSE;
1609 /* the number of lambdas is the number we've read in, which is either zero
1610 or the same for all */
1611 fep->n_lambda = max_n_lambda;
1613 /* if init_lambda is defined, we need to set lambda */
1614 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1616 ir->fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep] = TRUE;
1618 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1619 for (auto i : keysOf(nfep))
1621 fep->all_lambda[i].resize(fep->n_lambda);
1622 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1625 for (j = 0; j < fep->n_lambda; j++)
1627 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1632 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1633 internal bookkeeping -- for now, init_lambda */
1635 if ((nfep[FreeEnergyPerturbationCouplingType::Fep] == 0) && (fep->init_lambda >= 0))
1637 for (i = 0; i < fep->n_lambda; i++)
1639 fep->all_lambda[FreeEnergyPerturbationCouplingType::Fep][i] = fep->init_lambda;
1643 /* check to see if only a single component lambda is defined, and soft core is defined.
1644 In this case, turn on coulomb soft core */
1646 if (max_n_lambda == 0)
1652 for (auto i : keysOf(nfep))
1654 if ((nfep[i] != 0) && (i != FreeEnergyPerturbationCouplingType::Fep))
1660 if ((bOneLambda) && (fep->sc_alpha > 0))
1662 fep->bScCoul = TRUE;
1665 /* Fill in the others with the efptFEP if they are not explicitly
1666 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1667 they are all zero. */
1669 for (auto i : keysOf(nfep))
1671 if ((nfep[i] == 0) && (i != FreeEnergyPerturbationCouplingType::Fep))
1673 for (j = 0; j < fep->n_lambda; j++)
1675 fep->all_lambda[i][j] = fep->all_lambda[FreeEnergyPerturbationCouplingType::Fep][j];
1681 /* now read in the weights */
1682 expand->init_lambda_weights = parse_n_real(weights, &nweights, wi);
1685 expand->init_lambda_weights.resize(fep->n_lambda); /* initialize to zero */
1687 else if (nweights != fep->n_lambda)
1690 "Number of weights (%d) is not equal to number of lambda values (%d)",
1694 if ((expand->nstexpanded < 0) && (ir->efep != FreeEnergyPerturbationType::No))
1696 expand->nstexpanded = fep->nstdhdl;
1697 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1702 static void do_simtemp_params(t_inputrec* ir)
1704 ir->simtempvals->temperatures.resize(ir->fepvals->n_lambda);
1705 getSimTemps(ir->fepvals->n_lambda,
1706 ir->simtempvals.get(),
1707 ir->fepvals->all_lambda[FreeEnergyPerturbationCouplingType::Temperature]);
1710 template<typename T>
1711 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1714 for (const auto& input : inputs)
1718 outputs[i] = gmx::fromStdString<T>(input);
1720 catch (gmx::GromacsException&)
1722 auto message = gmx::formatString(
1723 "Invalid value for mdp option %s. %s should only consist of integers separated "
1727 warning_error(wi, message);
1733 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1736 for (const auto& input : inputs)
1740 outputs[i] = gmx::fromString<real>(input);
1742 catch (gmx::GromacsException&)
1744 auto message = gmx::formatString(
1745 "Invalid value for mdp option %s. %s should only consist of real numbers "
1746 "separated by spaces.",
1749 warning_error(wi, message);
1755 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1757 opts->wall_atomtype[0] = nullptr;
1758 opts->wall_atomtype[1] = nullptr;
1760 ir->wall_atomtype[0] = -1;
1761 ir->wall_atomtype[1] = -1;
1762 ir->wall_density[0] = 0;
1763 ir->wall_density[1] = 0;
1767 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1768 if (wallAtomTypes.size() != size_t(ir->nwall))
1771 "Expected %d elements for wall_atomtype, found %zu",
1773 wallAtomTypes.size());
1775 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1776 for (int i = 0; i < ir->nwall; i++)
1778 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1781 if (ir->wall_type == WallType::NineThree || ir->wall_type == WallType::TenFour)
1783 auto wallDensity = gmx::splitString(wall_density);
1784 if (wallDensity.size() != size_t(ir->nwall))
1787 "Expected %d elements for wall-density, found %zu",
1789 wallDensity.size());
1791 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1792 for (int i = 0; i < ir->nwall; i++)
1794 if (ir->wall_density[i] <= 0)
1796 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1803 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1807 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1808 for (int i = 0; i < nwall; i++)
1810 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1811 grps->emplace_back(groups->groupNames.size() - 1);
1816 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1818 /* read expanded ensemble parameters */
1819 printStringNewline(inp, "expanded ensemble variables");
1820 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1821 expand->elamstats = getEnum<LambdaWeightCalculation>(inp, "lmc-stats", wi);
1822 expand->elmcmove = getEnum<LambdaMoveCalculation>(inp, "lmc-move", wi);
1823 expand->elmceq = getEnum<LambdaWeightWillReachEquilibrium>(inp, "lmc-weights-equil", wi);
1824 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1825 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1826 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1827 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1828 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1829 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1830 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1831 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1832 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1833 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1834 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1835 expand->bSymmetrizedTMatrix =
1836 (getEnum<Boolean>(inp, "symmetrized-transition-matrix", wi) != Boolean::No);
1837 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1838 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1839 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1840 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1841 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1842 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1843 expand->bWLoneovert = (getEnum<Boolean>(inp, "wl-oneovert", wi) != Boolean::No);
1846 /*! \brief Return whether an end state with the given coupling-lambda
1847 * value describes fully-interacting VDW.
1849 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1850 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1852 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1854 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1860 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1863 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1865 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1867 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1870 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1871 std::string message = gmx::formatExceptionMessageToString(*ex);
1872 warning_error(wi_, message.c_str());
1877 std::string getOptionName(const gmx::KeyValueTreePath& context)
1879 if (mapping_ != nullptr)
1881 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1882 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1885 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1890 const gmx::IKeyValueTreeBackMapping* mapping_;
1895 void get_ir(const char* mdparin,
1896 const char* mdparout,
1897 gmx::MDModules* mdModules,
1900 WriteMdpHeader writeMdpHeader,
1904 double dumdub[2][6];
1906 char warn_buf[STRLEN];
1907 t_lambda* fep = ir->fepvals.get();
1908 t_expanded* expand = ir->expandedvals.get();
1910 const char* no_names[] = { "no", nullptr };
1912 init_inputrec_strings();
1913 gmx::TextInputFile stream(mdparin);
1914 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1916 snew(dumstr[0], STRLEN);
1917 snew(dumstr[1], STRLEN);
1919 /* ignore the following deprecated commands */
1920 replace_inp_entry(inp, "title", nullptr);
1921 replace_inp_entry(inp, "cpp", nullptr);
1922 replace_inp_entry(inp, "domain-decomposition", nullptr);
1923 replace_inp_entry(inp, "andersen-seed", nullptr);
1924 replace_inp_entry(inp, "dihre", nullptr);
1925 replace_inp_entry(inp, "dihre-fc", nullptr);
1926 replace_inp_entry(inp, "dihre-tau", nullptr);
1927 replace_inp_entry(inp, "nstdihreout", nullptr);
1928 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1929 replace_inp_entry(inp, "optimize-fft", nullptr);
1930 replace_inp_entry(inp, "adress_type", nullptr);
1931 replace_inp_entry(inp, "adress_const_wf", nullptr);
1932 replace_inp_entry(inp, "adress_ex_width", nullptr);
1933 replace_inp_entry(inp, "adress_hy_width", nullptr);
1934 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1935 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1936 replace_inp_entry(inp, "adress_site", nullptr);
1937 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1938 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1939 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1940 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1941 replace_inp_entry(inp, "rlistlong", nullptr);
1942 replace_inp_entry(inp, "nstcalclr", nullptr);
1943 replace_inp_entry(inp, "pull-print-com2", nullptr);
1944 replace_inp_entry(inp, "gb-algorithm", nullptr);
1945 replace_inp_entry(inp, "nstgbradii", nullptr);
1946 replace_inp_entry(inp, "rgbradii", nullptr);
1947 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1948 replace_inp_entry(inp, "gb-saltconc", nullptr);
1949 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1950 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1951 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1952 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1953 replace_inp_entry(inp, "sa-algorithm", nullptr);
1954 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1955 replace_inp_entry(inp, "ns-type", nullptr);
1957 /* replace the following commands with the clearer new versions*/
1958 replace_inp_entry(inp, "unconstrained-start", "continuation");
1959 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1960 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1961 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1962 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1963 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1964 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1966 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1967 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1968 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1969 setStringEntry(&inp, "include", opts->include, nullptr);
1970 printStringNoNewline(
1971 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1972 setStringEntry(&inp, "define", opts->define, nullptr);
1974 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1975 ir->eI = getEnum<IntegrationAlgorithm>(&inp, "integrator", wi);
1976 printStringNoNewline(&inp, "Start time and timestep in ps");
1977 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1978 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1979 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1980 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1981 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1982 printStringNoNewline(
1983 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1984 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1985 printStringNoNewline(&inp, "Multiple time-stepping");
1986 ir->useMts = (getEnum<Boolean>(&inp, "mts", wi) != Boolean::No);
1989 gmx::GromppMtsOpts& mtsOpts = opts->mtsOpts;
1990 mtsOpts.numLevels = get_eint(&inp, "mts-levels", 2, wi);
1991 mtsOpts.level2Forces = setStringEntry(&inp, "mts-level2-forces", "longrange-nonbonded");
1992 mtsOpts.level2Factor = get_eint(&inp, "mts-level2-factor", 2, wi);
1994 // We clear after reading without dynamics to not force the user to remove MTS mdp options
1995 if (!EI_DYNAMICS(ir->eI))
2000 printStringNoNewline(&inp, "mode for center of mass motion removal");
2001 ir->comm_mode = getEnum<ComRemovalAlgorithm>(&inp, "comm-mode", wi);
2002 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
2003 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
2004 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
2005 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
2007 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2008 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2009 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2010 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2013 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2014 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2015 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2016 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2017 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2018 ir->niter = get_eint(&inp, "niter", 20, wi);
2019 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2020 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2021 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2022 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2023 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2025 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2026 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2028 /* Output options */
2029 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2030 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2031 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2032 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2033 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2034 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2035 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2036 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2037 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2038 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2039 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2040 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2041 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2042 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2043 printStringNoNewline(&inp, "default, all atoms will be written.");
2044 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2045 printStringNoNewline(&inp, "Selection of energy groups");
2046 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2048 /* Neighbor searching */
2049 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2050 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2051 ir->cutoff_scheme = getEnum<CutoffScheme>(&inp, "cutoff-scheme", wi);
2052 printStringNoNewline(&inp, "nblist update frequency");
2053 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2054 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2055 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2056 std::vector<const char*> pbcTypesNamesChar;
2057 for (const auto& pbcTypeName : c_pbcTypeNames)
2059 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2061 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2062 ir->bPeriodicMols = getEnum<Boolean>(&inp, "periodic-molecules", wi) != Boolean::No;
2063 printStringNoNewline(&inp,
2064 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2065 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2066 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2067 printStringNoNewline(&inp, "nblist cut-off");
2068 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2069 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2071 /* Electrostatics */
2072 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2073 printStringNoNewline(&inp, "Method for doing electrostatics");
2074 ir->coulombtype = getEnum<CoulombInteractionType>(&inp, "coulombtype", wi);
2075 ir->coulomb_modifier = getEnum<InteractionModifiers>(&inp, "coulomb-modifier", wi);
2076 printStringNoNewline(&inp, "cut-off lengths");
2077 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2078 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2079 printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
2080 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2081 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2082 printStringNoNewline(&inp, "Method for doing Van der Waals");
2083 ir->vdwtype = getEnum<VanDerWaalsType>(&inp, "vdw-type", wi);
2084 ir->vdw_modifier = getEnum<InteractionModifiers>(&inp, "vdw-modifier", wi);
2085 printStringNoNewline(&inp, "cut-off lengths");
2086 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2087 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2088 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2089 ir->eDispCorr = getEnum<DispersionCorrectionType>(&inp, "DispCorr", wi);
2090 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2091 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2092 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2093 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2094 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2095 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2096 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2097 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2098 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2099 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2100 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2101 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2102 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2103 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2104 ir->ljpme_combination_rule = getEnum<LongRangeVdW>(&inp, "lj-pme-comb-rule", wi);
2105 ir->ewald_geometry = getEnum<EwaldGeometry>(&inp, "ewald-geometry", wi);
2106 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2108 /* Implicit solvation is no longer supported, but we need grompp
2109 to be able to refuse old .mdp files that would have built a tpr
2110 to run it. Thus, only "no" is accepted. */
2111 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2113 /* Coupling stuff */
2114 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2115 printStringNoNewline(&inp, "Temperature coupling");
2116 ir->etc = getEnum<TemperatureCoupling>(&inp, "tcoupl", wi);
2117 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2118 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2119 ir->bPrintNHChains = (getEnum<Boolean>(&inp, "print-nose-hoover-chain-variables", wi) != Boolean::No);
2120 printStringNoNewline(&inp, "Groups to couple separately");
2121 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2122 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2123 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2124 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2125 printStringNoNewline(&inp, "pressure coupling");
2126 ir->epc = getEnum<PressureCoupling>(&inp, "pcoupl", wi);
2127 ir->epct = getEnum<PressureCouplingType>(&inp, "pcoupltype", wi);
2128 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2129 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2130 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2131 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2132 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2133 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2134 ir->refcoord_scaling = getEnum<RefCoordScaling>(&inp, "refcoord-scaling", wi);
2137 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2138 ir->bQMMM = (getEnum<Boolean>(&inp, "QMMM", wi) != Boolean::No);
2139 printStringNoNewline(&inp, "Groups treated with MiMiC");
2140 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2142 /* Simulated annealing */
2143 printStringNewline(&inp, "SIMULATED ANNEALING");
2144 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2145 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2146 printStringNoNewline(&inp,
2147 "Number of time points to use for specifying annealing in each group");
2148 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2149 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2150 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2151 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2152 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2155 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2156 opts->bGenVel = (getEnum<Boolean>(&inp, "gen-vel", wi) != Boolean::No);
2157 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2158 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2161 printStringNewline(&inp, "OPTIONS FOR BONDS");
2162 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2163 printStringNoNewline(&inp, "Type of constraint algorithm");
2164 ir->eConstrAlg = getEnum<ConstraintAlgorithm>(&inp, "constraint-algorithm", wi);
2165 printStringNoNewline(&inp, "Do not constrain the start configuration");
2166 ir->bContinuation = (getEnum<Boolean>(&inp, "continuation", wi) != Boolean::No);
2167 printStringNoNewline(&inp,
2168 "Use successive overrelaxation to reduce the number of shake iterations");
2169 ir->bShakeSOR = (getEnum<Boolean>(&inp, "Shake-SOR", wi) != Boolean::No);
2170 printStringNoNewline(&inp, "Relative tolerance of shake");
2171 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2172 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2173 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2174 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2175 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2176 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2177 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2178 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2179 printStringNoNewline(&inp, "rotates over more degrees than");
2180 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2181 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2182 opts->bMorse = (getEnum<Boolean>(&inp, "morse", wi) != Boolean::No);
2184 /* Energy group exclusions */
2185 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2186 printStringNoNewline(
2187 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2188 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2191 printStringNewline(&inp, "WALLS");
2192 printStringNoNewline(
2193 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2194 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2195 ir->wall_type = getEnum<WallType>(&inp, "wall-type", wi);
2196 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2197 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2198 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2199 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2202 printStringNewline(&inp, "COM PULLING");
2203 ir->bPull = (getEnum<Boolean>(&inp, "pull", wi) != Boolean::No);
2206 ir->pull = std::make_unique<pull_params_t>();
2207 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
2211 for (int c = 0; c < ir->pull->ncoord; c++)
2213 if (ir->pull->coord[c].eType == PullingAlgorithm::Constraint)
2216 "Constraint COM pulling is not supported in combination with "
2217 "multiple time stepping");
2225 NOTE: needs COM pulling or free energy input */
2226 printStringNewline(&inp, "AWH biasing");
2227 ir->bDoAwh = (getEnum<Boolean>(&inp, "awh", wi) != Boolean::No);
2230 ir->awhParams = std::make_unique<gmx::AwhParams>(&inp, wi);
2233 /* Enforced rotation */
2234 printStringNewline(&inp, "ENFORCED ROTATION");
2235 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2236 ir->bRot = (getEnum<Boolean>(&inp, "rotation", wi) != Boolean::No);
2240 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2243 /* Interactive MD */
2245 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2246 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2247 if (inputrecStrings->imd_grp[0] != '\0')
2254 printStringNewline(&inp, "NMR refinement stuff");
2255 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2256 ir->eDisre = getEnum<DistanceRestraintRefinement>(&inp, "disre", wi);
2257 printStringNoNewline(
2258 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2259 ir->eDisreWeighting = getEnum<DistanceRestraintWeighting>(&inp, "disre-weighting", wi);
2260 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2261 ir->bDisreMixed = (getEnum<Boolean>(&inp, "disre-mixed", wi) != Boolean::No);
2262 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2263 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2264 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2265 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2266 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2267 opts->bOrire = (getEnum<Boolean>(&inp, "orire", wi) != Boolean::No);
2268 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2269 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2270 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2271 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2272 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2273 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2275 /* free energy variables */
2276 printStringNewline(&inp, "Free energy variables");
2277 ir->efep = getEnum<FreeEnergyPerturbationType>(&inp, "free-energy", wi);
2278 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2279 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2280 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2281 opts->bCoupleIntra = (getEnum<Boolean>(&inp, "couple-intramol", wi) != Boolean::No);
2283 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2285 it was not entered */
2286 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2287 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2288 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2289 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Fep] =
2290 setStringEntry(&inp, "fep-lambdas", "");
2291 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Mass] =
2292 setStringEntry(&inp, "mass-lambdas", "");
2293 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Coul] =
2294 setStringEntry(&inp, "coul-lambdas", "");
2295 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Vdw] =
2296 setStringEntry(&inp, "vdw-lambdas", "");
2297 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Bonded] =
2298 setStringEntry(&inp, "bonded-lambdas", "");
2299 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Restraint] =
2300 setStringEntry(&inp, "restraint-lambdas", "");
2301 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Temperature] =
2302 setStringEntry(&inp, "temperature-lambdas", "");
2303 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2304 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2305 fep->edHdLPrintEnergy = getEnum<FreeEnergyPrintEnergy>(&inp, "dhdl-print-energy", wi);
2306 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2307 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2308 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2309 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2310 fep->bScCoul = (getEnum<Boolean>(&inp, "sc-coul", wi) != Boolean::No);
2311 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2312 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2313 fep->separate_dhdl_file = getEnum<SeparateDhdlFile>(&inp, "separate-dhdl-file", wi);
2314 fep->dhdl_derivatives = getEnum<DhDlDerivativeCalculation>(&inp, "dhdl-derivatives", wi);
2315 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2316 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2318 /* Non-equilibrium MD stuff */
2319 printStringNewline(&inp, "Non-equilibrium MD stuff");
2320 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2321 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2322 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2323 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2325 /* simulated tempering variables */
2326 printStringNewline(&inp, "simulated tempering variables");
2327 ir->bSimTemp = (getEnum<Boolean>(&inp, "simulated-tempering", wi) != Boolean::No);
2328 ir->simtempvals->eSimTempScale = getEnum<SimulatedTempering>(&inp, "simulated-tempering-scaling", wi);
2329 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2330 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2332 /* expanded ensemble variables */
2333 if (ir->efep == FreeEnergyPerturbationType::Expanded || ir->bSimTemp)
2335 read_expandedparams(&inp, expand, wi);
2338 /* Electric fields */
2340 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2341 gmx::KeyValueTreeTransformer transform;
2342 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2343 mdModules->initMdpTransform(transform.rules());
2344 for (const auto& path : transform.mappedPaths())
2346 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2347 mark_einp_set(inp, path[0].c_str());
2349 MdpErrorHandler errorHandler(wi);
2350 auto result = transform.transform(convertedValues, &errorHandler);
2351 ir->params = new gmx::KeyValueTreeObject(result.object());
2352 mdModules->adjustInputrecBasedOnModules(ir);
2353 errorHandler.setBackMapping(result.backMapping());
2354 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2357 /* Ion/water position swapping ("computational electrophysiology") */
2358 printStringNewline(&inp,
2359 "Ion/water position swapping for computational electrophysiology setups");
2360 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2361 ir->eSwapCoords = getEnum<SwapType>(&inp, "swapcoords", wi);
2362 if (ir->eSwapCoords != SwapType::No)
2369 printStringNoNewline(&inp, "Swap attempt frequency");
2370 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2371 printStringNoNewline(&inp, "Number of ion types to be controlled");
2372 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2375 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2377 ir->swap->ngrp = nIonTypes + static_cast<int>(SwapGroupSplittingType::Count);
2378 snew(ir->swap->grp, ir->swap->ngrp);
2379 for (i = 0; i < ir->swap->ngrp; i++)
2381 snew(ir->swap->grp[i].molname, STRLEN);
2383 printStringNoNewline(&inp,
2384 "Two index groups that contain the compartment-partitioning atoms");
2385 setStringEntry(&inp,
2387 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
2389 setStringEntry(&inp,
2391 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname,
2393 printStringNoNewline(&inp,
2394 "Use center of mass of split groups (yes/no), otherwise center of "
2395 "geometry is used");
2396 ir->swap->massw_split[0] = (getEnum<Boolean>(&inp, "massw-split0", wi) != Boolean::No);
2397 ir->swap->massw_split[1] = (getEnum<Boolean>(&inp, "massw-split1", wi) != Boolean::No);
2399 printStringNoNewline(&inp, "Name of solvent molecules");
2400 setStringEntry(&inp,
2402 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].molname,
2405 printStringNoNewline(&inp,
2406 "Split cylinder: radius, upper and lower extension (nm) (this will "
2407 "define the channels)");
2408 printStringNoNewline(&inp,
2409 "Note that the split cylinder settings do not have an influence on "
2410 "the swapping protocol,");
2411 printStringNoNewline(
2413 "however, if correctly defined, the permeation events are recorded per channel");
2414 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2415 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2416 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2417 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2418 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2419 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2421 printStringNoNewline(
2423 "Average the number of ions per compartment over these many swap attempt steps");
2424 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2426 printStringNoNewline(
2427 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2428 printStringNoNewline(
2429 &inp, "and the requested number of ions of this type in compartments A and B");
2430 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2431 for (i = 0; i < nIonTypes; i++)
2433 int ig = static_cast<int>(SwapGroupSplittingType::Count) + i;
2435 sprintf(buf, "iontype%d-name", i);
2436 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2437 sprintf(buf, "iontype%d-in-A", i);
2438 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2439 sprintf(buf, "iontype%d-in-B", i);
2440 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2443 printStringNoNewline(
2445 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2446 printStringNoNewline(
2448 "at maximum distance (= bulk concentration) to the split group layers. However,");
2449 printStringNoNewline(&inp,
2450 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2451 "layer from the middle at 0.0");
2452 printStringNoNewline(&inp,
2453 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2454 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2455 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2456 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2457 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2459 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2462 printStringNoNewline(
2463 &inp, "Start to swap ions if threshold difference to requested count is reached");
2464 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2467 /* AdResS is no longer supported, but we need grompp to be able to
2468 refuse to process old .mdp files that used it. */
2469 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2471 /* User defined thingies */
2472 printStringNewline(&inp, "User defined thingies");
2473 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2474 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2475 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2476 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2477 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2478 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2479 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2480 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2481 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2482 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.
3260 auto names = gmx::splitString(val);
3261 if (names.size() % 2 != 0)
3263 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3265 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3267 for (size_t i = 0; i < names.size() / 2; i++)
3269 // TODO this needs to be replaced by a solution using std::find_if
3273 names[2 * i].c_str(),
3274 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3280 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3285 names[2 * i + 1].c_str(),
3286 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3292 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3294 if ((j < nr) && (k < nr))
3296 ir->opts.egp_flags[nr * j + k] |= flag;
3297 ir->opts.egp_flags[nr * k + j] |= flag;
3306 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3308 int ig = -1, i = 0, gind;
3312 /* Just a quick check here, more thorough checks are in mdrun */
3313 if (strcmp(swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
3314 swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname)
3318 "The split groups can not both be '%s'.",
3319 swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname);
3322 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3323 for (ig = 0; ig < swap->ngrp; ig++)
3325 swapg = &swap->grp[ig];
3326 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3327 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3332 "%s group '%s' contains %d atoms.\n",
3333 ig < 3 ? enumValueToString(static_cast<SwapGroupSplittingType>(ig)) : "Swap",
3334 swap->grp[ig].molname,
3336 snew(swapg->ind, swapg->nat);
3337 for (i = 0; i < swapg->nat; i++)
3339 swapg->ind[i] = grps->a[grps->index[gind] + i];
3344 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3350 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3355 ig = search_string(IMDgname, grps->nr, gnames);
3356 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3358 if (IMDgroup->nat > 0)
3361 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3365 snew(IMDgroup->ind, IMDgroup->nat);
3366 for (i = 0; i < IMDgroup->nat; i++)
3368 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3373 /* Checks whether atoms are both part of a COM removal group and frozen.
3374 * If a fully frozen atom is part of a COM removal group, it is removed
3375 * from the COM removal group. A note is issued if such atoms are present.
3376 * A warning is issued for atom with one or two dimensions frozen that
3377 * are part of a COM removal group (mdrun would need to compute COM mass
3378 * per dimension to handle this correctly).
3379 * Also issues a warning when non-frozen atoms are not part of a COM
3380 * removal group while COM removal is active.
3382 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3384 const t_grpopts& opts,
3387 const int vcmRestGroup =
3388 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3390 int numFullyFrozenVcmAtoms = 0;
3391 int numPartiallyFrozenVcmAtoms = 0;
3392 int numNonVcmAtoms = 0;
3393 for (int a = 0; a < numAtoms; a++)
3395 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3396 int numFrozenDims = 0;
3397 for (int d = 0; d < DIM; d++)
3399 numFrozenDims += opts.nFreeze[freezeGroup][d];
3402 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3403 if (vcmGroup < vcmRestGroup)
3405 if (numFrozenDims == DIM)
3407 /* Do not remove COM motion for this fully frozen atom */
3408 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3410 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3413 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3414 numFullyFrozenVcmAtoms++;
3416 else if (numFrozenDims > 0)
3418 numPartiallyFrozenVcmAtoms++;
3421 else if (numFrozenDims < DIM)
3427 if (numFullyFrozenVcmAtoms > 0)
3429 std::string warningText = gmx::formatString(
3430 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3431 "removing these atoms from the COMM removal group(s)",
3432 numFullyFrozenVcmAtoms);
3433 warning_note(wi, warningText.c_str());
3435 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3437 std::string warningText = gmx::formatString(
3438 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3439 "removal group(s), due to limitations in the code these still contribute to the "
3440 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3442 numPartiallyFrozenVcmAtoms,
3444 warning(wi, warningText.c_str());
3446 if (numNonVcmAtoms > 0)
3448 std::string warningText = gmx::formatString(
3449 "%d atoms are not part of any center of mass motion removal group.\n"
3450 "This may lead to artifacts.\n"
3451 "In most cases one should use one group for the whole system.",
3453 warning(wi, warningText.c_str());
3457 void do_index(const char* mdparin,
3461 const gmx::MDModulesNotifiers& mdModulesNotifiers,
3465 t_blocka* defaultIndexGroups;
3473 int i, j, k, restnm;
3474 bool bExcl, bTable, bAnneal;
3475 char warn_buf[STRLEN];
3479 fprintf(stderr, "processing index file...\n");
3483 snew(defaultIndexGroups, 1);
3484 snew(defaultIndexGroups->index, 1);
3486 atoms_all = gmx_mtop_global_atoms(*mtop);
3487 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3488 done_atom(&atoms_all);
3492 defaultIndexGroups = init_index(ndx, &gnames);
3495 SimulationGroups* groups = &mtop->groups;
3496 natoms = mtop->natoms;
3497 symtab = &mtop->symtab;
3499 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3501 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3503 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3504 restnm = groups->groupNames.size() - 1;
3505 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3506 srenew(gnames, defaultIndexGroups->nr + 1);
3507 gnames[restnm] = *(groups->groupNames.back());
3509 set_warning_line(wi, mdparin, -1);
3511 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3512 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3513 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3514 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3515 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3518 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3520 temperatureCouplingGroupNames.size(),
3521 temperatureCouplingReferenceValues.size(),
3522 temperatureCouplingTauValues.size());
3525 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3526 do_numbering(natoms,
3528 temperatureCouplingGroupNames,
3531 SimulationAtomGroupType::TemperatureCoupling,
3533 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST,
3536 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3538 snew(ir->opts.nrdf, nr);
3539 snew(ir->opts.tau_t, nr);
3540 snew(ir->opts.ref_t, nr);
3541 if (ir->eI == IntegrationAlgorithm::BD && ir->bd_fric == 0)
3543 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3546 if (useReferenceTemperature)
3548 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3550 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3554 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3555 for (i = 0; (i < nr); i++)
3557 if ((ir->eI == IntegrationAlgorithm::BD) && ir->opts.tau_t[i] <= 0)
3560 "With integrator %s tau-t should be larger than 0",
3561 enumValueToString(ir->eI));
3562 warning_error(wi, warn_buf);
3565 if (ir->etc != TemperatureCoupling::VRescale && ir->opts.tau_t[i] == 0)
3569 "tau-t = -1 is the value to signal that a group should not have "
3570 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3573 if (ir->opts.tau_t[i] >= 0)
3575 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3578 if (ir->etc != TemperatureCoupling::No && ir->nsttcouple == -1)
3580 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3585 if ((ir->etc == TemperatureCoupling::NoseHoover) && (ir->epc == PressureCoupling::Berendsen))
3588 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3589 "md-vv; use either vrescale temperature with berendsen pressure or "
3590 "Nose-Hoover temperature with MTTK pressure");
3592 if (ir->epc == PressureCoupling::Mttk)
3594 if (ir->etc != TemperatureCoupling::NoseHoover)
3597 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3602 if (ir->nstpcouple != ir->nsttcouple)
3604 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3605 ir->nstpcouple = ir->nsttcouple = mincouple;
3607 "for current Trotter decomposition methods with vv, nsttcouple and "
3608 "nstpcouple must be equal. Both have been reset to "
3609 "min(nsttcouple,nstpcouple) = %d",
3611 warning_note(wi, warn_buf);
3616 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3617 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3619 if (ETC_ANDERSEN(ir->etc))
3621 if (ir->nsttcouple != 1)
3625 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3626 "need for larger nsttcouple > 1, since no global parameters are computed. "
3627 "nsttcouple has been reset to 1");
3628 warning_note(wi, warn_buf);
3631 nstcmin = tcouple_min_integration_steps(ir->etc);
3634 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3637 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3638 "least %d times larger than nsttcouple*dt (%g)",
3639 enumValueToString(ir->etc),
3642 ir->nsttcouple * ir->delta_t);
3643 warning(wi, warn_buf);
3646 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3647 for (i = 0; (i < nr); i++)
3649 if (ir->opts.ref_t[i] < 0)
3651 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3654 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3655 if we are in this conditional) if mc_temp is negative */
3656 if (ir->expandedvals->mc_temp < 0)
3658 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3662 /* Simulated annealing for each group. There are nr groups */
3663 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3664 if (simulatedAnnealingGroupNames.size() == 1
3665 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3667 simulatedAnnealingGroupNames.resize(0);
3669 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3672 "Wrong number of annealing values: %zu (for %d groups)\n",
3673 simulatedAnnealingGroupNames.size(),
3678 snew(ir->opts.annealing, nr);
3679 snew(ir->opts.anneal_npoints, nr);
3680 snew(ir->opts.anneal_time, nr);
3681 snew(ir->opts.anneal_temp, nr);
3682 for (i = 0; i < nr; i++)
3684 ir->opts.annealing[i] = SimulatedAnnealing::No;
3685 ir->opts.anneal_npoints[i] = 0;
3686 ir->opts.anneal_time[i] = nullptr;
3687 ir->opts.anneal_temp[i] = nullptr;
3689 if (!simulatedAnnealingGroupNames.empty())
3692 for (i = 0; i < nr; i++)
3694 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3696 ir->opts.annealing[i] = SimulatedAnnealing::No;
3698 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3700 ir->opts.annealing[i] = SimulatedAnnealing::Single;
3703 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3705 ir->opts.annealing[i] = SimulatedAnnealing::Periodic;
3711 /* Read the other fields too */
3712 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3713 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3716 "Found %zu annealing-npoints values for %zu groups\n",
3717 simulatedAnnealingPoints.size(),
3718 simulatedAnnealingGroupNames.size());
3720 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3721 size_t numSimulatedAnnealingFields = 0;
3722 for (i = 0; i < nr; i++)
3724 if (ir->opts.anneal_npoints[i] == 1)
3728 "Please specify at least a start and an end point for annealing\n");
3730 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3731 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3732 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3735 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3737 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3740 "Found %zu annealing-time values, wanted %zu\n",
3741 simulatedAnnealingTimes.size(),
3742 numSimulatedAnnealingFields);
3744 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3745 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3748 "Found %zu annealing-temp values, wanted %zu\n",
3749 simulatedAnnealingTemperatures.size(),
3750 numSimulatedAnnealingFields);
3753 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3754 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3755 convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
3757 simulatedAnnealingTemperatures,
3759 allSimulatedAnnealingTemperatures.data());
3760 for (i = 0, k = 0; i < nr; i++)
3762 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3764 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3765 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3768 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3770 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3776 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3779 "Annealing timepoints out of order: t=%f comes after "
3781 ir->opts.anneal_time[i][j],
3782 ir->opts.anneal_time[i][j - 1]);
3785 if (ir->opts.anneal_temp[i][j] < 0)
3788 "Found negative temperature in annealing: %f\n",
3789 ir->opts.anneal_temp[i][j]);
3794 /* Print out some summary information, to make sure we got it right */
3795 for (i = 0; i < nr; i++)
3797 if (ir->opts.annealing[i] != SimulatedAnnealing::No)
3799 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3801 "Simulated annealing for group %s: %s, %d timepoints\n",
3802 *(groups->groupNames[j]),
3803 enumValueToString(ir->opts.annealing[i]),
3804 ir->opts.anneal_npoints[i]);
3805 fprintf(stderr, "Time (ps) Temperature (K)\n");
3806 /* All terms except the last one */
3807 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3811 ir->opts.anneal_time[i][j],
3812 ir->opts.anneal_temp[i][j]);
3815 /* Finally the last one */
3816 j = ir->opts.anneal_npoints[i] - 1;
3817 if (ir->opts.annealing[i] == SimulatedAnnealing::Single)
3821 ir->opts.anneal_time[i][j],
3822 ir->opts.anneal_temp[i][j]);
3828 ir->opts.anneal_time[i][j],
3829 ir->opts.anneal_temp[i][j]);
3830 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3833 "There is a temperature jump when your annealing "
3845 for (int i = 1; i < ir->pull->ngroup; i++)
3847 const int gid = search_string(
3848 inputrecStrings->pullGroupNames[i].c_str(), defaultIndexGroups->nr, gnames);
3849 GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
3850 atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
3853 process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3855 checkPullCoords(ir->pull->group, ir->pull->coord);
3860 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3863 if (ir->eSwapCoords != SwapType::No)
3865 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3868 /* Make indices for IMD session */
3871 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3874 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3875 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3876 mdModulesNotifiers.preProcessingNotifier_.notify(defaultIndexGroupsAndNames);
3878 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3879 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3880 if (freezeDims.size() != DIM * freezeGroupNames.size())
3883 "Invalid Freezing input: %zu groups and %zu freeze values",
3884 freezeGroupNames.size(),
3887 do_numbering(natoms,
3892 SimulationAtomGroupType::Freeze,
3897 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3898 ir->opts.ngfrz = nr;
3899 snew(ir->opts.nFreeze, nr);
3900 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3902 for (j = 0; (j < DIM); j++, k++)
3904 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3905 if (!ir->opts.nFreeze[i][j])
3907 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3910 "Please use Y(ES) or N(O) for freezedim only "
3912 freezeDims[k].c_str());
3913 warning(wi, warn_buf);
3918 for (; (i < nr); i++)
3920 for (j = 0; (j < DIM); j++)
3922 ir->opts.nFreeze[i][j] = 0;
3926 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3927 do_numbering(natoms,
3932 SimulationAtomGroupType::EnergyOutput,
3937 add_wall_energrps(groups, ir->nwall, symtab);
3938 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3939 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3940 do_numbering(natoms,
3945 SimulationAtomGroupType::MassCenterVelocityRemoval,
3947 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART,
3951 if (ir->comm_mode != ComRemovalAlgorithm::No)
3953 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3956 /* Now we have filled the freeze struct, so we can calculate NRDF */
3957 calc_nrdf(mtop, ir, gnames);
3959 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3960 do_numbering(natoms,
3965 SimulationAtomGroupType::User1,
3970 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3971 do_numbering(natoms,
3976 SimulationAtomGroupType::User2,
3981 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3982 do_numbering(natoms,
3984 compressedXGroupNames,
3987 SimulationAtomGroupType::CompressedPositionOutput,
3992 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3993 do_numbering(natoms,
3995 orirefFitGroupNames,
3998 SimulationAtomGroupType::OrientationRestraintsFit,
4004 /* MiMiC QMMM input processing */
4005 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
4006 if (qmGroupNames.size() > 1)
4008 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
4010 /* group rest, if any, is always MM! */
4011 do_numbering(natoms,
4016 SimulationAtomGroupType::QuantumMechanics,
4021 ir->opts.ngQM = qmGroupNames.size();
4023 /* end of MiMiC QMMM input */
4027 for (auto group : gmx::keysOf(groups->groups))
4029 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
4030 for (const auto& entry : groups->groups[group])
4032 fprintf(stderr, " %s", *(groups->groupNames[entry]));
4034 fprintf(stderr, "\n");
4038 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
4039 snew(ir->opts.egp_flags, nr * nr);
4041 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
4042 if (bExcl && ir->cutoff_scheme == CutoffScheme::Verlet)
4044 warning_error(wi, "Energy group exclusions are currently not supported");
4046 if (bExcl && EEL_FULL(ir->coulombtype))
4048 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
4051 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
4052 if (bTable && !(ir->vdwtype == VanDerWaalsType::User)
4053 && !(ir->coulombtype == CoulombInteractionType::User)
4054 && !(ir->coulombtype == CoulombInteractionType::PmeUser)
4055 && !(ir->coulombtype == CoulombInteractionType::PmeUserSwitch))
4058 "Can only have energy group pair tables in combination with user tables for VdW "
4062 /* final check before going out of scope if simulated tempering variables
4063 * need to be set to default values.
4065 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
4067 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
4070 "the value for nstexpanded was not specified for "
4071 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
4072 "by default, but it is recommended to set it to an explicit value!",
4073 ir->expandedvals->nstexpanded));
4075 for (i = 0; (i < defaultIndexGroups->nr); i++)
4080 done_blocka(defaultIndexGroups);
4081 sfree(defaultIndexGroups);
4085 static void check_disre(const gmx_mtop_t& mtop)
4087 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
4089 const gmx_ffparams_t& ffparams = mtop.ffparams;
4092 for (int i = 0; i < ffparams.numTypes(); i++)
4094 int ftype = ffparams.functype[i];
4095 if (ftype == F_DISRES)
4097 int label = ffparams.iparams[i].disres.label;
4098 if (label == old_label)
4100 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
4109 "Found %d double distance restraint indices,\n"
4110 "probably the parameters for multiple pairs in one restraint "
4111 "are not identical\n",
4117 //! Returns whether dimensions have an absolute reference due to walls, pbc or freezing
4118 static BasicVector<bool> haveAbsoluteReference(const t_inputrec& ir)
4120 BasicVector<bool> absRef = { false, false, false };
4122 /* Check the degrees of freedom of the COM (not taking COMM removal into account) */
4123 for (int d = 0; d < DIM; d++)
4125 absRef[d] = (d >= ndof_com(&ir));
4127 /* Check for freeze groups */
4128 for (int g = 0; g < ir.opts.ngfrz; g++)
4130 for (int d = 0; d < DIM; d++)
4132 if (ir.opts.nFreeze[g][d] != 0)
4142 //! Returns whether position restraints are used for dimensions
4143 static BasicVector<bool> havePositionRestraints(const gmx_mtop_t& sys)
4145 BasicVector<bool> havePosres = { false, false, false };
4147 for (const auto ilists : IListRange(sys))
4149 const auto& posResList = ilists.list()[F_POSRES];
4150 const auto& fbPosResList = ilists.list()[F_FBPOSRES];
4151 if (ilists.nmol() > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
4153 for (int i = 0; i < posResList.size(); i += 2)
4155 const t_iparams& pr = sys.ffparams.iparams[posResList.iatoms[i]];
4156 for (int d = 0; d < DIM; d++)
4158 if (pr.posres.fcA[d] != 0)
4160 havePosres[d] = true;
4164 for (int i = 0; i < fbPosResList.size(); i += 2)
4166 /* Check for flat-bottom posres */
4167 const t_iparams& pr = sys.ffparams.iparams[fbPosResList.iatoms[i]];
4168 if (pr.fbposres.k != 0)
4170 switch (pr.fbposres.geom)
4172 case efbposresSPHERE: havePosres = { true, true, true }; break;
4173 case efbposresCYLINDERX: havePosres[YY] = havePosres[ZZ] = true; break;
4174 case efbposresCYLINDERY: havePosres[XX] = havePosres[ZZ] = true; break;
4175 case efbposresCYLINDER:
4176 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4177 case efbposresCYLINDERZ: havePosres[XX] = havePosres[YY] = true; break;
4178 case efbposresX: /* d=XX */
4179 case efbposresY: /* d=YY */
4180 case efbposresZ: /* d=ZZ */
4181 havePosres[pr.fbposres.geom - efbposresX] = true;
4185 "Invalid geometry for flat-bottom position restraint.\n"
4186 "Expected nr between 1 and %d. Found %d\n",
4198 static void check_combination_rule_differences(const gmx_mtop_t& mtop,
4200 bool* bC6ParametersWorkWithGeometricRules,
4201 bool* bC6ParametersWorkWithLBRules,
4202 bool* bLBRulesPossible)
4204 int ntypes, tpi, tpj;
4207 double c6i, c6j, c12i, c12j;
4208 double c6, c6_geometric, c6_LB;
4209 double sigmai, sigmaj, epsi, epsj;
4210 bool bCanDoLBRules, bCanDoGeometricRules;
4213 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4214 * force-field floating point parameters.
4217 ptr = getenv("GMX_LJCOMB_TOL");
4221 double gmx_unused canary;
4223 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4226 FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4231 *bC6ParametersWorkWithLBRules = TRUE;
4232 *bC6ParametersWorkWithGeometricRules = TRUE;
4233 bCanDoLBRules = TRUE;
4234 ntypes = mtop.ffparams.atnr;
4235 snew(typecount, ntypes);
4236 gmx_mtop_count_atomtypes(mtop, state, typecount);
4237 *bLBRulesPossible = TRUE;
4238 for (tpi = 0; tpi < ntypes; ++tpi)
4240 c6i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4241 c12i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4242 for (tpj = tpi; tpj < ntypes; ++tpj)
4244 c6j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4245 c12j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4246 c6 = mtop.ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4247 c6_geometric = std::sqrt(c6i * c6j);
4248 if (!gmx_numzero(c6_geometric))
4250 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4252 sigmai = gmx::sixthroot(c12i / c6i);
4253 sigmaj = gmx::sixthroot(c12j / c6j);
4254 epsi = c6i * c6i / (4.0 * c12i);
4255 epsj = c6j * c6j / (4.0 * c12j);
4256 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4260 *bLBRulesPossible = FALSE;
4261 c6_LB = c6_geometric;
4263 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4268 *bC6ParametersWorkWithLBRules = FALSE;
4271 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4273 if (!bCanDoGeometricRules)
4275 *bC6ParametersWorkWithGeometricRules = FALSE;
4282 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t& mtop, warninp_t wi)
4284 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4286 check_combination_rule_differences(
4287 mtop, 0, &bC6ParametersWorkWithGeometricRules, &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4288 if (ir->ljpme_combination_rule == LongRangeVdW::LB)
4290 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4293 "You are using arithmetic-geometric combination rules "
4294 "in LJ-PME, but your non-bonded C6 parameters do not "
4295 "follow these rules.");
4300 if (!bC6ParametersWorkWithGeometricRules)
4302 if (ir->eDispCorr != DispersionCorrectionType::No)
4305 "You are using geometric combination rules in "
4306 "LJ-PME, but your non-bonded C6 parameters do "
4307 "not follow these rules. "
4308 "This will introduce very small errors in the forces and energies in "
4309 "your simulations. Dispersion correction will correct total energy "
4310 "and/or pressure for isotropic systems, but not forces or surface "
4316 "You are using geometric combination rules in "
4317 "LJ-PME, but your non-bonded C6 parameters do "
4318 "not follow these rules. "
4319 "This will introduce very small errors in the forces and energies in "
4320 "your simulations. If your system is homogeneous, consider using "
4321 "dispersion correction "
4322 "for the total energy and pressure.");
4328 static bool allTrue(const BasicVector<bool>& boolVector)
4330 return boolVector[0] && boolVector[1] && boolVector[2];
4333 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4335 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4336 GMX_ASSERT(gmx::checkMtsRequirements(*ir).empty(), "All MTS requirements should be met here");
4338 char err_buf[STRLEN];
4341 gmx_mtop_atomloop_block_t aloopb;
4342 char warn_buf[STRLEN];
4344 set_warning_line(wi, mdparin, -1);
4346 if (ir->comm_mode != ComRemovalAlgorithm::No && allTrue(havePositionRestraints(*sys)))
4349 "Removing center of mass motion in the presence of position restraints might "
4350 "cause artifacts. When you are using position restraints to equilibrate a "
4351 "macro-molecule, the artifacts are usually negligible.");
4354 if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0 && ir->nstlist > 1
4355 && ((EI_MD(ir->eI) || EI_SD(ir->eI))
4356 && (ir->etc == TemperatureCoupling::VRescale || ir->etc == TemperatureCoupling::Berendsen)))
4358 /* Check if a too small Verlet buffer might potentially
4359 * cause more drift than the thermostat can couple off.
4361 /* Temperature error fraction for warning and suggestion */
4362 const real T_error_warn = 0.002;
4363 const real T_error_suggest = 0.001;
4364 /* For safety: 2 DOF per atom (typical with constraints) */
4365 const real nrdf_at = 2;
4366 real T, tau, max_T_error;
4371 for (i = 0; i < ir->opts.ngtc; i++)
4373 T = std::max(T, ir->opts.ref_t[i]);
4374 tau = std::max(tau, ir->opts.tau_t[i]);
4378 /* This is a worst case estimate of the temperature error,
4379 * assuming perfect buffer estimation and no cancelation
4380 * of errors. The factor 0.5 is because energy distributes
4381 * equally over Ekin and Epot.
4383 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * gmx::c_boltz * T);
4384 if (max_T_error > T_error_warn)
4387 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4388 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4389 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4390 "%.0e or decrease tau_t.",
4395 100 * T_error_suggest,
4396 ir->verletbuf_tol * T_error_suggest / max_T_error);
4397 warning(wi, warn_buf);
4402 if (ETC_ANDERSEN(ir->etc))
4406 for (i = 0; i < ir->opts.ngtc; i++)
4409 "all tau_t must currently be equal using Andersen temperature control, "
4410 "violated for group %d",
4412 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4414 "all tau_t must be positive using Andersen temperature control, "
4418 CHECK(ir->opts.tau_t[i] < 0);
4421 if (ir->etc == TemperatureCoupling::AndersenMassive && ir->comm_mode != ComRemovalAlgorithm::No)
4423 for (i = 0; i < ir->opts.ngtc; i++)
4425 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4427 "tau_t/delta_t for group %d for temperature control method %s must be a "
4428 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4429 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4432 enumValueToString(ir->etc),
4436 CHECK(nsteps % ir->nstcomm != 0);
4441 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != IntegrationAlgorithm::BD
4442 && ir->comm_mode == ComRemovalAlgorithm::No
4443 && !(allTrue(haveAbsoluteReference(*ir)) || allTrue(havePositionRestraints(*sys)) || ir->nsteps <= 10)
4444 && !ETC_ANDERSEN(ir->etc))
4447 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4448 "rounding errors can lead to build up of kinetic energy of the center of mass");
4451 if (ir->epc == PressureCoupling::ParrinelloRahman && ir->etc == TemperatureCoupling::NoseHoover)
4454 for (int g = 0; g < ir->opts.ngtc; g++)
4456 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4458 if (ir->tau_p < 1.9 * tau_t_max)
4460 std::string message = gmx::formatString(
4461 "With %s T-coupling and %s p-coupling, "
4462 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4463 enumValueToString(ir->etc),
4464 enumValueToString(ir->epc),
4469 warning(wi, message.c_str());
4473 /* Check for pressure coupling with absolute position restraints */
4474 if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == RefCoordScaling::No)
4476 const BasicVector<bool> havePosres = havePositionRestraints(*sys);
4478 for (m = 0; m < DIM; m++)
4480 if (havePosres[m] && norm2(ir->compress[m]) > 0)
4483 "You are using pressure coupling with absolute position restraints, "
4484 "this will give artifacts. Use the refcoord_scaling option.");
4492 aloopb = gmx_mtop_atomloop_block_init(*sys);
4494 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4496 if (atom->q != 0 || atom->qB != 0)
4504 if (EEL_FULL(ir->coulombtype))
4507 "You are using full electrostatics treatment %s for a system without charges.\n"
4508 "This costs a lot of performance for just processing zeros, consider using %s "
4510 enumValueToString(ir->coulombtype),
4511 enumValueToString(CoulombInteractionType::Cut));
4512 warning(wi, err_buf);
4517 if (ir->coulombtype == CoulombInteractionType::Cut && ir->rcoulomb > 0)
4520 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4521 "You might want to consider using %s electrostatics.\n",
4522 enumValueToString(CoulombInteractionType::Pme));
4523 warning_note(wi, err_buf);
4527 /* Check if combination rules used in LJ-PME are the same as in the force field */
4528 if (EVDW_PME(ir->vdwtype))
4530 check_combination_rules(ir, *sys, wi);
4533 /* Generalized reaction field */
4534 if (ir->coulombtype == CoulombInteractionType::GRFNotused)
4537 "Generalized reaction-field electrostatics is no longer supported. "
4538 "You can use normal reaction-field instead and compute the reaction-field "
4539 "constant by hand.");
4542 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->sc_alpha != 0
4543 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4545 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4553 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4555 if (ir->pull->coord[i].eGeom != PullGroupGeometry::Transformation
4556 && (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);