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;
142 //! How to treat coverage of the whole system for a set of atom groupsx
143 enum class GroupCoverage
145 All, //!< All particles have to be a member of a group
146 AllGenerateRest, //<! A rest group with name is generated for particles not part of any group
147 Partial, //<! As \p AllGenerateRest, but no name for the rest group is generated
148 OneGroup //<! Merge all selected groups into one group, make a rest group for the remaining particles
151 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
152 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
153 "h-angles", "all-angles", nullptr };
155 // NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
156 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
158 static void getSimTemps(int ntemps, t_simtemp* simtemp, gmx::ArrayRef<double> temperature_lambdas)
163 for (i = 0; i < ntemps; i++)
165 /* simple linear scaling -- allows more control */
166 if (simtemp->eSimTempScale == SimulatedTempering::Linear)
168 simtemp->temperatures[i] =
170 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
172 else if (simtemp->eSimTempScale
173 == SimulatedTempering::Geometric) /* should give roughly equal acceptance for constant heat capacity . . . */
175 simtemp->temperatures[i] = simtemp->simtemp_low
176 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
177 static_cast<real>((1.0 * i) / (ntemps - 1)));
179 else if (simtemp->eSimTempScale == SimulatedTempering::Exponential)
181 simtemp->temperatures[i] = simtemp->simtemp_low
182 + (simtemp->simtemp_high - simtemp->simtemp_low)
183 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
188 sprintf(errorstr, "eSimTempScale=%s not defined", enumValueToString(simtemp->eSimTempScale));
189 gmx_fatal(FARGS, "%s", errorstr);
195 static void _low_check(bool b, const char* s, warninp_t wi)
199 warning_error(wi, s);
203 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
207 if (*p > 0 && *p % nst != 0)
209 /* Round up to the next multiple of nst */
210 *p = ((*p) / nst + 1) * nst;
211 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
216 //! Convert legacy mdp entries to modern ones.
217 static void process_interaction_modifier(InteractionModifiers* eintmod)
219 if (*eintmod == InteractionModifiers::PotShiftVerletUnsupported)
221 *eintmod = InteractionModifiers::PotShift;
225 void check_ir(const char* mdparin,
226 const gmx::MDModulesNotifiers& mdModulesNotifiers,
230 /* Check internal consistency.
231 * NOTE: index groups are not set here yet, don't check things
232 * like temperature coupling group options here, but in triple_check
235 /* Strange macro: first one fills the err_buf, and then one can check
236 * the condition, which will print the message and increase the error
239 #define CHECK(b) _low_check(b, err_buf, wi)
240 char err_buf[256], warn_buf[STRLEN];
243 t_lambda* fep = ir->fepvals.get();
244 t_expanded* expand = ir->expandedvals.get();
246 set_warning_line(wi, mdparin, -1);
248 /* We cannot check MTS requirements with an invalid MTS setup
249 * and we will already have generated errors with an invalid MTS setup.
251 if (gmx::haveValidMtsSetup(*ir))
253 std::vector<std::string> errorMessages = gmx::checkMtsRequirements(*ir);
255 for (const auto& errorMessage : errorMessages)
257 warning_error(wi, errorMessage.c_str());
261 if (ir->coulombtype == CoulombInteractionType::RFNecUnsupported)
263 std::string message =
264 gmx::formatString("%s electrostatics is no longer supported",
265 enumValueToString(CoulombInteractionType::RFNecUnsupported));
266 warning_error(wi, message);
269 /* BASIC CUT-OFF STUFF */
270 if (ir->rcoulomb < 0)
272 warning_error(wi, "rcoulomb should be >= 0");
276 warning_error(wi, "rvdw should be >= 0");
278 if (ir->rlist < 0 && !(ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0))
280 warning_error(wi, "rlist should be >= 0");
283 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
284 "neighbour-list update scheme for efficient buffering for improved energy "
285 "conservation, please use the Verlet cut-off scheme instead.)");
286 CHECK(ir->nstlist < 0);
288 process_interaction_modifier(&ir->coulomb_modifier);
289 process_interaction_modifier(&ir->vdw_modifier);
291 if (ir->cutoff_scheme == CutoffScheme::Group)
294 "The group cutoff scheme has been removed since GROMACS 2020. "
295 "Please use the Verlet cutoff scheme.");
297 if (ir->cutoff_scheme == CutoffScheme::Verlet)
301 /* Normal Verlet type neighbor-list, currently only limited feature support */
302 if (inputrec2nboundeddim(ir) < 3)
304 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
307 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
308 if (ir->rcoulomb != ir->rvdw)
310 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
311 // for PME load balancing, we can support this exception.
312 bool bUsesPmeTwinRangeKernel =
313 (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == VanDerWaalsType::Cut
314 && ir->rcoulomb > ir->rvdw);
315 if (!bUsesPmeTwinRangeKernel)
318 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
319 "rcoulomb>rvdw with PME electrostatics)");
323 if (ir->vdwtype == VanDerWaalsType::Shift || ir->vdwtype == VanDerWaalsType::Switch)
325 if (ir->vdw_modifier == InteractionModifiers::None
326 || ir->vdw_modifier == InteractionModifiers::PotShift)
329 (ir->vdwtype == VanDerWaalsType::Shift ? InteractionModifiers::ForceSwitch
330 : InteractionModifiers::PotSwitch);
333 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
335 enumValueToString(ir->vdwtype),
336 enumValueToString(VanDerWaalsType::Cut),
337 enumValueToString(ir->vdw_modifier));
338 warning_note(wi, warn_buf);
340 ir->vdwtype = VanDerWaalsType::Cut;
345 "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
346 enumValueToString(ir->vdwtype),
347 enumValueToString(ir->vdw_modifier));
348 warning_error(wi, warn_buf);
352 if (!(ir->vdwtype == VanDerWaalsType::Cut || ir->vdwtype == VanDerWaalsType::Pme))
355 "With Verlet lists only cut-off and PME LJ interactions are supported");
357 if (!(ir->coulombtype == CoulombInteractionType::Cut || EEL_RF(ir->coulombtype)
358 || EEL_PME(ir->coulombtype) || ir->coulombtype == CoulombInteractionType::Ewald))
361 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
362 "electrostatics are supported");
364 if (!(ir->coulomb_modifier == InteractionModifiers::None
365 || ir->coulomb_modifier == InteractionModifiers::PotShift))
367 sprintf(warn_buf, "coulomb_modifier=%s is not supported", enumValueToString(ir->coulomb_modifier));
368 warning_error(wi, warn_buf);
371 if (EEL_USER(ir->coulombtype))
374 "Coulomb type %s is not supported with the verlet scheme",
375 enumValueToString(ir->coulombtype));
376 warning_error(wi, warn_buf);
379 if (ir->nstlist <= 0)
381 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
384 if (ir->nstlist < 10)
387 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
388 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
392 rc_max = std::max(ir->rvdw, ir->rcoulomb);
396 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
397 ir->verletbuf_tol = 0;
400 else if (ir->verletbuf_tol <= 0)
402 if (ir->verletbuf_tol == 0)
404 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
407 if (ir->rlist < rc_max)
410 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
413 if (ir->rlist == rc_max && ir->nstlist > 1)
417 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
418 "buffer. The cluster pair list does have a buffering effect, but choosing "
419 "a larger rlist might be necessary for good energy conservation.");
424 if (ir->rlist > rc_max)
427 "You have set rlist larger than the interaction cut-off, but you also "
428 "have verlet-buffer-tolerance > 0. Will set rlist using "
429 "verlet-buffer-tolerance.");
432 if (ir->nstlist == 1)
434 /* No buffer required */
439 if (EI_DYNAMICS(ir->eI))
441 if (inputrec2nboundeddim(ir) < 3)
444 "The box volume is required for calculating rlist from the "
445 "energy drift with verlet-buffer-tolerance > 0. You are "
446 "using at least one unbounded dimension, so no volume can be "
447 "computed. Either use a finite box, or set rlist yourself "
448 "together with verlet-buffer-tolerance = -1.");
450 /* Set rlist temporarily so we can continue processing */
455 /* Set the buffer to 5% of the cut-off */
456 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
462 /* GENERAL INTEGRATOR STUFF */
465 if (ir->etc != TemperatureCoupling::No)
467 if (EI_RANDOM(ir->eI))
470 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
471 "implicitly. See the documentation for more information on which "
472 "parameters affect temperature for %s.",
473 enumValueToString(ir->etc),
474 enumValueToString(ir->eI),
475 enumValueToString(ir->eI));
480 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
482 enumValueToString(ir->etc),
483 enumValueToString(ir->eI));
485 warning_note(wi, warn_buf);
487 ir->etc = TemperatureCoupling::No;
489 if (ir->eI == IntegrationAlgorithm::VVAK)
492 "Integrator method %s is implemented primarily for validation purposes; for "
493 "molecular dynamics, you should probably be using %s or %s",
494 enumValueToString(IntegrationAlgorithm::VVAK),
495 enumValueToString(IntegrationAlgorithm::MD),
496 enumValueToString(IntegrationAlgorithm::VV));
497 warning_note(wi, warn_buf);
499 if (!EI_DYNAMICS(ir->eI))
501 if (ir->epc != PressureCoupling::No)
504 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
505 enumValueToString(ir->epc),
506 enumValueToString(ir->eI));
507 warning_note(wi, warn_buf);
509 ir->epc = PressureCoupling::No;
511 if (EI_DYNAMICS(ir->eI))
513 if (ir->nstcalcenergy < 0)
515 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
516 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
518 /* nstcalcenergy larger than nstener does not make sense.
519 * We ideally want nstcalcenergy=nstener.
523 ir->nstcalcenergy = std::gcd(ir->nstenergy, ir->nstlist);
527 ir->nstcalcenergy = ir->nstenergy;
531 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
532 || (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl > 0
533 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
536 const char* nsten = "nstenergy";
537 const char* nstdh = "nstdhdl";
538 const char* min_name = nsten;
539 int min_nst = ir->nstenergy;
541 /* find the smallest of ( nstenergy, nstdhdl ) */
542 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->nstdhdl > 0
543 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
545 min_nst = ir->fepvals->nstdhdl;
548 /* If the user sets nstenergy small, we should respect that */
549 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy, min_name, min_nst);
550 warning_note(wi, warn_buf);
551 ir->nstcalcenergy = min_nst;
554 if (ir->epc != PressureCoupling::No)
556 if (ir->nstpcouple < 0)
558 ir->nstpcouple = ir_optimal_nstpcouple(ir);
560 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
563 "With multiple time stepping, nstpcouple should be a mutiple of "
568 if (ir->nstcalcenergy > 0)
570 if (ir->efep != FreeEnergyPerturbationType::No)
572 /* nstdhdl should be a multiple of nstcalcenergy */
573 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
577 /* nstexpanded should be a multiple of nstcalcenergy */
578 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded", &ir->expandedvals->nstexpanded, wi);
580 /* for storing exact averages nstenergy should be
581 * a multiple of nstcalcenergy
583 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
586 // Inquire all MDModules, if their parameters match with the energy
587 // calculation frequency
588 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
589 mdModulesNotifiers.preProcessingNotifier_.notify(&energyCalculationFrequencyErrors);
591 // Emit all errors from the energy calculation frequency checks
592 for (const std::string& energyFrequencyErrorMessage :
593 energyCalculationFrequencyErrors.errorMessages())
595 warning_error(wi, energyFrequencyErrorMessage);
599 if (ir->nsteps == 0 && !ir->bContinuation)
602 "For a correct single-point energy evaluation with nsteps = 0, use "
603 "continuation = yes to avoid constraining the input coordinates.");
607 if ((EI_SD(ir->eI) || ir->eI == IntegrationAlgorithm::BD) && ir->bContinuation && ir->ld_seed != -1)
610 "You are doing a continuation with SD or BD, make sure that ld_seed is "
611 "different from the previous run (using ld_seed=-1 will ensure this)");
617 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
618 CHECK(ir->pbcType != PbcType::Xyz);
619 sprintf(err_buf, "with TPI nstlist should be larger than zero");
620 CHECK(ir->nstlist <= 0);
621 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
622 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
626 if ((opts->nshake > 0) && (opts->bMorse))
628 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
629 warning(wi, warn_buf);
632 if ((EI_SD(ir->eI) || ir->eI == IntegrationAlgorithm::BD) && ir->bContinuation && ir->ld_seed != -1)
635 "You are doing a continuation with SD or BD, make sure that ld_seed is "
636 "different from the previous run (using ld_seed=-1 will ensure this)");
638 /* verify simulated tempering options */
642 bool bAllTempZero = TRUE;
643 for (i = 0; i < fep->n_lambda; i++)
646 "Entry %d for %s must be between 0 and 1, instead is %g",
648 enumValueToString(FreeEnergyPerturbationCouplingType::Temperature),
649 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i]);
650 CHECK((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i] < 0)
651 || (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i]
653 if (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Temperature)][i] > 0)
655 bAllTempZero = FALSE;
658 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
659 CHECK(bAllTempZero == TRUE);
661 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
662 CHECK(ir->eI != IntegrationAlgorithm::VV);
664 /* check compatability of the temperature coupling with simulated tempering */
666 if (ir->etc == TemperatureCoupling::NoseHoover)
669 "Nose-Hoover based temperature control such as [%s] my not be "
670 "entirelyconsistent with simulated tempering",
671 enumValueToString(ir->etc));
672 warning_note(wi, warn_buf);
675 /* check that the temperatures make sense */
678 "Higher simulated tempering temperature (%g) must be >= than the simulated "
679 "tempering lower temperature (%g)",
680 ir->simtempvals->simtemp_high,
681 ir->simtempvals->simtemp_low);
682 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
685 "Higher simulated tempering temperature (%g) must be >= zero",
686 ir->simtempvals->simtemp_high);
687 CHECK(ir->simtempvals->simtemp_high <= 0);
690 "Lower simulated tempering temperature (%g) must be >= zero",
691 ir->simtempvals->simtemp_low);
692 CHECK(ir->simtempvals->simtemp_low <= 0);
695 /* verify free energy options */
697 if (ir->efep != FreeEnergyPerturbationType::No)
699 fep = ir->fepvals.get();
700 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
701 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
704 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
706 static_cast<int>(fep->sc_r_power));
707 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
710 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
713 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
716 "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
718 CHECK(fep->delta_lambda > 0 && (ir->efep == FreeEnergyPerturbationType::Expanded));
720 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
721 CHECK(!(EI_VV(ir->eI)) && (ir->efep == FreeEnergyPerturbationType::Expanded));
723 sprintf(err_buf, "Free-energy not implemented for Ewald");
724 CHECK(ir->coulombtype == CoulombInteractionType::Ewald);
726 /* check validty of lambda inputs */
727 if (fep->n_lambda == 0)
729 /* Clear output in case of no states:*/
730 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
731 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
736 "initial thermodynamic state %d does not exist, only goes to %d",
739 CHECK((fep->init_fep_state >= fep->n_lambda));
743 "Lambda state must be set, either with init-lambda-state or with init-lambda");
744 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
747 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
748 "init-lambda-state or with init-lambda, but not both",
750 fep->init_fep_state);
751 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
754 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
758 for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
760 if (fep->separate_dvdl[i])
765 if (n_lambda_terms > 1)
768 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
769 "use init-lambda to set lambda state (except for slow growth). Use "
770 "init-lambda-state instead.");
771 warning(wi, warn_buf);
774 if (n_lambda_terms < 2 && fep->n_lambda > 0)
777 "init-lambda is deprecated for setting lambda state (except for slow "
778 "growth). Use init-lambda-state instead.");
782 for (j = 0; j < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); j++)
784 for (i = 0; i < fep->n_lambda; i++)
786 auto enumValue = static_cast<FreeEnergyPerturbationCouplingType>(j);
788 "Entry %d for %s must be between 0 and 1, instead is %g",
790 enumValueToString(enumValue),
791 fep->all_lambda[j][i]);
792 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
796 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
798 for (i = 0; i < fep->n_lambda; i++)
801 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
802 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
803 "crashes, and is not supported.",
805 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i],
806 fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i]);
807 CHECK((fep->sc_alpha > 0)
808 && (((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i] > 0.0)
809 && (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Coul)][i] < 1.0))
810 && ((fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i] > 0.0)
811 && (fep->all_lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Vdw)][i]
816 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
818 real sigma, lambda, r_sc;
821 /* Maximum estimate for A and B charges equal with lambda power 1 */
823 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
824 1.0 / fep->sc_r_power);
826 "With PME there is a minor soft core effect present at the cut-off, "
827 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
828 "energy conservation, but usually other effects dominate. With a common sigma "
829 "value of %g nm the fraction of the particle-particle potential at the cut-off "
830 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
836 warning_note(wi, warn_buf);
839 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
840 be treated differently, but that's the next step */
842 for (i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
844 auto enumValue = static_cast<FreeEnergyPerturbationCouplingType>(i);
845 for (j = 0; j < fep->n_lambda; j++)
847 sprintf(err_buf, "%s[%d] must be between 0 and 1", enumValueToString(enumValue), j);
848 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
852 if (fep->softcoreFunction == SoftcoreType::Gapsys)
854 if (fep->scGapsysScaleLinpointQ < 0.0)
857 "sc_scale_linpoint_Q_gapsys is equal %g but must be >= 0",
858 fep->scGapsysScaleLinpointQ);
859 warning_note(wi, warn_buf);
862 if ((fep->scGapsysScaleLinpointLJ < 0.0) || (fep->scGapsysScaleLinpointLJ >= 1.0))
865 "sc_scale_linpoint_LJ_gapsys is equal %g but must be in [0,1) when used "
867 "sc_function=gapsys.",
868 fep->scGapsysScaleLinpointLJ);
869 warning_note(wi, warn_buf);
874 if ((ir->bSimTemp) || (ir->efep == FreeEnergyPerturbationType::Expanded))
876 fep = ir->fepvals.get();
878 /* checking equilibration of weights inputs for validity */
881 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
883 expand->equil_n_at_lam,
884 enumValueToString(LambdaWeightWillReachEquilibrium::NumAtLambda));
885 CHECK((expand->equil_n_at_lam > 0)
886 && (expand->elmceq != LambdaWeightWillReachEquilibrium::NumAtLambda));
889 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
891 expand->equil_samples,
892 enumValueToString(LambdaWeightWillReachEquilibrium::Samples));
893 CHECK((expand->equil_samples > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Samples));
896 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
898 enumValueToString(LambdaWeightWillReachEquilibrium::Steps));
899 CHECK((expand->equil_steps > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Steps));
902 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
903 expand->equil_samples,
904 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta));
905 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::WLDelta));
908 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
910 enumValueToString(LambdaWeightWillReachEquilibrium::Ratio));
911 CHECK((expand->equil_ratio > 0) && (expand->elmceq != LambdaWeightWillReachEquilibrium::Ratio));
914 "weight-equil-number-all-lambda (%d) must be a positive integer if "
915 "lmc-weights-equil=%s",
916 expand->equil_n_at_lam,
917 enumValueToString(LambdaWeightWillReachEquilibrium::NumAtLambda));
918 CHECK((expand->equil_n_at_lam <= 0)
919 && (expand->elmceq == LambdaWeightWillReachEquilibrium::NumAtLambda));
922 "weight-equil-number-samples (%d) must be a positive integer if "
923 "lmc-weights-equil=%s",
924 expand->equil_samples,
925 enumValueToString(LambdaWeightWillReachEquilibrium::Samples));
926 CHECK((expand->equil_samples <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Samples));
929 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
931 enumValueToString(LambdaWeightWillReachEquilibrium::Steps));
932 CHECK((expand->equil_steps <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Steps));
935 "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
936 expand->equil_wl_delta,
937 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta));
938 CHECK((expand->equil_wl_delta <= 0)
939 && (expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta));
942 "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
944 enumValueToString(LambdaWeightWillReachEquilibrium::Ratio));
945 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == LambdaWeightWillReachEquilibrium::Ratio));
948 "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
949 enumValueToString(LambdaWeightWillReachEquilibrium::WLDelta),
950 enumValueToString(LambdaWeightCalculation::WL),
951 enumValueToString(LambdaWeightCalculation::WWL));
952 CHECK((expand->elmceq == LambdaWeightWillReachEquilibrium::WLDelta) && (!EWL(expand->elamstats)));
954 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
955 CHECK((expand->lmc_repeats <= 0));
956 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
957 CHECK((expand->minvarmin <= 0));
958 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
959 CHECK((expand->c_range < 0));
961 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
964 expand->lmc_forced_nstart);
965 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
966 && (expand->elmcmove != LambdaMoveCalculation::No));
967 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
968 CHECK((expand->lmc_forced_nstart < 0));
970 "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
971 fep->init_fep_state);
972 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
974 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
975 CHECK((expand->init_wl_delta < 0));
976 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
977 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
978 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
979 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
981 /* if there is no temperature control, we need to specify an MC temperature */
982 if (!integratorHasReferenceTemperature(ir)
983 && (expand->elmcmove != LambdaMoveCalculation::No) && (expand->mc_temp <= 0.0))
986 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
987 "to a positive number");
988 warning_error(wi, err_buf);
990 if (expand->nstTij > 0)
992 sprintf(err_buf, "nstlog must be non-zero");
993 CHECK(ir->nstlog == 0);
994 // Avoid modulus by zero in the case that already triggered an error exit.
998 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
1001 CHECK((expand->nstTij % ir->nstlog) != 0);
1007 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
1008 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
1011 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
1013 if (ir->pbcType == PbcType::No)
1015 if (ir->epc != PressureCoupling::No)
1017 warning(wi, "Turning off pressure coupling for vacuum system");
1018 ir->epc = PressureCoupling::No;
1024 "Can not have pressure coupling with pbc=%s",
1025 c_pbcTypeNames[ir->pbcType].c_str());
1026 CHECK(ir->epc != PressureCoupling::No);
1028 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
1029 CHECK(EEL_FULL(ir->coulombtype));
1032 "Can not have dispersion correction with pbc=%s",
1033 c_pbcTypeNames[ir->pbcType].c_str());
1034 CHECK(ir->eDispCorr != DispersionCorrectionType::No);
1037 if (ir->rlist == 0.0)
1040 "can only have neighborlist cut-off zero (=infinite)\n"
1041 "with coulombtype = %s or coulombtype = %s\n"
1042 "without periodic boundary conditions (pbc = %s) and\n"
1043 "rcoulomb and rvdw set to zero",
1044 enumValueToString(CoulombInteractionType::Cut),
1045 enumValueToString(CoulombInteractionType::User),
1046 c_pbcTypeNames[PbcType::No].c_str());
1047 CHECK(((ir->coulombtype != CoulombInteractionType::Cut)
1048 && (ir->coulombtype != CoulombInteractionType::User))
1049 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1051 if (ir->nstlist > 0)
1054 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1055 "nstype=simple and only one MPI rank");
1060 if (ir->nstcomm == 0)
1062 // TODO Change this behaviour. There should be exactly one way
1063 // to turn off an algorithm.
1064 ir->comm_mode = ComRemovalAlgorithm::No;
1066 if (ir->comm_mode != ComRemovalAlgorithm::No)
1068 if (ir->nstcomm < 0)
1070 // TODO Such input was once valid. Now that we've been
1071 // helpful for a few years, we should reject such input,
1072 // lest we have to support every historical decision
1075 "If you want to remove the rotation around the center of mass, you should set "
1076 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1077 "its absolute value");
1078 ir->nstcomm = abs(ir->nstcomm);
1081 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy
1082 && ir->comm_mode != ComRemovalAlgorithm::LinearAccelerationCorrection)
1085 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, consider "
1086 "setting nstcomm equal to nstcalcenergy for less overhead");
1089 if (ir->comm_mode == ComRemovalAlgorithm::Angular)
1092 "Can not remove the rotation around the center of mass with periodic "
1094 CHECK(ir->bPeriodicMols);
1095 if (ir->pbcType != PbcType::No)
1098 "Removing the rotation around the center of mass in a periodic system, "
1099 "this can lead to artifacts. Only use this on a single (cluster of) "
1100 "molecules. This cluster should not cross periodic boundaries.");
1105 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No
1106 && ir->comm_mode != ComRemovalAlgorithm::Angular)
1109 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1110 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1112 enumValueToString(IntegrationAlgorithm::SD1));
1113 warning_note(wi, warn_buf);
1116 /* TEMPERATURE COUPLING */
1117 if (ir->etc == TemperatureCoupling::Yes)
1119 ir->etc = TemperatureCoupling::Berendsen;
1121 "Old option for temperature coupling given: "
1122 "changing \"yes\" to \"Berendsen\"\n");
1125 if ((ir->etc == TemperatureCoupling::NoseHoover) || (ir->epc == PressureCoupling::Mttk))
1127 if (ir->opts.nhchainlength < 1)
1130 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1132 ir->opts.nhchainlength);
1133 ir->opts.nhchainlength = 1;
1134 warning(wi, warn_buf);
1137 if (ir->etc == TemperatureCoupling::NoseHoover && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1141 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1142 ir->opts.nhchainlength = 1;
1147 ir->opts.nhchainlength = 0;
1150 if (ir->eI == IntegrationAlgorithm::VVAK)
1153 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1155 enumValueToString(IntegrationAlgorithm::VVAK));
1156 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1159 if (ETC_ANDERSEN(ir->etc))
1162 "%s temperature control not supported for integrator %s.",
1163 enumValueToString(ir->etc),
1164 enumValueToString(ir->eI));
1165 CHECK(!(EI_VV(ir->eI)));
1167 if (ir->nstcomm > 0 && (ir->etc == TemperatureCoupling::Andersen))
1170 "Center of mass removal not necessary for %s. All velocities of coupled "
1171 "groups are rerandomized periodically, so flying ice cube errors will not "
1173 enumValueToString(ir->etc));
1174 warning_note(wi, warn_buf);
1178 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1179 "randomized every time step",
1181 enumValueToString(ir->etc));
1182 CHECK(ir->nstcomm > 1 && (ir->etc == TemperatureCoupling::Andersen));
1185 if (ir->etc == TemperatureCoupling::Berendsen)
1188 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1189 "might want to consider using the %s thermostat.",
1190 enumValueToString(ir->etc),
1191 enumValueToString(TemperatureCoupling::VRescale));
1192 warning_note(wi, warn_buf);
1195 if ((ir->etc == TemperatureCoupling::NoseHoover || ETC_ANDERSEN(ir->etc))
1196 && ir->epc == PressureCoupling::Berendsen)
1199 "Using Berendsen pressure coupling invalidates the "
1200 "true ensemble for the thermostat");
1201 warning(wi, warn_buf);
1204 /* PRESSURE COUPLING */
1205 if (ir->epc == PressureCoupling::Isotropic)
1207 ir->epc = PressureCoupling::Berendsen;
1209 "Old option for pressure coupling given: "
1210 "changing \"Isotropic\" to \"Berendsen\"\n");
1213 if (ir->epc != PressureCoupling::No)
1215 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1217 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1218 CHECK(ir->tau_p <= 0);
1220 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1223 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1224 "times larger than nstpcouple*dt (%g)",
1225 enumValueToString(ir->epc),
1227 pcouple_min_integration_steps(ir->epc),
1229 warning(wi, warn_buf);
1233 "compressibility must be > 0 when using pressure"
1235 enumValueToString(ir->epc));
1236 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1237 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1238 && ir->compress[ZZ][YY] <= 0));
1240 if (PressureCoupling::ParrinelloRahman == ir->epc && opts->bGenVel)
1243 "You are generating velocities so I am assuming you "
1244 "are equilibrating a system. You are using "
1245 "%s pressure coupling, but this can be "
1246 "unstable for equilibration. If your system crashes, try "
1247 "equilibrating first with Berendsen pressure coupling. If "
1248 "you are not equilibrating the system, you can probably "
1249 "ignore this warning.",
1250 enumValueToString(ir->epc));
1251 warning(wi, warn_buf);
1257 if (ir->epc == PressureCoupling::Mttk)
1259 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1263 /* ELECTROSTATICS */
1264 /* More checks are in triple check (grompp.c) */
1266 if (ir->coulombtype == CoulombInteractionType::Switch)
1269 "coulombtype = %s is only for testing purposes and can lead to serious "
1270 "artifacts, advice: use coulombtype = %s",
1271 enumValueToString(ir->coulombtype),
1272 enumValueToString(CoulombInteractionType::RFZero));
1273 warning(wi, warn_buf);
1276 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1279 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1280 "format and exchanging epsilon-r and epsilon-rf",
1282 warning(wi, warn_buf);
1283 ir->epsilon_rf = ir->epsilon_r;
1284 ir->epsilon_r = 1.0;
1287 if (ir->epsilon_r == 0)
1290 "It is pointless to use long-range electrostatics with infinite relative "
1292 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1294 CHECK(EEL_FULL(ir->coulombtype));
1297 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1299 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1300 CHECK(ir->epsilon_r < 0);
1303 if (EEL_RF(ir->coulombtype))
1305 /* reaction field (at the cut-off) */
1307 if (ir->coulombtype == CoulombInteractionType::RFZero && ir->epsilon_rf != 0)
1310 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1311 enumValueToString(ir->coulombtype));
1312 warning(wi, warn_buf);
1313 ir->epsilon_rf = 0.0;
1316 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1317 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1318 if (ir->epsilon_rf == ir->epsilon_r)
1321 "Using epsilon-rf = epsilon-r with %s does not make sense",
1322 enumValueToString(ir->coulombtype));
1323 warning(wi, warn_buf);
1326 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1327 * means the interaction is zero outside rcoulomb, but it helps to
1328 * provide accurate energy conservation.
1330 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1332 if (ir_coulomb_switched(ir))
1335 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1336 "potential modifier options!",
1337 enumValueToString(ir->coulombtype));
1338 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1342 if (ir->coulombtype == CoulombInteractionType::Switch || ir->coulombtype == CoulombInteractionType::Shift)
1345 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1346 "secondary coulomb-modifier.");
1347 CHECK(ir->coulomb_modifier != InteractionModifiers::None);
1349 if (ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift)
1352 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1353 "secondary vdw-modifier.");
1354 CHECK(ir->vdw_modifier != InteractionModifiers::None);
1357 if (ir->coulombtype == CoulombInteractionType::Switch || ir->coulombtype == CoulombInteractionType::Shift
1358 || ir->vdwtype == VanDerWaalsType::Switch || ir->vdwtype == VanDerWaalsType::Shift)
1361 "The switch/shift interaction settings are just for compatibility; you will get "
1363 "performance from applying potential modifiers to your interactions!\n");
1364 warning_note(wi, warn_buf);
1367 if (ir->coulombtype == CoulombInteractionType::PmeSwitch
1368 || ir->coulomb_modifier == InteractionModifiers::PotSwitch)
1370 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1372 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1374 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1375 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1376 "will be good regardless, since ewald_rtol = %g.",
1378 ir->rcoulomb_switch,
1381 warning(wi, warn_buf);
1385 if (ir->vdwtype == VanDerWaalsType::Switch || ir->vdw_modifier == InteractionModifiers::PotSwitch)
1387 if (ir->rvdw_switch == 0)
1390 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1391 "potential. This suggests it was not set in the mdp, which can lead to large "
1392 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1393 "switching range.");
1394 warning(wi, warn_buf);
1398 if (EEL_FULL(ir->coulombtype))
1400 if (ir->coulombtype == CoulombInteractionType::PmeSwitch
1401 || ir->coulombtype == CoulombInteractionType::PmeUser
1402 || ir->coulombtype == CoulombInteractionType::PmeUserSwitch)
1405 "With coulombtype = %s, rcoulomb must be <= rlist",
1406 enumValueToString(ir->coulombtype));
1407 CHECK(ir->rcoulomb > ir->rlist);
1411 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1413 // TODO: Move these checks into the ewald module with the options class
1415 int orderMax = (ir->coulombtype == CoulombInteractionType::P3mAD ? 8 : 12);
1417 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1420 "With coulombtype = %s, you should have %d <= pme-order <= %d",
1421 enumValueToString(ir->coulombtype),
1424 warning_error(wi, warn_buf);
1428 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1430 if (ir->ewald_geometry == EwaldGeometry::ThreeD)
1433 "With pbc=%s you should use ewald-geometry=%s",
1434 c_pbcTypeNames[ir->pbcType].c_str(),
1435 enumValueToString(EwaldGeometry::ThreeDC));
1436 warning(wi, warn_buf);
1438 /* This check avoids extra pbc coding for exclusion corrections */
1439 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1440 CHECK(ir->wall_ewald_zfac < 2);
1442 if ((ir->ewald_geometry == EwaldGeometry::ThreeDC) && (ir->pbcType != PbcType::XY)
1443 && EEL_FULL(ir->coulombtype))
1446 "With %s and ewald_geometry = %s you should use pbc = %s",
1447 enumValueToString(ir->coulombtype),
1448 enumValueToString(EwaldGeometry::ThreeDC),
1449 c_pbcTypeNames[PbcType::XY].c_str());
1450 warning(wi, warn_buf);
1452 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1454 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1455 CHECK(ir->bPeriodicMols);
1456 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1457 warning_note(wi, warn_buf);
1459 "With epsilon_surface > 0 you can only use domain decomposition "
1460 "when there are only small molecules with all bonds constrained (mdrun will check "
1462 warning_note(wi, warn_buf);
1465 if (ir_vdw_switched(ir))
1467 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1468 CHECK(ir->rvdw_switch >= ir->rvdw);
1470 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1473 "You are applying a switch function to vdw forces or potentials from %g to %g "
1474 "nm, which is more than half the interaction range, whereas switch functions "
1475 "are intended to act only close to the cut-off.",
1478 warning_note(wi, warn_buf);
1482 if (ir->vdwtype == VanDerWaalsType::Pme)
1484 if (!(ir->vdw_modifier == InteractionModifiers::None
1485 || ir->vdw_modifier == InteractionModifiers::PotShift))
1488 "With vdwtype = %s, the only supported modifiers are %s and %s",
1489 enumValueToString(ir->vdwtype),
1490 enumValueToString(InteractionModifiers::PotShift),
1491 enumValueToString(InteractionModifiers::None));
1492 warning_error(wi, err_buf);
1496 if (ir->vdwtype == VanDerWaalsType::User && ir->eDispCorr != DispersionCorrectionType::No)
1499 "You have selected user tables with dispersion correction, the dispersion "
1500 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1501 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1502 "really want dispersion correction to -C6/r^6.");
1505 if (ir->eI == IntegrationAlgorithm::LBFGS
1506 && (ir->coulombtype == CoulombInteractionType::Cut || ir->vdwtype == VanDerWaalsType::Cut)
1509 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1512 if (ir->eI == IntegrationAlgorithm::LBFGS && ir->nbfgscorr <= 0)
1514 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1517 /* IMPLICIT SOLVENT */
1518 if (ir->coulombtype == CoulombInteractionType::GBNotused)
1520 sprintf(warn_buf, "Invalid option %s for coulombtype", enumValueToString(ir->coulombtype));
1521 warning_error(wi, warn_buf);
1526 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1531 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1534 // cosine acceleration is only supported in leap-frog
1535 if (ir->cos_accel != 0.0 && ir->eI != IntegrationAlgorithm::MD)
1537 warning_error(wi, "cos-acceleration is only supported by integrator = md");
1541 /* interpret a number of doubles from a string and put them in an array,
1542 after allocating space for them.
1543 str = the input string
1544 n = the (pre-allocated) number of doubles read
1545 r = the output array of doubles. */
1546 static std::vector<real> parse_n_real(const std::string& str, int* n, warninp_t wi)
1548 auto values = gmx::splitString(str);
1551 std::vector<real> r;
1552 for (int i = 0; i < *n; i++)
1556 r.emplace_back(gmx::fromString<real>(values[i]));
1558 catch (gmx::GromacsException&)
1561 "Invalid value " + values[i]
1562 + " in string in mdp file. Expected a real number.");
1569 static void do_fep_params(t_inputrec* ir, gmx::ArrayRef<std::string> fep_lambda, char weights[STRLEN], warninp_t wi)
1572 int i, j, max_n_lambda, nweights;
1573 t_lambda* fep = ir->fepvals.get();
1574 t_expanded* expand = ir->expandedvals.get();
1575 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, std::vector<real>> count_fep_lambdas;
1576 bool bOneLambda = TRUE;
1577 gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, int> nfep;
1579 /* FEP input processing */
1580 /* first, identify the number of lambda values for each type.
1581 All that are nonzero must have the same number */
1583 for (auto i : keysOf(nfep))
1585 count_fep_lambdas[i] = parse_n_real(fep_lambda[static_cast<int>(i)], &(nfep[i]), wi);
1588 /* now, determine the number of components. All must be either zero, or equal. */
1591 for (auto i : keysOf(nfep))
1593 if (nfep[i] > max_n_lambda)
1595 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1596 must have the same number if its not zero.*/
1601 for (auto i : keysOf(nfep))
1605 ir->fepvals->separate_dvdl[i] = FALSE;
1607 else if (nfep[i] == max_n_lambda)
1609 if (i != FreeEnergyPerturbationCouplingType::Temperature) /* we treat this differently -- not really a reason to compute
1610 the derivative with respect to the temperature currently */
1612 ir->fepvals->separate_dvdl[i] = TRUE;
1618 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1621 enumValueToString(i),
1625 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1626 ir->fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Temperature] = FALSE;
1628 /* the number of lambdas is the number we've read in, which is either zero
1629 or the same for all */
1630 fep->n_lambda = max_n_lambda;
1632 /* if init_lambda is defined, we need to set lambda */
1633 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1635 ir->fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep] = TRUE;
1637 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1638 for (auto i : keysOf(nfep))
1640 fep->all_lambda[i].resize(fep->n_lambda);
1641 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1644 for (j = 0; j < fep->n_lambda; j++)
1646 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1651 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1652 internal bookkeeping -- for now, init_lambda */
1654 if ((nfep[FreeEnergyPerturbationCouplingType::Fep] == 0) && (fep->init_lambda >= 0))
1656 for (i = 0; i < fep->n_lambda; i++)
1658 fep->all_lambda[FreeEnergyPerturbationCouplingType::Fep][i] = fep->init_lambda;
1662 /* check to see if only a single component lambda is defined, and soft core is defined.
1663 In this case, turn on coulomb soft core */
1665 if (max_n_lambda == 0)
1671 for (auto i : keysOf(nfep))
1673 if ((nfep[i] != 0) && (i != FreeEnergyPerturbationCouplingType::Fep))
1679 if ((bOneLambda) && (fep->sc_alpha > 0))
1681 fep->bScCoul = TRUE;
1684 /* Fill in the others with the efptFEP if they are not explicitly
1685 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1686 they are all zero. */
1688 for (auto i : keysOf(nfep))
1690 if ((nfep[i] == 0) && (i != FreeEnergyPerturbationCouplingType::Fep))
1692 for (j = 0; j < fep->n_lambda; j++)
1694 fep->all_lambda[i][j] = fep->all_lambda[FreeEnergyPerturbationCouplingType::Fep][j];
1700 /* now read in the weights */
1701 expand->init_lambda_weights = parse_n_real(weights, &nweights, wi);
1704 expand->init_lambda_weights.resize(fep->n_lambda); /* initialize to zero */
1706 else if (nweights != fep->n_lambda)
1709 "Number of weights (%d) is not equal to number of lambda values (%d)",
1713 if ((expand->nstexpanded < 0) && (ir->efep != FreeEnergyPerturbationType::No))
1715 expand->nstexpanded = fep->nstdhdl;
1716 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1721 static void do_simtemp_params(t_inputrec* ir)
1723 ir->simtempvals->temperatures.resize(ir->fepvals->n_lambda);
1724 getSimTemps(ir->fepvals->n_lambda,
1725 ir->simtempvals.get(),
1726 ir->fepvals->all_lambda[FreeEnergyPerturbationCouplingType::Temperature]);
1729 template<typename T>
1730 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1733 for (const auto& input : inputs)
1737 outputs[i] = gmx::fromStdString<T>(input);
1739 catch (gmx::GromacsException&)
1741 auto message = gmx::formatString(
1742 "Invalid value for mdp option %s. %s should only consist of integers separated "
1746 warning_error(wi, message);
1752 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1755 for (const auto& input : inputs)
1759 outputs[i] = gmx::fromString<real>(input);
1761 catch (gmx::GromacsException&)
1763 auto message = gmx::formatString(
1764 "Invalid value for mdp option %s. %s should only consist of real numbers "
1765 "separated by spaces.",
1768 warning_error(wi, message);
1774 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1776 opts->wall_atomtype[0] = nullptr;
1777 opts->wall_atomtype[1] = nullptr;
1779 ir->wall_atomtype[0] = -1;
1780 ir->wall_atomtype[1] = -1;
1781 ir->wall_density[0] = 0;
1782 ir->wall_density[1] = 0;
1786 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1787 if (wallAtomTypes.size() != size_t(ir->nwall))
1790 "Expected %d elements for wall_atomtype, found %zu",
1792 wallAtomTypes.size());
1794 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1795 for (int i = 0; i < ir->nwall; i++)
1797 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1800 if (ir->wall_type == WallType::NineThree || ir->wall_type == WallType::TenFour)
1802 auto wallDensity = gmx::splitString(wall_density);
1803 if (wallDensity.size() != size_t(ir->nwall))
1806 "Expected %d elements for wall-density, found %zu",
1808 wallDensity.size());
1810 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1811 for (int i = 0; i < ir->nwall; i++)
1813 if (ir->wall_density[i] <= 0)
1815 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1822 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1826 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1827 for (int i = 0; i < nwall; i++)
1829 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1830 grps->emplace_back(groups->groupNames.size() - 1);
1835 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1837 /* read expanded ensemble parameters */
1838 printStringNewline(inp, "expanded ensemble variables");
1839 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1840 expand->elamstats = getEnum<LambdaWeightCalculation>(inp, "lmc-stats", wi);
1841 expand->elmcmove = getEnum<LambdaMoveCalculation>(inp, "lmc-move", wi);
1842 expand->elmceq = getEnum<LambdaWeightWillReachEquilibrium>(inp, "lmc-weights-equil", wi);
1843 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1844 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1845 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1846 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1847 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1848 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1849 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1850 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1851 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1852 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1853 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1854 expand->bSymmetrizedTMatrix =
1855 (getEnum<Boolean>(inp, "symmetrized-transition-matrix", wi) != Boolean::No);
1856 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1857 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1858 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1859 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1860 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1861 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1862 expand->bWLoneovert = (getEnum<Boolean>(inp, "wl-oneovert", wi) != Boolean::No);
1865 /*! \brief Return whether an end state with the given coupling-lambda
1866 * value describes fully-interacting VDW.
1868 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1869 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1871 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1873 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1879 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1882 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1884 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1886 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1889 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1890 std::string message = gmx::formatExceptionMessageToString(*ex);
1891 warning_error(wi_, message.c_str());
1896 std::string getOptionName(const gmx::KeyValueTreePath& context)
1898 if (mapping_ != nullptr)
1900 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1901 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1904 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1909 const gmx::IKeyValueTreeBackMapping* mapping_;
1914 void get_ir(const char* mdparin,
1915 const char* mdparout,
1916 gmx::MDModules* mdModules,
1919 WriteMdpHeader writeMdpHeader,
1923 double dumdub[2][6];
1925 char warn_buf[STRLEN];
1926 t_lambda* fep = ir->fepvals.get();
1927 t_expanded* expand = ir->expandedvals.get();
1929 const char* no_names[] = { "no", nullptr };
1931 init_inputrec_strings();
1932 gmx::TextInputFile stream(mdparin);
1933 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1935 snew(dumstr[0], STRLEN);
1936 snew(dumstr[1], STRLEN);
1938 /* ignore the following deprecated commands */
1939 replace_inp_entry(inp, "title", nullptr);
1940 replace_inp_entry(inp, "cpp", nullptr);
1941 replace_inp_entry(inp, "domain-decomposition", nullptr);
1942 replace_inp_entry(inp, "andersen-seed", nullptr);
1943 replace_inp_entry(inp, "dihre", nullptr);
1944 replace_inp_entry(inp, "dihre-fc", nullptr);
1945 replace_inp_entry(inp, "dihre-tau", nullptr);
1946 replace_inp_entry(inp, "nstdihreout", nullptr);
1947 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1948 replace_inp_entry(inp, "optimize-fft", nullptr);
1949 replace_inp_entry(inp, "adress_type", nullptr);
1950 replace_inp_entry(inp, "adress_const_wf", nullptr);
1951 replace_inp_entry(inp, "adress_ex_width", nullptr);
1952 replace_inp_entry(inp, "adress_hy_width", nullptr);
1953 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1954 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1955 replace_inp_entry(inp, "adress_site", nullptr);
1956 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1957 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1958 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1959 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1960 replace_inp_entry(inp, "rlistlong", nullptr);
1961 replace_inp_entry(inp, "nstcalclr", nullptr);
1962 replace_inp_entry(inp, "pull-print-com2", nullptr);
1963 replace_inp_entry(inp, "gb-algorithm", nullptr);
1964 replace_inp_entry(inp, "nstgbradii", nullptr);
1965 replace_inp_entry(inp, "rgbradii", nullptr);
1966 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1967 replace_inp_entry(inp, "gb-saltconc", nullptr);
1968 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1969 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1970 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1971 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1972 replace_inp_entry(inp, "sa-algorithm", nullptr);
1973 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1974 replace_inp_entry(inp, "ns-type", nullptr);
1976 /* replace the following commands with the clearer new versions*/
1977 replace_inp_entry(inp, "unconstrained-start", "continuation");
1978 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1979 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1980 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1981 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1982 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1983 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1985 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1986 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1987 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1988 setStringEntry(&inp, "include", opts->include, nullptr);
1989 printStringNoNewline(
1990 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1991 setStringEntry(&inp, "define", opts->define, nullptr);
1993 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1994 ir->eI = getEnum<IntegrationAlgorithm>(&inp, "integrator", wi);
1995 printStringNoNewline(&inp, "Start time and timestep in ps");
1996 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1997 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1998 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1999 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
2000 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
2001 printStringNoNewline(
2002 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
2003 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
2004 printStringNoNewline(&inp, "Multiple time-stepping");
2005 ir->useMts = (getEnum<Boolean>(&inp, "mts", wi) != Boolean::No);
2008 gmx::GromppMtsOpts& mtsOpts = opts->mtsOpts;
2009 mtsOpts.numLevels = get_eint(&inp, "mts-levels", 2, wi);
2010 mtsOpts.level2Forces = setStringEntry(&inp, "mts-level2-forces", "longrange-nonbonded");
2011 mtsOpts.level2Factor = get_eint(&inp, "mts-level2-factor", 2, wi);
2013 // We clear after reading without dynamics to not force the user to remove MTS mdp options
2014 if (!EI_DYNAMICS(ir->eI))
2019 printStringNoNewline(&inp, "mode for center of mass motion removal");
2020 ir->comm_mode = getEnum<ComRemovalAlgorithm>(&inp, "comm-mode", wi);
2021 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
2022 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
2023 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
2024 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
2026 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2027 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2028 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2029 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2032 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2033 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2034 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2035 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2036 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2037 ir->niter = get_eint(&inp, "niter", 20, wi);
2038 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2039 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2040 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2041 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2042 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2044 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2045 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2047 /* Output options */
2048 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2049 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2050 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2051 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2052 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2053 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2054 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2055 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2056 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2057 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2058 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2059 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2060 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2061 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2062 printStringNoNewline(&inp, "default, all atoms will be written.");
2063 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2064 printStringNoNewline(&inp, "Selection of energy groups");
2065 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2067 /* Neighbor searching */
2068 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2069 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2070 ir->cutoff_scheme = getEnum<CutoffScheme>(&inp, "cutoff-scheme", wi);
2071 printStringNoNewline(&inp, "nblist update frequency");
2072 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2073 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2074 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2075 std::vector<const char*> pbcTypesNamesChar;
2076 for (const auto& pbcTypeName : c_pbcTypeNames)
2078 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2080 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2081 ir->bPeriodicMols = getEnum<Boolean>(&inp, "periodic-molecules", wi) != Boolean::No;
2082 printStringNoNewline(&inp,
2083 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2084 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2085 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2086 printStringNoNewline(&inp, "nblist cut-off");
2087 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2088 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2090 /* Electrostatics */
2091 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2092 printStringNoNewline(&inp, "Method for doing electrostatics");
2093 ir->coulombtype = getEnum<CoulombInteractionType>(&inp, "coulombtype", wi);
2094 ir->coulomb_modifier = getEnum<InteractionModifiers>(&inp, "coulomb-modifier", wi);
2095 printStringNoNewline(&inp, "cut-off lengths");
2096 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2097 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2098 printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
2099 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2100 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2101 printStringNoNewline(&inp, "Method for doing Van der Waals");
2102 ir->vdwtype = getEnum<VanDerWaalsType>(&inp, "vdw-type", wi);
2103 ir->vdw_modifier = getEnum<InteractionModifiers>(&inp, "vdw-modifier", wi);
2104 printStringNoNewline(&inp, "cut-off lengths");
2105 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2106 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2107 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2108 ir->eDispCorr = getEnum<DispersionCorrectionType>(&inp, "DispCorr", wi);
2109 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2110 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2111 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2112 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2113 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2114 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2115 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2116 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2117 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2118 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2119 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2120 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2121 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2122 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2123 ir->ljpme_combination_rule = getEnum<LongRangeVdW>(&inp, "lj-pme-comb-rule", wi);
2124 ir->ewald_geometry = getEnum<EwaldGeometry>(&inp, "ewald-geometry", wi);
2125 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2127 /* Implicit solvation is no longer supported, but we need grompp
2128 to be able to refuse old .mdp files that would have built a tpr
2129 to run it. Thus, only "no" is accepted. */
2130 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2132 /* Coupling stuff */
2133 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2134 printStringNoNewline(&inp, "Temperature coupling");
2135 ir->etc = getEnum<TemperatureCoupling>(&inp, "tcoupl", wi);
2136 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2137 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2138 ir->bPrintNHChains = (getEnum<Boolean>(&inp, "print-nose-hoover-chain-variables", wi) != Boolean::No);
2139 printStringNoNewline(&inp, "Groups to couple separately");
2140 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2141 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2142 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2143 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2144 printStringNoNewline(&inp, "pressure coupling");
2145 ir->epc = getEnum<PressureCoupling>(&inp, "pcoupl", wi);
2146 ir->epct = getEnum<PressureCouplingType>(&inp, "pcoupltype", wi);
2147 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2148 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2149 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2150 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2151 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2152 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2153 ir->refcoord_scaling = getEnum<RefCoordScaling>(&inp, "refcoord-scaling", wi);
2156 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2157 ir->bQMMM = (getEnum<Boolean>(&inp, "QMMM", wi) != Boolean::No);
2158 printStringNoNewline(&inp, "Groups treated with MiMiC");
2159 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2161 /* Simulated annealing */
2162 printStringNewline(&inp, "SIMULATED ANNEALING");
2163 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2164 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2165 printStringNoNewline(&inp,
2166 "Number of time points to use for specifying annealing in each group");
2167 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2168 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2169 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2170 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2171 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2174 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2175 opts->bGenVel = (getEnum<Boolean>(&inp, "gen-vel", wi) != Boolean::No);
2176 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2177 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2180 printStringNewline(&inp, "OPTIONS FOR BONDS");
2181 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2182 printStringNoNewline(&inp, "Type of constraint algorithm");
2183 ir->eConstrAlg = getEnum<ConstraintAlgorithm>(&inp, "constraint-algorithm", wi);
2184 printStringNoNewline(&inp, "Do not constrain the start configuration");
2185 ir->bContinuation = (getEnum<Boolean>(&inp, "continuation", wi) != Boolean::No);
2186 printStringNoNewline(&inp,
2187 "Use successive overrelaxation to reduce the number of shake iterations");
2188 ir->bShakeSOR = (getEnum<Boolean>(&inp, "Shake-SOR", wi) != Boolean::No);
2189 printStringNoNewline(&inp, "Relative tolerance of shake");
2190 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2191 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2192 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2193 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2194 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2195 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2196 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2197 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2198 printStringNoNewline(&inp, "rotates over more degrees than");
2199 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2200 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2201 opts->bMorse = (getEnum<Boolean>(&inp, "morse", wi) != Boolean::No);
2203 /* Energy group exclusions */
2204 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2205 printStringNoNewline(
2206 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2207 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2210 printStringNewline(&inp, "WALLS");
2211 printStringNoNewline(
2212 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2213 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2214 ir->wall_type = getEnum<WallType>(&inp, "wall-type", wi);
2215 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2216 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2217 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2218 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2221 printStringNewline(&inp, "COM PULLING");
2222 ir->bPull = (getEnum<Boolean>(&inp, "pull", wi) != Boolean::No);
2225 ir->pull = std::make_unique<pull_params_t>();
2226 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
2230 for (int c = 0; c < ir->pull->ncoord; c++)
2232 if (ir->pull->coord[c].eType == PullingAlgorithm::Constraint)
2235 "Constraint COM pulling is not supported in combination with "
2236 "multiple time stepping");
2244 NOTE: needs COM pulling or free energy input */
2245 printStringNewline(&inp, "AWH biasing");
2246 ir->bDoAwh = (getEnum<Boolean>(&inp, "awh", wi) != Boolean::No);
2249 ir->awhParams = std::make_unique<gmx::AwhParams>(&inp, wi);
2252 /* Enforced rotation */
2253 printStringNewline(&inp, "ENFORCED ROTATION");
2254 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2255 ir->bRot = (getEnum<Boolean>(&inp, "rotation", wi) != Boolean::No);
2259 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2262 /* Interactive MD */
2264 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2265 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2266 if (inputrecStrings->imd_grp[0] != '\0')
2273 printStringNewline(&inp, "NMR refinement stuff");
2274 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2275 ir->eDisre = getEnum<DistanceRestraintRefinement>(&inp, "disre", wi);
2276 printStringNoNewline(
2277 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2278 ir->eDisreWeighting = getEnum<DistanceRestraintWeighting>(&inp, "disre-weighting", wi);
2279 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2280 ir->bDisreMixed = (getEnum<Boolean>(&inp, "disre-mixed", wi) != Boolean::No);
2281 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2282 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2283 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2284 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2285 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2286 opts->bOrire = (getEnum<Boolean>(&inp, "orire", wi) != Boolean::No);
2287 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2288 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2289 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2290 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2291 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2292 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2294 /* free energy variables */
2295 printStringNewline(&inp, "Free energy variables");
2296 ir->efep = getEnum<FreeEnergyPerturbationType>(&inp, "free-energy", wi);
2297 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2298 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2299 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2300 opts->bCoupleIntra = (getEnum<Boolean>(&inp, "couple-intramol", wi) != Boolean::No);
2302 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2304 it was not entered */
2305 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2306 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2307 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2308 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Fep] =
2309 setStringEntry(&inp, "fep-lambdas", "");
2310 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Mass] =
2311 setStringEntry(&inp, "mass-lambdas", "");
2312 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Coul] =
2313 setStringEntry(&inp, "coul-lambdas", "");
2314 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Vdw] =
2315 setStringEntry(&inp, "vdw-lambdas", "");
2316 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Bonded] =
2317 setStringEntry(&inp, "bonded-lambdas", "");
2318 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Restraint] =
2319 setStringEntry(&inp, "restraint-lambdas", "");
2320 inputrecStrings->fep_lambda[FreeEnergyPerturbationCouplingType::Temperature] =
2321 setStringEntry(&inp, "temperature-lambdas", "");
2322 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2323 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2324 fep->edHdLPrintEnergy = getEnum<FreeEnergyPrintEnergy>(&inp, "dhdl-print-energy", wi);
2325 fep->softcoreFunction = getEnum<SoftcoreType>(&inp, "sc-function", wi);
2326 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2327 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2328 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2329 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2330 fep->bScCoul = (getEnum<Boolean>(&inp, "sc-coul", wi) != Boolean::No);
2331 fep->scGapsysScaleLinpointLJ = get_ereal(&inp, "sc-gapsys-scale-linpoint-lj", 0.85, wi);
2332 fep->scGapsysScaleLinpointQ = get_ereal(&inp, "sc-gapsys-scale-linpoint-q", 0.3, wi);
2333 fep->scGapsysSigmaLJ = get_ereal(&inp, "sc-gapsys-sigma-lj", 0.3, wi);
2334 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2335 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2336 fep->separate_dhdl_file = getEnum<SeparateDhdlFile>(&inp, "separate-dhdl-file", wi);
2337 fep->dhdl_derivatives = getEnum<DhDlDerivativeCalculation>(&inp, "dhdl-derivatives", wi);
2338 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2339 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2341 /* Non-equilibrium MD stuff */
2342 printStringNewline(&inp, "Non-equilibrium MD stuff");
2343 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2344 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2345 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2346 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2348 /* simulated tempering variables */
2349 printStringNewline(&inp, "simulated tempering variables");
2350 ir->bSimTemp = (getEnum<Boolean>(&inp, "simulated-tempering", wi) != Boolean::No);
2351 ir->simtempvals->eSimTempScale = getEnum<SimulatedTempering>(&inp, "simulated-tempering-scaling", wi);
2352 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2353 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2355 /* expanded ensemble variables */
2356 if (ir->efep == FreeEnergyPerturbationType::Expanded || ir->bSimTemp)
2358 read_expandedparams(&inp, expand, wi);
2361 /* Electric fields */
2363 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2364 gmx::KeyValueTreeTransformer transform;
2365 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2366 mdModules->initMdpTransform(transform.rules());
2367 for (const auto& path : transform.mappedPaths())
2369 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2370 mark_einp_set(inp, path[0].c_str());
2372 MdpErrorHandler errorHandler(wi);
2373 auto result = transform.transform(convertedValues, &errorHandler);
2374 ir->params = new gmx::KeyValueTreeObject(result.object());
2375 mdModules->adjustInputrecBasedOnModules(ir);
2376 errorHandler.setBackMapping(result.backMapping());
2377 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2380 /* Ion/water position swapping ("computational electrophysiology") */
2381 printStringNewline(&inp,
2382 "Ion/water position swapping for computational electrophysiology setups");
2383 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2384 ir->eSwapCoords = getEnum<SwapType>(&inp, "swapcoords", wi);
2385 if (ir->eSwapCoords != SwapType::No)
2392 printStringNoNewline(&inp, "Swap attempt frequency");
2393 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2394 printStringNoNewline(&inp, "Number of ion types to be controlled");
2395 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2398 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2400 ir->swap->ngrp = nIonTypes + static_cast<int>(SwapGroupSplittingType::Count);
2401 snew(ir->swap->grp, ir->swap->ngrp);
2402 for (i = 0; i < ir->swap->ngrp; i++)
2404 snew(ir->swap->grp[i].molname, STRLEN);
2406 printStringNoNewline(&inp,
2407 "Two index groups that contain the compartment-partitioning atoms");
2408 setStringEntry(&inp,
2410 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
2412 setStringEntry(&inp,
2414 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname,
2416 printStringNoNewline(&inp,
2417 "Use center of mass of split groups (yes/no), otherwise center of "
2418 "geometry is used");
2419 ir->swap->massw_split[0] = (getEnum<Boolean>(&inp, "massw-split0", wi) != Boolean::No);
2420 ir->swap->massw_split[1] = (getEnum<Boolean>(&inp, "massw-split1", wi) != Boolean::No);
2422 printStringNoNewline(&inp, "Name of solvent molecules");
2423 setStringEntry(&inp,
2425 ir->swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].molname,
2428 printStringNoNewline(&inp,
2429 "Split cylinder: radius, upper and lower extension (nm) (this will "
2430 "define the channels)");
2431 printStringNoNewline(&inp,
2432 "Note that the split cylinder settings do not have an influence on "
2433 "the swapping protocol,");
2434 printStringNoNewline(
2436 "however, if correctly defined, the permeation events are recorded per channel");
2437 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2438 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2439 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2440 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2441 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2442 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2444 printStringNoNewline(
2446 "Average the number of ions per compartment over these many swap attempt steps");
2447 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2449 printStringNoNewline(
2450 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2451 printStringNoNewline(
2452 &inp, "and the requested number of ions of this type in compartments A and B");
2453 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2454 for (i = 0; i < nIonTypes; i++)
2456 int ig = static_cast<int>(SwapGroupSplittingType::Count) + i;
2458 sprintf(buf, "iontype%d-name", i);
2459 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2460 sprintf(buf, "iontype%d-in-A", i);
2461 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2462 sprintf(buf, "iontype%d-in-B", i);
2463 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2466 printStringNoNewline(
2468 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2469 printStringNoNewline(
2471 "at maximum distance (= bulk concentration) to the split group layers. However,");
2472 printStringNoNewline(&inp,
2473 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2474 "layer from the middle at 0.0");
2475 printStringNoNewline(&inp,
2476 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2477 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2478 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2479 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2480 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2482 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2485 printStringNoNewline(
2486 &inp, "Start to swap ions if threshold difference to requested count is reached");
2487 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2490 /* AdResS is no longer supported, but we need grompp to be able to
2491 refuse to process old .mdp files that used it. */
2492 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2494 /* User defined thingies */
2495 printStringNewline(&inp, "User defined thingies");
2496 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2497 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2498 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2499 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2500 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2501 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2502 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2503 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2504 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2505 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2510 gmx::TextOutputFile stream(mdparout);
2511 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2513 // Transform module data into a flat key-value tree for output.
2514 gmx::KeyValueTreeBuilder builder;
2515 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2516 mdModules->buildMdpOutput(&builderObject);
2518 gmx::TextWriter writer(&stream);
2519 writeKeyValueTreeAsMdp(&writer, builder.build());
2524 /* Process options if necessary */
2525 for (m = 0; m < 2; m++)
2527 for (i = 0; i < 2 * DIM; i++)
2531 if (ir->epc != PressureCoupling::No)
2535 case PressureCouplingType::Isotropic:
2536 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2540 "Pressure coupling incorrect number of values (I need exactly 1)");
2542 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2544 case PressureCouplingType::SemiIsotropic:
2545 case PressureCouplingType::SurfaceTension:
2546 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2550 "Pressure coupling incorrect number of values (I need exactly 2)");
2552 dumdub[m][YY] = dumdub[m][XX];
2554 case PressureCouplingType::Anisotropic:
2555 if (sscanf(dumstr[m],
2556 "%lf%lf%lf%lf%lf%lf",
2567 "Pressure coupling incorrect number of values (I need exactly 6)");
2572 "Pressure coupling type %s not implemented yet",
2573 enumValueToString(ir->epct));
2577 clear_mat(ir->ref_p);
2578 clear_mat(ir->compress);
2579 for (i = 0; i < DIM; i++)
2581 ir->ref_p[i][i] = dumdub[1][i];
2582 ir->compress[i][i] = dumdub[0][i];
2584 if (ir->epct == PressureCouplingType::Anisotropic)
2586 ir->ref_p[XX][YY] = dumdub[1][3];
2587 ir->ref_p[XX][ZZ] = dumdub[1][4];
2588 ir->ref_p[YY][ZZ] = dumdub[1][5];
2589 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2592 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2593 "apply a threefold shear stress?\n");
2595 ir->compress[XX][YY] = dumdub[0][3];
2596 ir->compress[XX][ZZ] = dumdub[0][4];
2597 ir->compress[YY][ZZ] = dumdub[0][5];
2598 for (i = 0; i < DIM; i++)
2600 for (m = 0; m < i; m++)
2602 ir->ref_p[i][m] = ir->ref_p[m][i];
2603 ir->compress[i][m] = ir->compress[m][i];
2608 if (ir->comm_mode == ComRemovalAlgorithm::No)
2613 opts->couple_moltype = nullptr;
2614 if (strlen(inputrecStrings->couple_moltype) > 0)
2616 if (ir->efep != FreeEnergyPerturbationType::No)
2618 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2619 if (opts->couple_lam0 == opts->couple_lam1)
2621 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2623 if (ir->eI == IntegrationAlgorithm::MD
2624 && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2628 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2635 "Free energy is turned off, so we will not decouple the molecule listed "
2639 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2640 if (ir->efep != FreeEnergyPerturbationType::No)
2642 if (fep->delta_lambda != 0)
2644 ir->efep = FreeEnergyPerturbationType::SlowGrowth;
2648 if (fep->edHdLPrintEnergy == FreeEnergyPrintEnergy::Yes)
2650 fep->edHdLPrintEnergy = FreeEnergyPrintEnergy::Total;
2652 "Old option for dhdl-print-energy given: "
2653 "changing \"yes\" to \"total\"\n");
2656 if (ir->bSimTemp && (fep->edHdLPrintEnergy == FreeEnergyPrintEnergy::No))
2658 /* always print out the energy to dhdl if we are doing
2659 expanded ensemble, since we need the total energy for
2660 analysis if the temperature is changing. In some
2661 conditions one may only want the potential energy, so
2662 we will allow that if the appropriate mdp setting has
2663 been enabled. Otherwise, total it is:
2665 fep->edHdLPrintEnergy = FreeEnergyPrintEnergy::Total;
2668 if ((ir->efep != FreeEnergyPerturbationType::No) || ir->bSimTemp)
2670 ir->bExpanded = FALSE;
2671 if ((ir->efep == FreeEnergyPerturbationType::Expanded) || ir->bSimTemp)
2673 ir->bExpanded = TRUE;
2675 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2676 if (ir->bSimTemp) /* done after fep params */
2678 do_simtemp_params(ir);
2681 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2682 * setup and not on the old way of specifying the free-energy setup,
2683 * we should check for using soft-core when not needed, since that
2684 * can complicate the sampling significantly.
2685 * Note that we only check for the automated coupling setup.
2686 * If the (advanced) user does FEP through manual topology changes,
2687 * this check will not be triggered.
2689 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->n_lambda == 0
2690 && ir->fepvals->sc_alpha != 0
2691 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2694 "You are using soft-core interactions while the Van der Waals interactions are "
2695 "not decoupled (note that the sc-coul option is only active when using lambda "
2696 "states). Although this will not lead to errors, you will need much more "
2697 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2702 ir->fepvals->n_lambda = 0;
2705 /* WALL PARAMETERS */
2707 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2709 /* ORIENTATION RESTRAINT PARAMETERS */
2711 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2713 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2716 /* DEFORMATION PARAMETERS */
2718 clear_mat(ir->deform);
2719 for (i = 0; i < 6; i++)
2724 double gmx_unused canary;
2725 int ndeform = sscanf(inputrecStrings->deform,
2726 "%lf %lf %lf %lf %lf %lf %lf",
2735 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2739 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2740 inputrecStrings->deform)
2743 for (i = 0; i < 3; i++)
2745 ir->deform[i][i] = dumdub[0][i];
2747 ir->deform[YY][XX] = dumdub[0][3];
2748 ir->deform[ZZ][XX] = dumdub[0][4];
2749 ir->deform[ZZ][YY] = dumdub[0][5];
2750 if (ir->epc != PressureCoupling::No)
2752 for (i = 0; i < 3; i++)
2754 for (j = 0; j <= i; j++)
2756 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2758 warning_error(wi, "A box element has deform set and compressibility > 0");
2762 for (i = 0; i < 3; i++)
2764 for (j = 0; j < i; j++)
2766 if (ir->deform[i][j] != 0)
2768 for (m = j; m < DIM; m++)
2770 if (ir->compress[m][j] != 0)
2773 "An off-diagonal box element has deform set while "
2774 "compressibility > 0 for the same component of another box "
2775 "vector, this might lead to spurious periodicity effects.");
2776 warning(wi, warn_buf);
2784 /* Ion/water position swapping checks */
2785 if (ir->eSwapCoords != SwapType::No)
2787 if (ir->swap->nstswap < 1)
2789 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2791 if (ir->swap->nAverage < 1)
2793 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2795 if (ir->swap->threshold < 1.0)
2797 warning_error(wi, "Ion count threshold must be at least 1.\n");
2801 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2804 std::vector<std::string> errorMessages;
2805 ir->mtsLevels = gmx::setupMtsLevels(opts->mtsOpts, &errorMessages);
2807 for (const auto& errorMessage : errorMessages)
2809 warning_error(wi, errorMessage.c_str());
2815 gmx::checkAwhParams(*ir->awhParams, *ir, wi);
2822 int search_string(const char* s, int ng, char* const gn[])
2826 for (i = 0; (i < ng); i++)
2828 if (gmx_strcasecmp(s, gn[i]) == 0)
2835 "Group %s referenced in the .mdp file was not found in the index file.\n"
2836 "Group names must match either [moleculetype] names or custom index group\n"
2837 "names, in which case you must supply an index file to the '-n' option\n"
2842 static void atomGroupRangeValidation(int natoms, int groupIndex, const t_blocka& block)
2844 /* Now go over the atoms in the group */
2845 for (int j = block.index[groupIndex]; (j < block.index[groupIndex + 1]); j++)
2847 int aj = block.a[j];
2849 /* Range checking */
2850 if ((aj < 0) || (aj >= natoms))
2852 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2857 /*! Creates the groups of atom indices for group type \p gtype
2859 * \param[in] natoms The total number of atoms in the system
2860 * \param[in,out] groups Index \p gtype in this list of list of groups will be set
2861 * \param[in] groupsFromMdpFile The list of group names set for \p gtype in the mdp file
2862 * \param[in] block The list of atom indices for all available index groups
2863 * \param[in] gnames The list of names for all available index groups
2864 * \param[in] gtype The group type to creates groups for
2865 * \param[in] restnm The index of rest group name in \p gnames
2866 * \param[in] coverage How to treat coverage of all atoms in the system
2867 * \param[in] bVerbose Whether to print when we make a rest group
2868 * \param[in,out] wi List of warnings
2870 static void do_numbering(const int natoms,
2871 SimulationGroups* groups,
2872 gmx::ArrayRef<const std::string> groupsFromMdpFile,
2873 const t_blocka* block,
2874 char* const gnames[],
2875 const SimulationAtomGroupType gtype,
2877 const GroupCoverage coverage,
2878 const bool bVerbose,
2881 unsigned short* cbuf;
2882 AtomGroupIndices* grps = &(groups->groups[gtype]);
2885 char warn_buf[STRLEN];
2887 title = shortName(gtype);
2890 /* Mark all id's as not set */
2891 for (int i = 0; (i < natoms); i++)
2896 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2898 /* Lookup the group name in the block structure */
2899 const int gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2900 if ((coverage != GroupCoverage::OneGroup) || (i == 0))
2902 grps->emplace_back(gid);
2904 GMX_ASSERT(block, "Can't have a nullptr block");
2905 atomGroupRangeValidation(natoms, gid, *block);
2906 /* Now go over the atoms in the group */
2907 for (int j = block->index[gid]; (j < block->index[gid + 1]); j++)
2909 const int aj = block->a[j];
2910 /* Lookup up the old group number */
2911 const int ognr = cbuf[aj];
2914 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title, ognr + 1, i + 1);
2918 /* Store the group number in buffer */
2919 if (coverage == GroupCoverage::OneGroup)
2932 /* Now check whether we have done all atoms */
2935 if (coverage == GroupCoverage::All)
2937 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2939 else if (coverage == GroupCoverage::Partial)
2941 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2942 warning_note(wi, warn_buf);
2944 /* Assign all atoms currently unassigned to a rest group */
2945 for (int j = 0; (j < natoms); j++)
2947 if (cbuf[j] == NOGID)
2949 cbuf[j] = grps->size();
2952 if (coverage != GroupCoverage::Partial)
2956 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title, natoms - ntot);
2958 /* Add group name "rest" */
2959 grps->emplace_back(restnm);
2961 /* Assign the rest name to all atoms not currently assigned to a group */
2962 for (int j = 0; (j < natoms); j++)
2964 if (cbuf[j] == NOGID)
2966 // group size was not updated before this here, so need to use -1.
2967 cbuf[j] = grps->size() - 1;
2973 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2975 /* All atoms are part of one (or no) group, no index required */
2976 groups->groupNumbers[gtype].clear();
2980 for (int j = 0; (j < natoms); j++)
2982 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2989 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2992 pull_params_t* pull;
2993 int natoms, imin, jmin;
2994 int * nrdf2, *na_vcm, na_tot;
2995 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
3000 * First calc 3xnr-atoms for each group
3001 * then subtract half a degree of freedom for each constraint
3003 * Only atoms and nuclei contribute to the degrees of freedom...
3008 const SimulationGroups& groups = mtop->groups;
3009 natoms = mtop->natoms;
3011 /* Allocate one more for a possible rest group */
3012 /* We need to sum degrees of freedom into doubles,
3013 * since floats give too low nrdf's above 3 million atoms.
3015 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
3016 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
3017 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
3018 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
3019 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
3021 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
3025 for (gmx::index i = 0;
3026 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3030 clear_ivec(dof_vcm[i]);
3032 nrdf_vcm_sub[i] = 0;
3034 snew(nrdf2, natoms);
3035 for (const AtomProxy atomP : AtomRange(*mtop))
3037 const t_atom& local = atomP.atom();
3038 int i = atomP.globalAtomNumber();
3040 if (local.ptype == ParticleType::Atom || local.ptype == ParticleType::Nucleus)
3042 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
3043 for (int d = 0; d < DIM; d++)
3045 if (opts->nFreeze[g][d] == 0)
3047 /* Add one DOF for particle i (counted as 2*1) */
3049 /* VCM group i has dim d as a DOF */
3050 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
3054 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
3056 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
3062 for (const gmx_molblock_t& molb : mtop->molblock)
3064 const gmx_moltype_t& molt = mtop->moltype[molb.type];
3065 const t_atom* atom = molt.atoms.atom;
3066 for (int mol = 0; mol < molb.nmol; mol++)
3068 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
3070 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
3071 for (int i = 0; i < molt.ilist[ftype].size();)
3073 /* Subtract degrees of freedom for the constraints,
3074 * if the particles still have degrees of freedom left.
3075 * If one of the particles is a vsite or a shell, then all
3076 * constraint motion will go there, but since they do not
3077 * contribute to the constraints the degrees of freedom do not
3080 int ai = as + ia[i + 1];
3081 int aj = as + ia[i + 2];
3082 if (((atom[ia[i + 1]].ptype == ParticleType::Nucleus)
3083 || (atom[ia[i + 1]].ptype == ParticleType::Atom))
3084 && ((atom[ia[i + 2]].ptype == ParticleType::Nucleus)
3085 || (atom[ia[i + 2]].ptype == ParticleType::Atom)))
3103 imin = std::min(imin, nrdf2[ai]);
3104 jmin = std::min(jmin, nrdf2[aj]);
3107 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3109 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3111 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3113 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3116 i += interaction_function[ftype].nratoms + 1;
3119 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3120 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3122 /* Subtract 1 dof from every atom in the SETTLE */
3123 for (int j = 0; j < 3; j++)
3125 int ai = as + ia[i + 1 + j];
3126 imin = std::min(2, nrdf2[ai]);
3128 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3130 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3135 as += molt.atoms.nr;
3141 /* Correct nrdf for the COM constraints.
3142 * We correct using the TC and VCM group of the first atom
3143 * in the reference and pull group. If atoms in one pull group
3144 * belong to different TC or VCM groups it is anyhow difficult
3145 * to determine the optimal nrdf assignment.
3147 pull = ir->pull.get();
3149 for (int i = 0; i < pull->ncoord; i++)
3151 if (pull->coord[i].eType != PullingAlgorithm::Constraint)
3158 for (int j = 0; j < 2; j++)
3160 const t_pull_group* pgrp;
3162 pgrp = &pull->group[pull->coord[i].group[j]];
3164 if (!pgrp->ind.empty())
3166 /* Subtract 1/2 dof from each group */
3167 int ai = pgrp->ind[0];
3168 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3170 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3172 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3175 "Center of mass pulling constraints caused the number of degrees "
3176 "of freedom for temperature coupling group %s to be negative",
3177 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3178 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3183 /* We need to subtract the whole DOF from group j=1 */
3190 if (ir->nstcomm != 0)
3192 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3193 "Expect at least one group when removing COM motion");
3195 /* We remove COM motion up to dim ndof_com() */
3196 const int ndim_rm_vcm = ndof_com(ir);
3198 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3199 * the number of degrees of freedom in each vcm group when COM
3200 * translation is removed and 6 when rotation is removed as well.
3201 * Note that we do not and should not include the rest group here.
3203 for (gmx::index j = 0;
3204 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]);
3207 switch (ir->comm_mode)
3209 case ComRemovalAlgorithm::Linear:
3210 case ComRemovalAlgorithm::LinearAccelerationCorrection:
3211 nrdf_vcm_sub[j] = 0;
3212 for (int d = 0; d < ndim_rm_vcm; d++)
3220 case ComRemovalAlgorithm::Angular: nrdf_vcm_sub[j] = 6; break;
3221 default: gmx_incons("Checking comm_mode");
3225 for (gmx::index i = 0;
3226 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]);
3229 /* Count the number of atoms of TC group i for every VCM group */
3230 for (gmx::index j = 0;
3231 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3237 for (int ai = 0; ai < natoms; ai++)
3239 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3241 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3245 /* Correct for VCM removal according to the fraction of each VCM
3246 * group present in this TC group.
3248 nrdf_uc = nrdf_tc[i];
3250 for (gmx::index j = 0;
3251 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3254 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3256 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3257 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3262 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3264 opts->nrdf[i] = nrdf_tc[i];
3265 if (opts->nrdf[i] < 0)
3270 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3271 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]],
3280 sfree(nrdf_vcm_sub);
3283 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3285 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3286 * But since this is much larger than STRLEN, such a line can not be parsed.
3287 * The real maximum is the number of names that fit in a string: STRLEN/2.
3292 auto names = gmx::splitString(val);
3293 if (names.size() % 2 != 0)
3295 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3297 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3299 for (size_t i = 0; i < names.size() / 2; i++)
3301 // TODO this needs to be replaced by a solution using std::find_if
3305 names[2 * i].c_str(),
3306 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3312 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3317 names[2 * i + 1].c_str(),
3318 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3324 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3326 if ((j < nr) && (k < nr))
3328 ir->opts.egp_flags[nr * j + k] |= flag;
3329 ir->opts.egp_flags[nr * k + j] |= flag;
3338 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3340 int ig = -1, i = 0, gind;
3344 /* Just a quick check here, more thorough checks are in mdrun */
3345 if (strcmp(swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname,
3346 swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname)
3350 "The split groups can not both be '%s'.",
3351 swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname);
3354 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3355 for (ig = 0; ig < swap->ngrp; ig++)
3357 swapg = &swap->grp[ig];
3358 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3359 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3364 "%s group '%s' contains %d atoms.\n",
3365 ig < 3 ? enumValueToString(static_cast<SwapGroupSplittingType>(ig)) : "Swap",
3366 swap->grp[ig].molname,
3368 snew(swapg->ind, swapg->nat);
3369 for (i = 0; i < swapg->nat; i++)
3371 swapg->ind[i] = grps->a[grps->index[gind] + i];
3376 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3382 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3387 ig = search_string(IMDgname, grps->nr, gnames);
3388 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3390 if (IMDgroup->nat > 0)
3393 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3397 snew(IMDgroup->ind, IMDgroup->nat);
3398 for (i = 0; i < IMDgroup->nat; i++)
3400 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3405 /* Checks whether atoms are both part of a COM removal group and frozen.
3406 * If a fully frozen atom is part of a COM removal group, it is removed
3407 * from the COM removal group. A note is issued if such atoms are present.
3408 * A warning is issued for atom with one or two dimensions frozen that
3409 * are part of a COM removal group (mdrun would need to compute COM mass
3410 * per dimension to handle this correctly).
3411 * Also issues a warning when non-frozen atoms are not part of a COM
3412 * removal group while COM removal is active.
3414 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3416 const t_grpopts& opts,
3419 const int vcmRestGroup =
3420 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3422 int numFullyFrozenVcmAtoms = 0;
3423 int numPartiallyFrozenVcmAtoms = 0;
3424 int numNonVcmAtoms = 0;
3425 for (int a = 0; a < numAtoms; a++)
3427 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3428 int numFrozenDims = 0;
3429 for (int d = 0; d < DIM; d++)
3431 numFrozenDims += opts.nFreeze[freezeGroup][d];
3434 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3435 if (vcmGroup < vcmRestGroup)
3437 if (numFrozenDims == DIM)
3439 /* Do not remove COM motion for this fully frozen atom */
3440 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3442 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3445 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3446 numFullyFrozenVcmAtoms++;
3448 else if (numFrozenDims > 0)
3450 numPartiallyFrozenVcmAtoms++;
3453 else if (numFrozenDims < DIM)
3459 if (numFullyFrozenVcmAtoms > 0)
3461 std::string warningText = gmx::formatString(
3462 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3463 "removing these atoms from the COMM removal group(s)",
3464 numFullyFrozenVcmAtoms);
3465 warning_note(wi, warningText.c_str());
3467 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3469 std::string warningText = gmx::formatString(
3470 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3471 "removal group(s), due to limitations in the code these still contribute to the "
3472 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3474 numPartiallyFrozenVcmAtoms,
3476 warning(wi, warningText.c_str());
3478 if (numNonVcmAtoms > 0)
3480 std::string warningText = gmx::formatString(
3481 "%d atoms are not part of any center of mass motion removal group.\n"
3482 "This may lead to artifacts.\n"
3483 "In most cases one should use one group for the whole system.",
3485 warning(wi, warningText.c_str());
3489 void do_index(const char* mdparin,
3493 const gmx::MDModulesNotifiers& mdModulesNotifiers,
3497 t_blocka* defaultIndexGroups;
3505 int i, j, k, restnm;
3506 bool bExcl, bTable, bAnneal;
3507 char warn_buf[STRLEN];
3511 fprintf(stderr, "processing index file...\n");
3515 snew(defaultIndexGroups, 1);
3516 snew(defaultIndexGroups->index, 1);
3518 atoms_all = gmx_mtop_global_atoms(*mtop);
3519 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3520 done_atom(&atoms_all);
3524 defaultIndexGroups = init_index(ndx, &gnames);
3527 SimulationGroups* groups = &mtop->groups;
3528 natoms = mtop->natoms;
3529 symtab = &mtop->symtab;
3531 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3533 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3535 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3536 restnm = groups->groupNames.size() - 1;
3537 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3538 srenew(gnames, defaultIndexGroups->nr + 1);
3539 gnames[restnm] = *(groups->groupNames.back());
3541 set_warning_line(wi, mdparin, -1);
3543 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3544 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3545 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3546 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3547 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3550 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3552 temperatureCouplingGroupNames.size(),
3553 temperatureCouplingReferenceValues.size(),
3554 temperatureCouplingTauValues.size());
3557 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3558 do_numbering(natoms,
3560 temperatureCouplingGroupNames,
3563 SimulationAtomGroupType::TemperatureCoupling,
3565 useReferenceTemperature ? GroupCoverage::All : GroupCoverage::AllGenerateRest,
3568 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3570 snew(ir->opts.nrdf, nr);
3571 snew(ir->opts.tau_t, nr);
3572 snew(ir->opts.ref_t, nr);
3573 if (ir->eI == IntegrationAlgorithm::BD && ir->bd_fric == 0)
3575 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3578 if (useReferenceTemperature)
3580 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3582 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3586 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3587 for (i = 0; (i < nr); i++)
3589 if ((ir->eI == IntegrationAlgorithm::BD) && ir->opts.tau_t[i] <= 0)
3592 "With integrator %s tau-t should be larger than 0",
3593 enumValueToString(ir->eI));
3594 warning_error(wi, warn_buf);
3597 if (ir->etc != TemperatureCoupling::VRescale && ir->opts.tau_t[i] == 0)
3601 "tau-t = -1 is the value to signal that a group should not have "
3602 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3605 if (ir->opts.tau_t[i] >= 0)
3607 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3610 if (ir->etc != TemperatureCoupling::No && ir->nsttcouple == -1)
3612 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3617 if ((ir->etc == TemperatureCoupling::NoseHoover) && (ir->epc == PressureCoupling::Berendsen))
3620 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3621 "md-vv; use either vrescale temperature with berendsen pressure or "
3622 "Nose-Hoover temperature with MTTK pressure");
3624 if (ir->epc == PressureCoupling::Mttk)
3626 if (ir->etc != TemperatureCoupling::NoseHoover)
3629 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3634 if (ir->nstpcouple != ir->nsttcouple)
3636 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3637 ir->nstpcouple = ir->nsttcouple = mincouple;
3639 "for current Trotter decomposition methods with vv, nsttcouple and "
3640 "nstpcouple must be equal. Both have been reset to "
3641 "min(nsttcouple,nstpcouple) = %d",
3643 warning_note(wi, warn_buf);
3648 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3649 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3651 if (ETC_ANDERSEN(ir->etc))
3653 if (ir->nsttcouple != 1)
3657 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3658 "need for larger nsttcouple > 1, since no global parameters are computed. "
3659 "nsttcouple has been reset to 1");
3660 warning_note(wi, warn_buf);
3663 nstcmin = tcouple_min_integration_steps(ir->etc);
3666 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3669 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3670 "least %d times larger than nsttcouple*dt (%g)",
3671 enumValueToString(ir->etc),
3674 ir->nsttcouple * ir->delta_t);
3675 warning(wi, warn_buf);
3678 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3679 for (i = 0; (i < nr); i++)
3681 if (ir->opts.ref_t[i] < 0)
3683 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3686 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3687 if we are in this conditional) if mc_temp is negative */
3688 if (ir->expandedvals->mc_temp < 0)
3690 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3694 /* Simulated annealing for each group. There are nr groups */
3695 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3696 if (simulatedAnnealingGroupNames.size() == 1
3697 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3699 simulatedAnnealingGroupNames.resize(0);
3701 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3704 "Wrong number of annealing values: %zu (for %d groups)\n",
3705 simulatedAnnealingGroupNames.size(),
3710 snew(ir->opts.annealing, nr);
3711 snew(ir->opts.anneal_npoints, nr);
3712 snew(ir->opts.anneal_time, nr);
3713 snew(ir->opts.anneal_temp, nr);
3714 for (i = 0; i < nr; i++)
3716 ir->opts.annealing[i] = SimulatedAnnealing::No;
3717 ir->opts.anneal_npoints[i] = 0;
3718 ir->opts.anneal_time[i] = nullptr;
3719 ir->opts.anneal_temp[i] = nullptr;
3721 if (!simulatedAnnealingGroupNames.empty())
3724 for (i = 0; i < nr; i++)
3726 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3728 ir->opts.annealing[i] = SimulatedAnnealing::No;
3730 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3732 ir->opts.annealing[i] = SimulatedAnnealing::Single;
3735 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3737 ir->opts.annealing[i] = SimulatedAnnealing::Periodic;
3743 /* Read the other fields too */
3744 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3745 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3748 "Found %zu annealing-npoints values for %zu groups\n",
3749 simulatedAnnealingPoints.size(),
3750 simulatedAnnealingGroupNames.size());
3752 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3753 size_t numSimulatedAnnealingFields = 0;
3754 for (i = 0; i < nr; i++)
3756 if (ir->opts.anneal_npoints[i] == 1)
3760 "Please specify at least a start and an end point for annealing\n");
3762 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3763 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3764 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3767 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3769 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3772 "Found %zu annealing-time values, wanted %zu\n",
3773 simulatedAnnealingTimes.size(),
3774 numSimulatedAnnealingFields);
3776 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3777 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3780 "Found %zu annealing-temp values, wanted %zu\n",
3781 simulatedAnnealingTemperatures.size(),
3782 numSimulatedAnnealingFields);
3785 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3786 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3787 convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
3789 simulatedAnnealingTemperatures,
3791 allSimulatedAnnealingTemperatures.data());
3792 for (i = 0, k = 0; i < nr; i++)
3794 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3796 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3797 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3800 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3802 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3808 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3811 "Annealing timepoints out of order: t=%f comes after "
3813 ir->opts.anneal_time[i][j],
3814 ir->opts.anneal_time[i][j - 1]);
3817 if (ir->opts.anneal_temp[i][j] < 0)
3820 "Found negative temperature in annealing: %f\n",
3821 ir->opts.anneal_temp[i][j]);
3826 /* Print out some summary information, to make sure we got it right */
3827 for (i = 0; i < nr; i++)
3829 if (ir->opts.annealing[i] != SimulatedAnnealing::No)
3831 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3833 "Simulated annealing for group %s: %s, %d timepoints\n",
3834 *(groups->groupNames[j]),
3835 enumValueToString(ir->opts.annealing[i]),
3836 ir->opts.anneal_npoints[i]);
3837 fprintf(stderr, "Time (ps) Temperature (K)\n");
3838 /* All terms except the last one */
3839 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3843 ir->opts.anneal_time[i][j],
3844 ir->opts.anneal_temp[i][j]);
3847 /* Finally the last one */
3848 j = ir->opts.anneal_npoints[i] - 1;
3849 if (ir->opts.annealing[i] == SimulatedAnnealing::Single)
3853 ir->opts.anneal_time[i][j],
3854 ir->opts.anneal_temp[i][j]);
3860 ir->opts.anneal_time[i][j],
3861 ir->opts.anneal_temp[i][j]);
3862 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3865 "There is a temperature jump when your annealing "
3877 for (int i = 1; i < ir->pull->ngroup; i++)
3879 const int gid = search_string(
3880 inputrecStrings->pullGroupNames[i].c_str(), defaultIndexGroups->nr, gnames);
3881 GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
3882 atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
3885 process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3887 checkPullCoords(ir->pull->group, ir->pull->coord);
3892 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3895 if (ir->eSwapCoords != SwapType::No)
3897 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3900 /* Make indices for IMD session */
3903 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3906 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3907 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3908 mdModulesNotifiers.preProcessingNotifier_.notify(defaultIndexGroupsAndNames);
3910 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3911 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3912 if (freezeDims.size() != DIM * freezeGroupNames.size())
3915 "Invalid Freezing input: %zu groups and %zu freeze values",
3916 freezeGroupNames.size(),
3919 do_numbering(natoms,
3924 SimulationAtomGroupType::Freeze,
3926 GroupCoverage::AllGenerateRest,
3929 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3930 ir->opts.ngfrz = nr;
3931 snew(ir->opts.nFreeze, nr);
3932 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3934 for (j = 0; (j < DIM); j++, k++)
3936 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3937 if (!ir->opts.nFreeze[i][j])
3939 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3942 "Please use Y(ES) or N(O) for freezedim only "
3944 freezeDims[k].c_str());
3945 warning(wi, warn_buf);
3950 for (; (i < nr); i++)
3952 for (j = 0; (j < DIM); j++)
3954 ir->opts.nFreeze[i][j] = 0;
3958 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3959 do_numbering(natoms,
3964 SimulationAtomGroupType::EnergyOutput,
3966 GroupCoverage::AllGenerateRest,
3969 add_wall_energrps(groups, ir->nwall, symtab);
3970 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3971 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3972 do_numbering(natoms,
3977 SimulationAtomGroupType::MassCenterVelocityRemoval,
3979 vcmGroupNames.empty() ? GroupCoverage::AllGenerateRest : GroupCoverage::Partial,
3983 if (ir->comm_mode != ComRemovalAlgorithm::No)
3985 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3988 /* Now we have filled the freeze struct, so we can calculate NRDF */
3989 calc_nrdf(mtop, ir, gnames);
3991 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3992 do_numbering(natoms,
3997 SimulationAtomGroupType::User1,
3999 GroupCoverage::AllGenerateRest,
4002 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
4003 do_numbering(natoms,
4008 SimulationAtomGroupType::User2,
4010 GroupCoverage::AllGenerateRest,
4013 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
4014 do_numbering(natoms,
4016 compressedXGroupNames,
4019 SimulationAtomGroupType::CompressedPositionOutput,
4021 GroupCoverage::OneGroup,
4024 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
4025 do_numbering(natoms,
4027 orirefFitGroupNames,
4030 SimulationAtomGroupType::OrientationRestraintsFit,
4032 GroupCoverage::AllGenerateRest,
4036 /* MiMiC QMMM input processing */
4037 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
4038 if (qmGroupNames.size() > 1)
4040 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
4042 /* group rest, if any, is always MM! */
4043 do_numbering(natoms,
4048 SimulationAtomGroupType::QuantumMechanics,
4050 GroupCoverage::AllGenerateRest,
4053 ir->opts.ngQM = qmGroupNames.size();
4055 /* end of MiMiC QMMM input */
4059 for (auto group : gmx::keysOf(groups->groups))
4061 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
4062 for (const auto& entry : groups->groups[group])
4064 fprintf(stderr, " %s", *(groups->groupNames[entry]));
4066 fprintf(stderr, "\n");
4070 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
4071 snew(ir->opts.egp_flags, nr * nr);
4073 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
4074 if (bExcl && ir->cutoff_scheme == CutoffScheme::Verlet)
4076 warning_error(wi, "Energy group exclusions are currently not supported");
4078 if (bExcl && EEL_FULL(ir->coulombtype))
4080 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
4083 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
4084 if (bTable && !(ir->vdwtype == VanDerWaalsType::User)
4085 && !(ir->coulombtype == CoulombInteractionType::User)
4086 && !(ir->coulombtype == CoulombInteractionType::PmeUser)
4087 && !(ir->coulombtype == CoulombInteractionType::PmeUserSwitch))
4090 "Can only have energy group pair tables in combination with user tables for VdW "
4094 /* final check before going out of scope if simulated tempering variables
4095 * need to be set to default values.
4097 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
4099 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
4102 "the value for nstexpanded was not specified for "
4103 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
4104 "by default, but it is recommended to set it to an explicit value!",
4105 ir->expandedvals->nstexpanded));
4107 for (i = 0; (i < defaultIndexGroups->nr); i++)
4112 done_blocka(defaultIndexGroups);
4113 sfree(defaultIndexGroups);
4117 static void check_disre(const gmx_mtop_t& mtop)
4119 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
4121 const gmx_ffparams_t& ffparams = mtop.ffparams;
4124 for (int i = 0; i < ffparams.numTypes(); i++)
4126 int ftype = ffparams.functype[i];
4127 if (ftype == F_DISRES)
4129 int label = ffparams.iparams[i].disres.label;
4130 if (label == old_label)
4132 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
4141 "Found %d double distance restraint indices,\n"
4142 "probably the parameters for multiple pairs in one restraint "
4143 "are not identical\n",
4149 //! Returns whether dimensions have an absolute reference due to walls, pbc or freezing
4150 static BasicVector<bool> haveAbsoluteReference(const t_inputrec& ir)
4152 BasicVector<bool> absRef = { false, false, false };
4154 /* Check the degrees of freedom of the COM (not taking COMM removal into account) */
4155 for (int d = 0; d < DIM; d++)
4157 absRef[d] = (d >= ndof_com(&ir));
4159 /* Check for freeze groups */
4160 for (int g = 0; g < ir.opts.ngfrz; g++)
4162 for (int d = 0; d < DIM; d++)
4164 if (ir.opts.nFreeze[g][d] != 0)
4174 //! Returns whether position restraints are used for dimensions
4175 static BasicVector<bool> havePositionRestraints(const gmx_mtop_t& sys)
4177 BasicVector<bool> havePosres = { false, false, false };
4179 for (const auto ilists : IListRange(sys))
4181 const auto& posResList = ilists.list()[F_POSRES];
4182 const auto& fbPosResList = ilists.list()[F_FBPOSRES];
4183 if (ilists.nmol() > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
4185 for (int i = 0; i < posResList.size(); i += 2)
4187 const t_iparams& pr = sys.ffparams.iparams[posResList.iatoms[i]];
4188 for (int d = 0; d < DIM; d++)
4190 if (pr.posres.fcA[d] != 0)
4192 havePosres[d] = true;
4196 for (int i = 0; i < fbPosResList.size(); i += 2)
4198 /* Check for flat-bottom posres */
4199 const t_iparams& pr = sys.ffparams.iparams[fbPosResList.iatoms[i]];
4200 if (pr.fbposres.k != 0)
4202 switch (pr.fbposres.geom)
4204 case efbposresSPHERE: havePosres = { true, true, true }; break;
4205 case efbposresCYLINDERX: havePosres[YY] = havePosres[ZZ] = true; break;
4206 case efbposresCYLINDERY: havePosres[XX] = havePosres[ZZ] = true; break;
4207 case efbposresCYLINDER:
4208 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4209 case efbposresCYLINDERZ: havePosres[XX] = havePosres[YY] = true; break;
4210 case efbposresX: /* d=XX */
4211 case efbposresY: /* d=YY */
4212 case efbposresZ: /* d=ZZ */
4213 havePosres[pr.fbposres.geom - efbposresX] = true;
4217 "Invalid geometry for flat-bottom position restraint.\n"
4218 "Expected nr between 1 and %d. Found %d\n",
4230 static void check_combination_rule_differences(const gmx_mtop_t& mtop,
4232 bool* bC6ParametersWorkWithGeometricRules,
4233 bool* bC6ParametersWorkWithLBRules,
4234 bool* bLBRulesPossible)
4236 int ntypes, tpi, tpj;
4239 double c6i, c6j, c12i, c12j;
4240 double c6, c6_geometric, c6_LB;
4241 double sigmai, sigmaj, epsi, epsj;
4242 bool bCanDoLBRules, bCanDoGeometricRules;
4245 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4246 * force-field floating point parameters.
4249 ptr = getenv("GMX_LJCOMB_TOL");
4253 double gmx_unused canary;
4255 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4258 FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4263 *bC6ParametersWorkWithLBRules = TRUE;
4264 *bC6ParametersWorkWithGeometricRules = TRUE;
4265 bCanDoLBRules = TRUE;
4266 ntypes = mtop.ffparams.atnr;
4267 snew(typecount, ntypes);
4268 gmx_mtop_count_atomtypes(mtop, state, typecount);
4269 *bLBRulesPossible = TRUE;
4270 for (tpi = 0; tpi < ntypes; ++tpi)
4272 c6i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4273 c12i = mtop.ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4274 for (tpj = tpi; tpj < ntypes; ++tpj)
4276 c6j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4277 c12j = mtop.ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4278 c6 = mtop.ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4279 c6_geometric = std::sqrt(c6i * c6j);
4280 if (!gmx_numzero(c6_geometric))
4282 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4284 sigmai = gmx::sixthroot(c12i / c6i);
4285 sigmaj = gmx::sixthroot(c12j / c6j);
4286 epsi = c6i * c6i / (4.0 * c12i);
4287 epsj = c6j * c6j / (4.0 * c12j);
4288 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4292 *bLBRulesPossible = FALSE;
4293 c6_LB = c6_geometric;
4295 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4300 *bC6ParametersWorkWithLBRules = FALSE;
4303 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4305 if (!bCanDoGeometricRules)
4307 *bC6ParametersWorkWithGeometricRules = FALSE;
4314 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t& mtop, warninp_t wi)
4316 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4318 check_combination_rule_differences(
4319 mtop, 0, &bC6ParametersWorkWithGeometricRules, &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4320 if (ir->ljpme_combination_rule == LongRangeVdW::LB)
4322 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4325 "You are using arithmetic-geometric combination rules "
4326 "in LJ-PME, but your non-bonded C6 parameters do not "
4327 "follow these rules.");
4332 if (!bC6ParametersWorkWithGeometricRules)
4334 if (ir->eDispCorr != DispersionCorrectionType::No)
4337 "You are using geometric combination rules in "
4338 "LJ-PME, but your non-bonded C6 parameters do "
4339 "not follow these rules. "
4340 "This will introduce very small errors in the forces and energies in "
4341 "your simulations. Dispersion correction will correct total energy "
4342 "and/or pressure for isotropic systems, but not forces or surface "
4348 "You are using geometric combination rules in "
4349 "LJ-PME, but your non-bonded C6 parameters do "
4350 "not follow these rules. "
4351 "This will introduce very small errors in the forces and energies in "
4352 "your simulations. If your system is homogeneous, consider using "
4353 "dispersion correction "
4354 "for the total energy and pressure.");
4360 static bool allTrue(const BasicVector<bool>& boolVector)
4362 return boolVector[0] && boolVector[1] && boolVector[2];
4365 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4367 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4368 GMX_ASSERT(gmx::checkMtsRequirements(*ir).empty(), "All MTS requirements should be met here");
4370 char err_buf[STRLEN];
4373 gmx_mtop_atomloop_block_t aloopb;
4374 char warn_buf[STRLEN];
4376 set_warning_line(wi, mdparin, -1);
4378 if (ir->comm_mode != ComRemovalAlgorithm::No && allTrue(havePositionRestraints(*sys)))
4381 "Removing center of mass motion in the presence of position restraints might "
4382 "cause artifacts. When you are using position restraints to equilibrate a "
4383 "macro-molecule, the artifacts are usually negligible.");
4386 if (ir->cutoff_scheme == CutoffScheme::Verlet && ir->verletbuf_tol > 0 && ir->nstlist > 1
4387 && ((EI_MD(ir->eI) || EI_SD(ir->eI))
4388 && (ir->etc == TemperatureCoupling::VRescale || ir->etc == TemperatureCoupling::Berendsen)))
4390 /* Check if a too small Verlet buffer might potentially
4391 * cause more drift than the thermostat can couple off.
4393 /* Temperature error fraction for warning and suggestion */
4394 const real T_error_warn = 0.002;
4395 const real T_error_suggest = 0.001;
4396 /* For safety: 2 DOF per atom (typical with constraints) */
4397 const real nrdf_at = 2;
4398 real T, tau, max_T_error;
4403 for (i = 0; i < ir->opts.ngtc; i++)
4405 T = std::max(T, ir->opts.ref_t[i]);
4406 tau = std::max(tau, ir->opts.tau_t[i]);
4410 /* This is a worst case estimate of the temperature error,
4411 * assuming perfect buffer estimation and no cancelation
4412 * of errors. The factor 0.5 is because energy distributes
4413 * equally over Ekin and Epot.
4415 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * gmx::c_boltz * T);
4416 if (max_T_error > T_error_warn)
4419 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4420 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4421 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4422 "%.0e or decrease tau_t.",
4427 100 * T_error_suggest,
4428 ir->verletbuf_tol * T_error_suggest / max_T_error);
4429 warning(wi, warn_buf);
4434 if (ETC_ANDERSEN(ir->etc))
4438 for (i = 0; i < ir->opts.ngtc; i++)
4441 "all tau_t must currently be equal using Andersen temperature control, "
4442 "violated for group %d",
4444 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4446 "all tau_t must be positive using Andersen temperature control, "
4450 CHECK(ir->opts.tau_t[i] < 0);
4453 if (ir->etc == TemperatureCoupling::AndersenMassive && ir->comm_mode != ComRemovalAlgorithm::No)
4455 for (i = 0; i < ir->opts.ngtc; i++)
4457 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4459 "tau_t/delta_t for group %d for temperature control method %s must be a "
4460 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4461 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4464 enumValueToString(ir->etc),
4468 CHECK(nsteps % ir->nstcomm != 0);
4473 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != IntegrationAlgorithm::BD
4474 && ir->comm_mode == ComRemovalAlgorithm::No
4475 && !(allTrue(haveAbsoluteReference(*ir)) || allTrue(havePositionRestraints(*sys)) || ir->nsteps <= 10)
4476 && !ETC_ANDERSEN(ir->etc))
4479 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4480 "rounding errors can lead to build up of kinetic energy of the center of mass");
4483 if (ir->epc == PressureCoupling::ParrinelloRahman && ir->etc == TemperatureCoupling::NoseHoover)
4486 for (int g = 0; g < ir->opts.ngtc; g++)
4488 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4490 if (ir->tau_p < 1.9 * tau_t_max)
4492 std::string message = gmx::formatString(
4493 "With %s T-coupling and %s p-coupling, "
4494 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4495 enumValueToString(ir->etc),
4496 enumValueToString(ir->epc),
4501 warning(wi, message.c_str());
4505 /* Check for pressure coupling with absolute position restraints */
4506 if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == RefCoordScaling::No)
4508 const BasicVector<bool> havePosres = havePositionRestraints(*sys);
4510 for (m = 0; m < DIM; m++)
4512 if (havePosres[m] && norm2(ir->compress[m]) > 0)
4515 "You are using pressure coupling with absolute position restraints, "
4516 "this will give artifacts. Use the refcoord_scaling option.");
4524 aloopb = gmx_mtop_atomloop_block_init(*sys);
4526 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4528 if (atom->q != 0 || atom->qB != 0)
4536 if (EEL_FULL(ir->coulombtype))
4539 "You are using full electrostatics treatment %s for a system without charges.\n"
4540 "This costs a lot of performance for just processing zeros, consider using %s "
4542 enumValueToString(ir->coulombtype),
4543 enumValueToString(CoulombInteractionType::Cut));
4544 warning(wi, err_buf);
4549 if (ir->coulombtype == CoulombInteractionType::Cut && ir->rcoulomb > 0)
4552 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4553 "You might want to consider using %s electrostatics.\n",
4554 enumValueToString(CoulombInteractionType::Pme));
4555 warning_note(wi, err_buf);
4559 /* Check if combination rules used in LJ-PME are the same as in the force field */
4560 if (EVDW_PME(ir->vdwtype))
4562 check_combination_rules(ir, *sys, wi);
4565 /* Generalized reaction field */
4566 if (ir->coulombtype == CoulombInteractionType::GRFNotused)
4569 "Generalized reaction-field electrostatics is no longer supported. "
4570 "You can use normal reaction-field instead and compute the reaction-field "
4571 "constant by hand.");
4574 if (ir->efep != FreeEnergyPerturbationType::No && ir->fepvals->sc_alpha != 0
4575 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4577 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4585 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4587 if (ir->pull->coord[i].eGeom != PullGroupGeometry::Transformation
4588 && (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0))
4590 const auto absRef = haveAbsoluteReference(*ir);
4591 const auto havePosres = havePositionRestraints(*sys);
4592 for (m = 0; m < DIM; m++)
4594 if (ir->pull->coord[i].dim[m] && !(absRef[m] || havePosres[m]))
4597 "You are using an absolute reference for pulling, but the rest of "
4598 "the system does not have an absolute reference. This will lead to "
4607 for (i = 0; i < 3; i++)
4609 for (m = 0; m <= i; m++)
4611 if ((ir->epc != PressureCoupling::No && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4613 for (c = 0; c < ir->pull->ncoord; c++)
4615 if (ir->pull->coord[c].eGeom == PullGroupGeometry::DirectionPBC
4616 && ir->pull->coord[c].vec[m] != 0)
4619 "Can not have dynamic box while using pull geometry '%s' "
4621 enumValueToString(ir->pull->coord[c].eGeom),
4633 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4635 char warn_buf[STRLEN];
4638 ptr = check_box(ir->pbcType, box);
4641 warning_error(wi, ptr);
4644 if (bHasNormalConstraints && ir->eConstrAlg == ConstraintAlgorithm::Shake)
4646 if (ir->shake_tol <= 0.0)
4648 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4649 warning_error(wi, warn_buf);
4653 if ((ir->eConstrAlg == ConstraintAlgorithm::Lincs) && bHasNormalConstraints)
4655 /* If we have Lincs constraints: */
4656 if (ir->eI == IntegrationAlgorithm::MD && ir->etc == TemperatureCoupling::No
4657 && ir->eConstrAlg == ConstraintAlgorithm::Lincs && ir->nLincsIter == 1)
4660 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4661 warning_note(wi, warn_buf);
4664 if ((ir->eI == IntegrationAlgorithm::CG || ir->eI == IntegrationAlgorithm::LBFGS)
4665 && (ir->nProjOrder < 8))
4668 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4669 enumValueToString(ir->eI));
4670 warning_note(wi, warn_buf);
4672 if (ir->epc == PressureCoupling::Mttk)
4674 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4678 if (bHasAnyConstraints && ir->epc == PressureCoupling::Mttk)
4680 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4683 if (ir->LincsWarnAngle > 90.0)
4685 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4686 warning(wi, warn_buf);
4687 ir->LincsWarnAngle = 90.0;
4690 if (ir->pbcType != PbcType::No)
4692 if (ir->nstlist == 0)
4695 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4696 "atoms might cause the simulation to crash.");
4698 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4701 "ERROR: The cut-off length is longer than half the shortest box vector or "
4702 "longer than the smallest box diagonal element. Increase the box size or "
4703 "decrease rlist.\n");
4704 warning_error(wi, warn_buf);