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);
2486 gmx::TextOutputFile stream(mdparout);
2487 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2489 // Transform module data into a flat key-value tree for output.
2490 gmx::KeyValueTreeBuilder builder;
2491 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2492 mdModules->buildMdpOutput(&builderObject);
2494 gmx::TextWriter writer(&stream);
2495 writeKeyValueTreeAsMdp(&writer, builder.build());
2500 /* Process options if necessary */
2501 for (m = 0; m < 2; m++)
2503 for (i = 0; i < 2 * DIM; i++)
2507 if (ir->epc != PressureCoupling::No)
2511 case PressureCouplingType::Isotropic:
2512 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2516 "Pressure coupling incorrect number of values (I need exactly 1)");
2518 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2520 case PressureCouplingType::SemiIsotropic:
2521 case PressureCouplingType::SurfaceTension:
2522 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2526 "Pressure coupling incorrect number of values (I need exactly 2)");
2528 dumdub[m][YY] = dumdub[m][XX];
2530 case PressureCouplingType::Anisotropic:
2531 if (sscanf(dumstr[m],
2532 "%lf%lf%lf%lf%lf%lf",
2543 "Pressure coupling incorrect number of values (I need exactly 6)");
2548 "Pressure coupling type %s not implemented yet",
2549 enumValueToString(ir->epct));
2553 clear_mat(ir->ref_p);
2554 clear_mat(ir->compress);
2555 for (i = 0; i < DIM; i++)
2557 ir->ref_p[i][i] = dumdub[1][i];
2558 ir->compress[i][i] = dumdub[0][i];
2560 if (ir->epct == PressureCouplingType::Anisotropic)
2562 ir->ref_p[XX][YY] = dumdub[1][3];
2563 ir->ref_p[XX][ZZ] = dumdub[1][4];
2564 ir->ref_p[YY][ZZ] = dumdub[1][5];
2565 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2568 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2569 "apply a threefold shear stress?\n");
2571 ir->compress[XX][YY] = dumdub[0][3];
2572 ir->compress[XX][ZZ] = dumdub[0][4];
2573 ir->compress[YY][ZZ] = dumdub[0][5];
2574 for (i = 0; i < DIM; i++)
2576 for (m = 0; m < i; m++)
2578 ir->ref_p[i][m] = ir->ref_p[m][i];
2579 ir->compress[i][m] = ir->compress[m][i];
2584 if (ir->comm_mode == ComRemovalAlgorithm::No)
2589 opts->couple_moltype = nullptr;
2590 if (strlen(inputrecStrings->couple_moltype) > 0)
2592 if (ir->efep != FreeEnergyPerturbationType::No)
2594 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2595 if (opts->couple_lam0 == opts->couple_lam1)
2597 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2599 if (ir->eI == IntegrationAlgorithm::MD
2600 && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2604 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2611 "Free energy is turned off, so we will not decouple the molecule listed "
2615 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2616 if (ir->efep != FreeEnergyPerturbationType::No)
2618 if (fep->delta_lambda != 0)
2620 ir->efep = FreeEnergyPerturbationType::SlowGrowth;
2624 if (fep->edHdLPrintEnergy == FreeEnergyPrintEnergy::Yes)
2626 fep->edHdLPrintEnergy = FreeEnergyPrintEnergy::Total;
2628 "Old option for dhdl-print-energy given: "
2629 "changing \"yes\" to \"total\"\n");
2632 if (ir->bSimTemp && (fep->edHdLPrintEnergy == FreeEnergyPrintEnergy::No))
2634 /* always print out the energy to dhdl if we are doing
2635 expanded ensemble, since we need the total energy for
2636 analysis if the temperature is changing. In some
2637 conditions one may only want the potential energy, so
2638 we will allow that if the appropriate mdp setting has
2639 been enabled. Otherwise, total it is:
2641 fep->edHdLPrintEnergy = FreeEnergyPrintEnergy::Total;
2644 if ((ir->efep != FreeEnergyPerturbationType::No) || ir->bSimTemp)
2646 ir->bExpanded = FALSE;
2647 if ((ir->efep == FreeEnergyPerturbationType::Expanded) || ir->bSimTemp)
2649 ir->bExpanded = TRUE;
2651 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2652 if (ir->bSimTemp) /* done after fep params */
2654 do_simtemp_params(ir);
2657 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2658 * setup and not on the old way of specifying the free-energy setup,
2659 * we should check for using soft-core when not needed, since that
2660 * can complicate the sampling significantly.
2661 * Note that we only check for the automated coupling setup.
2662 * If the (advanced) user does FEP through manual topology changes,
2663 * this check will not be triggered.
2665 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->n_lambda == 0
2666 && ir->fepvals->sc_alpha != 0
2667 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2670 "You are using soft-core interactions while the Van der Waals interactions are "
2671 "not decoupled (note that the sc-coul option is only active when using lambda "
2672 "states). Although this will not lead to errors, you will need much more "
2673 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2678 ir->fepvals->n_lambda = 0;
2681 /* WALL PARAMETERS */
2683 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2685 /* ORIENTATION RESTRAINT PARAMETERS */
2687 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2689 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2692 /* DEFORMATION PARAMETERS */
2694 clear_mat(ir->deform);
2695 for (i = 0; i < 6; i++)
2700 double gmx_unused canary;
2701 int ndeform = sscanf(inputrecStrings->deform,
2702 "%lf %lf %lf %lf %lf %lf %lf",
2711 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2715 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2716 inputrecStrings->deform)
2719 for (i = 0; i < 3; i++)
2721 ir->deform[i][i] = dumdub[0][i];
2723 ir->deform[YY][XX] = dumdub[0][3];
2724 ir->deform[ZZ][XX] = dumdub[0][4];
2725 ir->deform[ZZ][YY] = dumdub[0][5];
2726 if (ir->epc != PressureCoupling::No)
2728 for (i = 0; i < 3; i++)
2730 for (j = 0; j <= i; j++)
2732 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2734 warning_error(wi, "A box element has deform set and compressibility > 0");
2738 for (i = 0; i < 3; i++)
2740 for (j = 0; j < i; j++)
2742 if (ir->deform[i][j] != 0)
2744 for (m = j; m < DIM; m++)
2746 if (ir->compress[m][j] != 0)
2749 "An off-diagonal box element has deform set while "
2750 "compressibility > 0 for the same component of another box "
2751 "vector, this might lead to spurious periodicity effects.");
2752 warning(wi, warn_buf);
2760 /* Ion/water position swapping checks */
2761 if (ir->eSwapCoords != SwapType::No)
2763 if (ir->swap->nstswap < 1)
2765 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2767 if (ir->swap->nAverage < 1)
2769 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2771 if (ir->swap->threshold < 1.0)
2773 warning_error(wi, "Ion count threshold must be at least 1.\n");
2777 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2780 std::vector<std::string> errorMessages;
2781 ir->mtsLevels = gmx::setupMtsLevels(opts->mtsOpts, &errorMessages);
2783 for (const auto& errorMessage : errorMessages)
2785 warning_error(wi, errorMessage.c_str());
2791 gmx::checkAwhParams(*ir->awhParams, *ir, wi);
2798 /* We would like gn to be const as well, but C doesn't allow this */
2799 /* TODO this is utility functionality (search for the index of a
2800 string in a collection), so should be refactored and located more
2802 int search_string(const char* s, int ng, char* gn[])
2806 for (i = 0; (i < ng); i++)
2808 if (gmx_strcasecmp(s, gn[i]) == 0)
2815 "Group %s referenced in the .mdp file was not found in the index file.\n"
2816 "Group names must match either [moleculetype] names or custom index group\n"
2817 "names, in which case you must supply an index file to the '-n' option\n"
2822 static void atomGroupRangeValidation(int natoms, int groupIndex, const t_blocka& block)
2824 /* Now go over the atoms in the group */
2825 for (int j = block.index[groupIndex]; (j < block.index[groupIndex + 1]); j++)
2827 int aj = block.a[j];
2829 /* Range checking */
2830 if ((aj < 0) || (aj >= natoms))
2832 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2837 static void do_numbering(int natoms,
2838 SimulationGroups* groups,
2839 gmx::ArrayRef<std::string> groupsFromMdpFile,
2842 SimulationAtomGroupType gtype,
2848 unsigned short* cbuf;
2849 AtomGroupIndices* grps = &(groups->groups[gtype]);
2852 char warn_buf[STRLEN];
2854 title = shortName(gtype);
2857 /* Mark all id's as not set */
2858 for (int i = 0; (i < natoms); i++)
2863 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2865 /* Lookup the group name in the block structure */
2866 const int gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2867 if ((grptp != egrptpONE) || (i == 0))
2869 grps->emplace_back(gid);
2871 GMX_ASSERT(block, "Can't have a nullptr block");
2872 atomGroupRangeValidation(natoms, gid, *block);
2873 /* Now go over the atoms in the group */
2874 for (int j = block->index[gid]; (j < block->index[gid + 1]); j++)
2876 const int aj = block->a[j];
2877 /* Lookup up the old group number */
2878 const int ognr = cbuf[aj];
2881 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title, ognr + 1, i + 1);
2885 /* Store the group number in buffer */
2886 if (grptp == egrptpONE)
2899 /* Now check whether we have done all atoms */
2902 if (grptp == egrptpALL)
2904 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2906 else if (grptp == egrptpPART)
2908 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2909 warning_note(wi, warn_buf);
2911 /* Assign all atoms currently unassigned to a rest group */
2912 for (int j = 0; (j < natoms); j++)
2914 if (cbuf[j] == NOGID)
2916 cbuf[j] = grps->size();
2919 if (grptp != egrptpPART)
2923 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title, natoms - ntot);
2925 /* Add group name "rest" */
2926 grps->emplace_back(restnm);
2928 /* Assign the rest name to all atoms not currently assigned to a group */
2929 for (int j = 0; (j < natoms); j++)
2931 if (cbuf[j] == NOGID)
2933 // group size was not updated before this here, so need to use -1.
2934 cbuf[j] = grps->size() - 1;
2940 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2942 /* All atoms are part of one (or no) group, no index required */
2943 groups->groupNumbers[gtype].clear();
2947 for (int j = 0; (j < natoms); j++)
2949 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2956 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2959 pull_params_t* pull;
2960 int natoms, imin, jmin;
2961 int * nrdf2, *na_vcm, na_tot;
2962 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2967 * First calc 3xnr-atoms for each group
2968 * then subtract half a degree of freedom for each constraint
2970 * Only atoms and nuclei contribute to the degrees of freedom...
2975 const SimulationGroups& groups = mtop->groups;
2976 natoms = mtop->natoms;
2978 /* Allocate one more for a possible rest group */
2979 /* We need to sum degrees of freedom into doubles,
2980 * since floats give too low nrdf's above 3 million atoms.
2982 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
2983 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2984 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2985 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2986 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2988 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2992 for (gmx::index i = 0;
2993 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
2997 clear_ivec(dof_vcm[i]);
2999 nrdf_vcm_sub[i] = 0;
3001 snew(nrdf2, natoms);
3002 for (const AtomProxy atomP : AtomRange(*mtop))
3004 const t_atom& local = atomP.atom();
3005 int i = atomP.globalAtomNumber();
3007 if (local.ptype == ParticleType::Atom || local.ptype == ParticleType::Nucleus)
3009 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
3010 for (int d = 0; d < DIM; d++)
3012 if (opts->nFreeze[g][d] == 0)
3014 /* Add one DOF for particle i (counted as 2*1) */
3016 /* VCM group i has dim d as a DOF */
3017 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
3021 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
3023 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
3029 for (const gmx_molblock_t& molb : mtop->molblock)
3031 const gmx_moltype_t& molt = mtop->moltype[molb.type];
3032 const t_atom* atom = molt.atoms.atom;
3033 for (int mol = 0; mol < molb.nmol; mol++)
3035 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
3037 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
3038 for (int i = 0; i < molt.ilist[ftype].size();)
3040 /* Subtract degrees of freedom for the constraints,
3041 * if the particles still have degrees of freedom left.
3042 * If one of the particles is a vsite or a shell, then all
3043 * constraint motion will go there, but since they do not
3044 * contribute to the constraints the degrees of freedom do not
3047 int ai = as + ia[i + 1];
3048 int aj = as + ia[i + 2];
3049 if (((atom[ia[i + 1]].ptype == ParticleType::Nucleus)
3050 || (atom[ia[i + 1]].ptype == ParticleType::Atom))
3051 && ((atom[ia[i + 2]].ptype == ParticleType::Nucleus)
3052 || (atom[ia[i + 2]].ptype == ParticleType::Atom)))
3070 imin = std::min(imin, nrdf2[ai]);
3071 jmin = std::min(jmin, nrdf2[aj]);
3074 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3076 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3078 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3080 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3083 i += interaction_function[ftype].nratoms + 1;
3086 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3087 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3089 /* Subtract 1 dof from every atom in the SETTLE */
3090 for (int j = 0; j < 3; j++)
3092 int ai = as + ia[i + 1 + j];
3093 imin = std::min(2, nrdf2[ai]);
3095 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3097 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3102 as += molt.atoms.nr;
3108 /* Correct nrdf for the COM constraints.
3109 * We correct using the TC and VCM group of the first atom
3110 * in the reference and pull group. If atoms in one pull group
3111 * belong to different TC or VCM groups it is anyhow difficult
3112 * to determine the optimal nrdf assignment.
3114 pull = ir->pull.get();
3116 for (int i = 0; i < pull->ncoord; i++)
3118 if (pull->coord[i].eType != PullingAlgorithm::Constraint)
3125 for (int j = 0; j < 2; j++)
3127 const t_pull_group* pgrp;
3129 pgrp = &pull->group[pull->coord[i].group[j]];
3131 if (!pgrp->ind.empty())
3133 /* Subtract 1/2 dof from each group */
3134 int ai = pgrp->ind[0];
3135 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3137 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3139 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3142 "Center of mass pulling constraints caused the number of degrees "
3143 "of freedom for temperature coupling group %s to be negative",
3144 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3145 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3150 /* We need to subtract the whole DOF from group j=1 */
3157 if (ir->nstcomm != 0)
3159 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3160 "Expect at least one group when removing COM motion");
3162 /* We remove COM motion up to dim ndof_com() */
3163 const int ndim_rm_vcm = ndof_com(ir);
3165 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3166 * the number of degrees of freedom in each vcm group when COM
3167 * translation is removed and 6 when rotation is removed as well.
3168 * Note that we do not and should not include the rest group here.
3170 for (gmx::index j = 0;
3171 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]);
3174 switch (ir->comm_mode)
3176 case ComRemovalAlgorithm::Linear:
3177 case ComRemovalAlgorithm::LinearAccelerationCorrection:
3178 nrdf_vcm_sub[j] = 0;
3179 for (int d = 0; d < ndim_rm_vcm; d++)
3187 case ComRemovalAlgorithm::Angular: nrdf_vcm_sub[j] = 6; break;
3188 default: gmx_incons("Checking comm_mode");
3192 for (gmx::index i = 0;
3193 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]);
3196 /* Count the number of atoms of TC group i for every VCM group */
3197 for (gmx::index j = 0;
3198 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3204 for (int ai = 0; ai < natoms; ai++)
3206 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3208 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3212 /* Correct for VCM removal according to the fraction of each VCM
3213 * group present in this TC group.
3215 nrdf_uc = nrdf_tc[i];
3217 for (gmx::index j = 0;
3218 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3221 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3223 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3224 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3229 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3231 opts->nrdf[i] = nrdf_tc[i];
3232 if (opts->nrdf[i] < 0)
3237 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3238 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]],
3247 sfree(nrdf_vcm_sub);
3250 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3252 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3253 * But since this is much larger than STRLEN, such a line can not be parsed.
3254 * The real maximum is the number of names that fit in a string: STRLEN/2.
3259 auto names = gmx::splitString(val);
3260 if (names.size() % 2 != 0)
3262 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3264 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3266 for (size_t i = 0; i < names.size() / 2; i++)
3268 // TODO this needs to be replaced by a solution using std::find_if
3272 names[2 * i].c_str(),
3273 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3279 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3284 names[2 * i + 1].c_str(),
3285 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3291 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3293 if ((j < nr) && (k < nr))
3295 ir->opts.egp_flags[nr * j + k] |= flag;
3296 ir->opts.egp_flags[nr * k + j] |= flag;
3305 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3307 int ig = -1, i = 0, gind;
3311 /* Just a quick check here, more thorough checks are in mdrun */
3312 if (strcmp(swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
3313 swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname)
3317 "The split groups can not both be '%s'.",
3318 swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname);
3321 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3322 for (ig = 0; ig < swap->ngrp; ig++)
3324 swapg = &swap->grp[ig];
3325 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3326 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3331 "%s group '%s' contains %d atoms.\n",
3332 ig < 3 ? enumValueToString(static_cast<SwapGroupSplittingType>(ig)) : "Swap",
3333 swap->grp[ig].molname,
3335 snew(swapg->ind, swapg->nat);
3336 for (i = 0; i < swapg->nat; i++)
3338 swapg->ind[i] = grps->a[grps->index[gind] + i];
3343 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3349 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3354 ig = search_string(IMDgname, grps->nr, gnames);
3355 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3357 if (IMDgroup->nat > 0)
3360 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3364 snew(IMDgroup->ind, IMDgroup->nat);
3365 for (i = 0; i < IMDgroup->nat; i++)
3367 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3372 /* Checks whether atoms are both part of a COM removal group and frozen.
3373 * If a fully frozen atom is part of a COM removal group, it is removed
3374 * from the COM removal group. A note is issued if such atoms are present.
3375 * A warning is issued for atom with one or two dimensions frozen that
3376 * are part of a COM removal group (mdrun would need to compute COM mass
3377 * per dimension to handle this correctly).
3378 * Also issues a warning when non-frozen atoms are not part of a COM
3379 * removal group while COM removal is active.
3381 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3383 const t_grpopts& opts,
3386 const int vcmRestGroup =
3387 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3389 int numFullyFrozenVcmAtoms = 0;
3390 int numPartiallyFrozenVcmAtoms = 0;
3391 int numNonVcmAtoms = 0;
3392 for (int a = 0; a < numAtoms; a++)
3394 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3395 int numFrozenDims = 0;
3396 for (int d = 0; d < DIM; d++)
3398 numFrozenDims += opts.nFreeze[freezeGroup][d];
3401 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3402 if (vcmGroup < vcmRestGroup)
3404 if (numFrozenDims == DIM)
3406 /* Do not remove COM motion for this fully frozen atom */
3407 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3409 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3412 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3413 numFullyFrozenVcmAtoms++;
3415 else if (numFrozenDims > 0)
3417 numPartiallyFrozenVcmAtoms++;
3420 else if (numFrozenDims < DIM)
3426 if (numFullyFrozenVcmAtoms > 0)
3428 std::string warningText = gmx::formatString(
3429 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3430 "removing these atoms from the COMM removal group(s)",
3431 numFullyFrozenVcmAtoms);
3432 warning_note(wi, warningText.c_str());
3434 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3436 std::string warningText = gmx::formatString(
3437 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3438 "removal group(s), due to limitations in the code these still contribute to the "
3439 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3441 numPartiallyFrozenVcmAtoms,
3443 warning(wi, warningText.c_str());
3445 if (numNonVcmAtoms > 0)
3447 std::string warningText = gmx::formatString(
3448 "%d atoms are not part of any center of mass motion removal group.\n"
3449 "This may lead to artifacts.\n"
3450 "In most cases one should use one group for the whole system.",
3452 warning(wi, warningText.c_str());
3456 void do_index(const char* mdparin,
3460 const gmx::MDModulesNotifiers& mdModulesNotifiers,
3464 t_blocka* defaultIndexGroups;
3472 int i, j, k, restnm;
3473 bool bExcl, bTable, bAnneal;
3474 char warn_buf[STRLEN];
3478 fprintf(stderr, "processing index file...\n");
3482 snew(defaultIndexGroups, 1);
3483 snew(defaultIndexGroups->index, 1);
3485 atoms_all = gmx_mtop_global_atoms(*mtop);
3486 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3487 done_atom(&atoms_all);
3491 defaultIndexGroups = init_index(ndx, &gnames);
3494 SimulationGroups* groups = &mtop->groups;
3495 natoms = mtop->natoms;
3496 symtab = &mtop->symtab;
3498 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3500 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3502 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3503 restnm = groups->groupNames.size() - 1;
3504 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3505 srenew(gnames, defaultIndexGroups->nr + 1);
3506 gnames[restnm] = *(groups->groupNames.back());
3508 set_warning_line(wi, mdparin, -1);
3510 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3511 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3512 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3513 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3514 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3517 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3519 temperatureCouplingGroupNames.size(),
3520 temperatureCouplingReferenceValues.size(),
3521 temperatureCouplingTauValues.size());
3524 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3525 do_numbering(natoms,
3527 temperatureCouplingGroupNames,
3530 SimulationAtomGroupType::TemperatureCoupling,
3532 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST,
3535 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3537 snew(ir->opts.nrdf, nr);
3538 snew(ir->opts.tau_t, nr);
3539 snew(ir->opts.ref_t, nr);
3540 if (ir->eI == IntegrationAlgorithm::BD && ir->bd_fric == 0)
3542 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3545 if (useReferenceTemperature)
3547 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3549 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3553 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3554 for (i = 0; (i < nr); i++)
3556 if ((ir->eI == IntegrationAlgorithm::BD) && ir->opts.tau_t[i] <= 0)
3559 "With integrator %s tau-t should be larger than 0",
3560 enumValueToString(ir->eI));
3561 warning_error(wi, warn_buf);
3564 if (ir->etc != TemperatureCoupling::VRescale && ir->opts.tau_t[i] == 0)
3568 "tau-t = -1 is the value to signal that a group should not have "
3569 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3572 if (ir->opts.tau_t[i] >= 0)
3574 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3577 if (ir->etc != TemperatureCoupling::No && ir->nsttcouple == -1)
3579 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3584 if ((ir->etc == TemperatureCoupling::NoseHoover) && (ir->epc == PressureCoupling::Berendsen))
3587 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3588 "md-vv; use either vrescale temperature with berendsen pressure or "
3589 "Nose-Hoover temperature with MTTK pressure");
3591 if (ir->epc == PressureCoupling::Mttk)
3593 if (ir->etc != TemperatureCoupling::NoseHoover)
3596 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3601 if (ir->nstpcouple != ir->nsttcouple)
3603 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3604 ir->nstpcouple = ir->nsttcouple = mincouple;
3606 "for current Trotter decomposition methods with vv, nsttcouple and "
3607 "nstpcouple must be equal. Both have been reset to "
3608 "min(nsttcouple,nstpcouple) = %d",
3610 warning_note(wi, warn_buf);
3615 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3616 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3618 if (ETC_ANDERSEN(ir->etc))
3620 if (ir->nsttcouple != 1)
3624 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3625 "need for larger nsttcouple > 1, since no global parameters are computed. "
3626 "nsttcouple has been reset to 1");
3627 warning_note(wi, warn_buf);
3630 nstcmin = tcouple_min_integration_steps(ir->etc);
3633 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3636 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3637 "least %d times larger than nsttcouple*dt (%g)",
3638 enumValueToString(ir->etc),
3641 ir->nsttcouple * ir->delta_t);
3642 warning(wi, warn_buf);
3645 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3646 for (i = 0; (i < nr); i++)
3648 if (ir->opts.ref_t[i] < 0)
3650 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3653 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3654 if we are in this conditional) if mc_temp is negative */
3655 if (ir->expandedvals->mc_temp < 0)
3657 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3661 /* Simulated annealing for each group. There are nr groups */
3662 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3663 if (simulatedAnnealingGroupNames.size() == 1
3664 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3666 simulatedAnnealingGroupNames.resize(0);
3668 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3671 "Wrong number of annealing values: %zu (for %d groups)\n",
3672 simulatedAnnealingGroupNames.size(),
3677 snew(ir->opts.annealing, nr);
3678 snew(ir->opts.anneal_npoints, nr);
3679 snew(ir->opts.anneal_time, nr);
3680 snew(ir->opts.anneal_temp, nr);
3681 for (i = 0; i < nr; i++)
3683 ir->opts.annealing[i] = SimulatedAnnealing::No;
3684 ir->opts.anneal_npoints[i] = 0;
3685 ir->opts.anneal_time[i] = nullptr;
3686 ir->opts.anneal_temp[i] = nullptr;
3688 if (!simulatedAnnealingGroupNames.empty())
3691 for (i = 0; i < nr; i++)
3693 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3695 ir->opts.annealing[i] = SimulatedAnnealing::No;
3697 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3699 ir->opts.annealing[i] = SimulatedAnnealing::Single;
3702 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3704 ir->opts.annealing[i] = SimulatedAnnealing::Periodic;
3710 /* Read the other fields too */
3711 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3712 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3715 "Found %zu annealing-npoints values for %zu groups\n",
3716 simulatedAnnealingPoints.size(),
3717 simulatedAnnealingGroupNames.size());
3719 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3720 size_t numSimulatedAnnealingFields = 0;
3721 for (i = 0; i < nr; i++)
3723 if (ir->opts.anneal_npoints[i] == 1)
3727 "Please specify at least a start and an end point for annealing\n");
3729 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3730 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3731 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3734 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3736 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3739 "Found %zu annealing-time values, wanted %zu\n",
3740 simulatedAnnealingTimes.size(),
3741 numSimulatedAnnealingFields);
3743 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3744 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3747 "Found %zu annealing-temp values, wanted %zu\n",
3748 simulatedAnnealingTemperatures.size(),
3749 numSimulatedAnnealingFields);
3752 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3753 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3754 convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
3756 simulatedAnnealingTemperatures,
3758 allSimulatedAnnealingTemperatures.data());
3759 for (i = 0, k = 0; i < nr; i++)
3761 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3763 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3764 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3767 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3769 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3775 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3778 "Annealing timepoints out of order: t=%f comes after "
3780 ir->opts.anneal_time[i][j],
3781 ir->opts.anneal_time[i][j - 1]);
3784 if (ir->opts.anneal_temp[i][j] < 0)
3787 "Found negative temperature in annealing: %f\n",
3788 ir->opts.anneal_temp[i][j]);
3793 /* Print out some summary information, to make sure we got it right */
3794 for (i = 0; i < nr; i++)
3796 if (ir->opts.annealing[i] != SimulatedAnnealing::No)
3798 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3800 "Simulated annealing for group %s: %s, %d timepoints\n",
3801 *(groups->groupNames[j]),
3802 enumValueToString(ir->opts.annealing[i]),
3803 ir->opts.anneal_npoints[i]);
3804 fprintf(stderr, "Time (ps) Temperature (K)\n");
3805 /* All terms except the last one */
3806 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3810 ir->opts.anneal_time[i][j],
3811 ir->opts.anneal_temp[i][j]);
3814 /* Finally the last one */
3815 j = ir->opts.anneal_npoints[i] - 1;
3816 if (ir->opts.annealing[i] == SimulatedAnnealing::Single)
3820 ir->opts.anneal_time[i][j],
3821 ir->opts.anneal_temp[i][j]);
3827 ir->opts.anneal_time[i][j],
3828 ir->opts.anneal_temp[i][j]);
3829 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3832 "There is a temperature jump when your annealing "
3844 for (int i = 1; i < ir->pull->ngroup; i++)
3846 const int gid = search_string(
3847 inputrecStrings->pullGroupNames[i].c_str(), defaultIndexGroups->nr, gnames);
3848 GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
3849 atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
3852 process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3854 checkPullCoords(ir->pull->group, ir->pull->coord);
3859 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3862 if (ir->eSwapCoords != SwapType::No)
3864 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3867 /* Make indices for IMD session */
3870 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3873 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3874 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3875 mdModulesNotifiers.preProcessingNotifier_.notify(defaultIndexGroupsAndNames);
3877 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3878 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3879 if (freezeDims.size() != DIM * freezeGroupNames.size())
3882 "Invalid Freezing input: %zu groups and %zu freeze values",
3883 freezeGroupNames.size(),
3886 do_numbering(natoms,
3891 SimulationAtomGroupType::Freeze,
3896 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3897 ir->opts.ngfrz = nr;
3898 snew(ir->opts.nFreeze, nr);
3899 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3901 for (j = 0; (j < DIM); j++, k++)
3903 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3904 if (!ir->opts.nFreeze[i][j])
3906 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3909 "Please use Y(ES) or N(O) for freezedim only "
3911 freezeDims[k].c_str());
3912 warning(wi, warn_buf);
3917 for (; (i < nr); i++)
3919 for (j = 0; (j < DIM); j++)
3921 ir->opts.nFreeze[i][j] = 0;
3925 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3926 do_numbering(natoms,
3931 SimulationAtomGroupType::EnergyOutput,
3936 add_wall_energrps(groups, ir->nwall, symtab);
3937 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3938 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3939 do_numbering(natoms,
3944 SimulationAtomGroupType::MassCenterVelocityRemoval,
3946 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART,
3950 if (ir->comm_mode != ComRemovalAlgorithm::No)
3952 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3955 /* Now we have filled the freeze struct, so we can calculate NRDF */
3956 calc_nrdf(mtop, ir, gnames);
3958 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3959 do_numbering(natoms,
3964 SimulationAtomGroupType::User1,
3969 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3970 do_numbering(natoms,
3975 SimulationAtomGroupType::User2,
3980 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3981 do_numbering(natoms,
3983 compressedXGroupNames,
3986 SimulationAtomGroupType::CompressedPositionOutput,
3991 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3992 do_numbering(natoms,
3994 orirefFitGroupNames,
3997 SimulationAtomGroupType::OrientationRestraintsFit,
4003 /* MiMiC QMMM input processing */
4004 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
4005 if (qmGroupNames.size() > 1)
4007 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
4009 /* group rest, if any, is always MM! */
4010 do_numbering(natoms,
4015 SimulationAtomGroupType::QuantumMechanics,
4020 ir->opts.ngQM = qmGroupNames.size();
4022 /* end of MiMiC QMMM input */
4026 for (auto group : gmx::keysOf(groups->groups))
4028 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
4029 for (const auto& entry : groups->groups[group])
4031 fprintf(stderr, " %s", *(groups->groupNames[entry]));
4033 fprintf(stderr, "\n");
4037 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
4038 snew(ir->opts.egp_flags, nr * nr);
4040 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
4041 if (bExcl && ir->cutoff_scheme == CutoffScheme::Verlet)
4043 warning_error(wi, "Energy group exclusions are currently not supported");
4045 if (bExcl && EEL_FULL(ir->coulombtype))
4047 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
4050 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
4051 if (bTable && !(ir->vdwtype == VanDerWaalsType::User)
4052 && !(ir->coulombtype == CoulombInteractionType::User)
4053 && !(ir->coulombtype == CoulombInteractionType::PmeUser)
4054 && !(ir->coulombtype == CoulombInteractionType::PmeUserSwitch))
4057 "Can only have energy group pair tables in combination with user tables for VdW "
4061 /* final check before going out of scope if simulated tempering variables
4062 * need to be set to default values.
4064 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
4066 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
4069 "the value for nstexpanded was not specified for "
4070 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
4071 "by default, but it is recommended to set it to an explicit value!",
4072 ir->expandedvals->nstexpanded));
4074 for (i = 0; (i < defaultIndexGroups->nr); i++)
4079 done_blocka(defaultIndexGroups);
4080 sfree(defaultIndexGroups);
4084 static void check_disre(const gmx_mtop_t& mtop)
4086 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
4088 const gmx_ffparams_t& ffparams = mtop.ffparams;
4091 for (int i = 0; i < ffparams.numTypes(); i++)
4093 int ftype = ffparams.functype[i];
4094 if (ftype == F_DISRES)
4096 int label = ffparams.iparams[i].disres.label;
4097 if (label == old_label)
4099 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
4108 "Found %d double distance restraint indices,\n"
4109 "probably the parameters for multiple pairs in one restraint "
4110 "are not identical\n",
4116 //! Returns whether dimensions have an absolute reference due to walls, pbc or freezing
4117 static BasicVector<bool> haveAbsoluteReference(const t_inputrec& ir)
4119 BasicVector<bool> absRef = { false, false, false };
4121 /* Check the degrees of freedom of the COM (not taking COMM removal into account) */
4122 for (int d = 0; d < DIM; d++)
4124 absRef[d] = (d >= ndof_com(&ir));
4126 /* Check for freeze groups */
4127 for (int g = 0; g < ir.opts.ngfrz; g++)
4129 for (int d = 0; d < DIM; d++)
4131 if (ir.opts.nFreeze[g][d] != 0)
4141 //! Returns whether position restraints are used for dimensions
4142 static BasicVector<bool> havePositionRestraints(const gmx_mtop_t& sys)
4144 BasicVector<bool> havePosres = { false, false, false };
4146 for (const auto ilists : IListRange(sys))
4148 const auto& posResList = ilists.list()[F_POSRES];
4149 const auto& fbPosResList = ilists.list()[F_FBPOSRES];
4150 if (ilists.nmol() > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
4152 for (int i = 0; i < posResList.size(); i += 2)
4154 const t_iparams& pr = sys.ffparams.iparams[posResList.iatoms[i]];
4155 for (int d = 0; d < DIM; d++)
4157 if (pr.posres.fcA[d] != 0)
4159 havePosres[d] = true;
4163 for (int i = 0; i < fbPosResList.size(); i += 2)
4165 /* Check for flat-bottom posres */
4166 const t_iparams& pr = sys.ffparams.iparams[fbPosResList.iatoms[i]];
4167 if (pr.fbposres.k != 0)
4169 switch (pr.fbposres.geom)
4171 case efbposresSPHERE: havePosres = { true, true, true }; break;
4172 case efbposresCYLINDERX: havePosres[YY] = havePosres[ZZ] = true; break;
4173 case efbposresCYLINDERY: havePosres[XX] = havePosres[ZZ] = true; break;
4174 case efbposresCYLINDER:
4175 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4176 case efbposresCYLINDERZ: havePosres[XX] = havePosres[YY] = true; break;
4177 case efbposresX: /* d=XX */
4178 case efbposresY: /* d=YY */
4179 case efbposresZ: /* d=ZZ */
4180 havePosres[pr.fbposres.geom - efbposresX] = true;
4184 "Invalid geometry for flat-bottom position restraint.\n"
4185 "Expected nr between 1 and %d. Found %d\n",
4197 static void check_combination_rule_differences(const gmx_mtop_t& mtop,
4199 bool* bC6ParametersWorkWithGeometricRules,
4200 bool* bC6ParametersWorkWithLBRules,
4201 bool* bLBRulesPossible)
4203 int ntypes, tpi, tpj;
4206 double c6i, c6j, c12i, c12j;
4207 double c6, c6_geometric, c6_LB;
4208 double sigmai, sigmaj, epsi, epsj;
4209 bool bCanDoLBRules, bCanDoGeometricRules;
4212 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4213 * force-field floating point parameters.
4216 ptr = getenv("GMX_LJCOMB_TOL");
4220 double gmx_unused canary;
4222 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4225 FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4230 *bC6ParametersWorkWithLBRules = TRUE;
4231 *bC6ParametersWorkWithGeometricRules = TRUE;
4232 bCanDoLBRules = TRUE;
4233 ntypes = mtop.ffparams.atnr;
4234 snew(typecount, ntypes);
4235 gmx_mtop_count_atomtypes(mtop, state, typecount);
4236 *bLBRulesPossible = TRUE;
4237 for (tpi = 0; tpi < ntypes; ++tpi)
4239 c6i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4240 c12i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4241 for (tpj = tpi; tpj < ntypes; ++tpj)
4243 c6j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4244 c12j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4245 c6 = mtop.ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4246 c6_geometric = std::sqrt(c6i * c6j);
4247 if (!gmx_numzero(c6_geometric))
4249 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4251 sigmai = gmx::sixthroot(c12i / c6i);
4252 sigmaj = gmx::sixthroot(c12j / c6j);
4253 epsi = c6i * c6i / (4.0 * c12i);
4254 epsj = c6j * c6j / (4.0 * c12j);
4255 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4259 *bLBRulesPossible = FALSE;
4260 c6_LB = c6_geometric;
4262 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4267 *bC6ParametersWorkWithLBRules = FALSE;
4270 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4272 if (!bCanDoGeometricRules)
4274 *bC6ParametersWorkWithGeometricRules = FALSE;
4281 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t& mtop, warninp_t wi)
4283 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4285 check_combination_rule_differences(
4286 mtop, 0, &bC6ParametersWorkWithGeometricRules, &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4287 if (ir->ljpme_combination_rule == LongRangeVdW::LB)
4289 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4292 "You are using arithmetic-geometric combination rules "
4293 "in LJ-PME, but your non-bonded C6 parameters do not "
4294 "follow these rules.");
4299 if (!bC6ParametersWorkWithGeometricRules)
4301 if (ir->eDispCorr != DispersionCorrectionType::No)
4304 "You are using geometric combination rules in "
4305 "LJ-PME, but your non-bonded C6 parameters do "
4306 "not follow these rules. "
4307 "This will introduce very small errors in the forces and energies in "
4308 "your simulations. Dispersion correction will correct total energy "
4309 "and/or pressure for isotropic systems, but not forces or surface "
4315 "You are using geometric combination rules in "
4316 "LJ-PME, but your non-bonded C6 parameters do "
4317 "not follow these rules. "
4318 "This will introduce very small errors in the forces and energies in "
4319 "your simulations. If your system is homogeneous, consider using "
4320 "dispersion correction "
4321 "for the total energy and pressure.");
4327 static bool allTrue(const BasicVector<bool>& boolVector)
4329 return boolVector[0] && boolVector[1] && boolVector[2];
4332 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4334 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4335 GMX_ASSERT(gmx::checkMtsRequirements(*ir).empty(), "All MTS requirements should be met here");
4337 char err_buf[STRLEN];
4340 gmx_mtop_atomloop_block_t aloopb;
4341 char warn_buf[STRLEN];
4343 set_warning_line(wi, mdparin, -1);
4345 if (allTrue(haveAbsoluteReference(*ir)) && allTrue(havePositionRestraints(*sys)))
4348 "Removing center of mass motion in the presence of position restraints might "
4349 "cause artifacts. When you are using position restraints to equilibrate a "
4350 "macro-molecule, the artifacts are usually negligible.");
4353 if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0 && ir->nstlist > 1
4354 && ((EI_MD(ir->eI) || EI_SD(ir->eI))
4355 && (ir->etc == TemperatureCoupling::VRescale || ir->etc == TemperatureCoupling::Berendsen)))
4357 /* Check if a too small Verlet buffer might potentially
4358 * cause more drift than the thermostat can couple off.
4360 /* Temperature error fraction for warning and suggestion */
4361 const real T_error_warn = 0.002;
4362 const real T_error_suggest = 0.001;
4363 /* For safety: 2 DOF per atom (typical with constraints) */
4364 const real nrdf_at = 2;
4365 real T, tau, max_T_error;
4370 for (i = 0; i < ir->opts.ngtc; i++)
4372 T = std::max(T, ir->opts.ref_t[i]);
4373 tau = std::max(tau, ir->opts.tau_t[i]);
4377 /* This is a worst case estimate of the temperature error,
4378 * assuming perfect buffer estimation and no cancelation
4379 * of errors. The factor 0.5 is because energy distributes
4380 * equally over Ekin and Epot.
4382 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * gmx::c_boltz * T);
4383 if (max_T_error > T_error_warn)
4386 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4387 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4388 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4389 "%.0e or decrease tau_t.",
4394 100 * T_error_suggest,
4395 ir->verletbuf_tol * T_error_suggest / max_T_error);
4396 warning(wi, warn_buf);
4401 if (ETC_ANDERSEN(ir->etc))
4405 for (i = 0; i < ir->opts.ngtc; i++)
4408 "all tau_t must currently be equal using Andersen temperature control, "
4409 "violated for group %d",
4411 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4413 "all tau_t must be positive using Andersen temperature control, "
4417 CHECK(ir->opts.tau_t[i] < 0);
4420 if (ir->etc == TemperatureCoupling::AndersenMassive && ir->comm_mode != ComRemovalAlgorithm::No)
4422 for (i = 0; i < ir->opts.ngtc; i++)
4424 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4426 "tau_t/delta_t for group %d for temperature control method %s must be a "
4427 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4428 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4431 enumValueToString(ir->etc),
4435 CHECK(nsteps % ir->nstcomm != 0);
4440 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != IntegrationAlgorithm::BD
4441 && ir->comm_mode == ComRemovalAlgorithm::No
4442 && !(allTrue(haveAbsoluteReference(*ir)) || allTrue(havePositionRestraints(*sys)) || ir->nsteps <= 10)
4443 && !ETC_ANDERSEN(ir->etc))
4446 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4447 "rounding errors can lead to build up of kinetic energy of the center of mass");
4450 if (ir->epc == PressureCoupling::ParrinelloRahman && ir->etc == TemperatureCoupling::NoseHoover)
4453 for (int g = 0; g < ir->opts.ngtc; g++)
4455 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4457 if (ir->tau_p < 1.9 * tau_t_max)
4459 std::string message = gmx::formatString(
4460 "With %s T-coupling and %s p-coupling, "
4461 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4462 enumValueToString(ir->etc),
4463 enumValueToString(ir->epc),
4468 warning(wi, message.c_str());
4472 /* Check for pressure coupling with absolute position restraints */
4473 if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == RefCoordScaling::No)
4475 const BasicVector<bool> havePosres = havePositionRestraints(*sys);
4477 for (m = 0; m < DIM; m++)
4479 if (havePosres[m] && norm2(ir->compress[m]) > 0)
4482 "You are using pressure coupling with absolute position restraints, "
4483 "this will give artifacts. Use the refcoord_scaling option.");
4491 aloopb = gmx_mtop_atomloop_block_init(*sys);
4493 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4495 if (atom->q != 0 || atom->qB != 0)
4503 if (EEL_FULL(ir->coulombtype))
4506 "You are using full electrostatics treatment %s for a system without charges.\n"
4507 "This costs a lot of performance for just processing zeros, consider using %s "
4509 enumValueToString(ir->coulombtype),
4510 enumValueToString(CoulombInteractionType::Cut));
4511 warning(wi, err_buf);
4516 if (ir->coulombtype == CoulombInteractionType::Cut && ir->rcoulomb > 0)
4519 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4520 "You might want to consider using %s electrostatics.\n",
4521 enumValueToString(CoulombInteractionType::Pme));
4522 warning_note(wi, err_buf);
4526 /* Check if combination rules used in LJ-PME are the same as in the force field */
4527 if (EVDW_PME(ir->vdwtype))
4529 check_combination_rules(ir, *sys, wi);
4532 /* Generalized reaction field */
4533 if (ir->coulombtype == CoulombInteractionType::GRFNotused)
4536 "Generalized reaction-field electrostatics is no longer supported. "
4537 "You can use normal reaction-field instead and compute the reaction-field "
4538 "constant by hand.");
4541 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->sc_alpha != 0
4542 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4544 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4552 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4554 if (ir->pull->coord[i].eGeom != PullGroupGeometry::Transformation
4555 && (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0))
4557 const auto absRef = haveAbsoluteReference(*ir);
4558 const auto havePosres = havePositionRestraints(*sys);
4559 for (m = 0; m < DIM; m++)
4561 if (ir->pull->coord[i].dim[m] && !(absRef[m] || havePosres[m]))
4564 "You are using an absolute reference for pulling, but the rest of "
4565 "the system does not have an absolute reference. This will lead to "
4574 for (i = 0; i < 3; i++)
4576 for (m = 0; m <= i; m++)
4578 if ((ir->epc != PressureCoupling::No && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4580 for (c = 0; c < ir->pull->ncoord; c++)
4582 if (ir->pull->coord[c].eGeom == PullGroupGeometry::DirectionPBC
4583 && ir->pull->coord[c].vec[m] != 0)
4586 "Can not have dynamic box while using pull geometry '%s' "
4588 enumValueToString(ir->pull->coord[c].eGeom),
4600 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4602 char warn_buf[STRLEN];
4605 ptr = check_box(ir->pbcType, box);
4608 warning_error(wi, ptr);
4611 if (bHasNormalConstraints && ir->eConstrAlg == ConstraintAlgorithm::Shake)
4613 if (ir->shake_tol <= 0.0)
4615 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4616 warning_error(wi, warn_buf);
4620 if ((ir->eConstrAlg == ConstraintAlgorithm::Lincs) && bHasNormalConstraints)
4622 /* If we have Lincs constraints: */
4623 if (ir->eI == IntegrationAlgorithm::MD && ir->etc == TemperatureCoupling::No
4624 && ir->eConstrAlg == ConstraintAlgorithm::Lincs && ir->nLincsIter == 1)
4627 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4628 warning_note(wi, warn_buf);
4631 if ((ir->eI == IntegrationAlgorithm::CG || ir->eI == IntegrationAlgorithm::LBFGS)
4632 && (ir->nProjOrder < 8))
4635 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4636 enumValueToString(ir->eI));
4637 warning_note(wi, warn_buf);
4639 if (ir->epc == PressureCoupling::Mttk)
4641 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4645 if (bHasAnyConstraints && ir->epc == PressureCoupling::Mttk)
4647 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4650 if (ir->LincsWarnAngle > 90.0)
4652 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4653 warning(wi, warn_buf);
4654 ir->LincsWarnAngle = 90.0;
4657 if (ir->pbcType != PbcType::No)
4659 if (ir->nstlist == 0)
4662 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4663 "atoms might cause the simulation to crash.");
4665 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4668 "ERROR: The cut-off length is longer than half the shortest box vector or "
4669 "longer than the smallest box diagonal element. Increase the box size or "
4670 "decrease rlist.\n");
4671 warning_error(wi, warn_buf);