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, 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.
50 #include "gromacs/applied_forces/awh/read_params.h"
51 #include "gromacs/fileio/readinp.h"
52 #include "gromacs/fileio/warninp.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/gmxpreprocess/toputil.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/calc_verletbuf.h"
59 #include "gromacs/mdrun/mdmodules.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/mdtypes/multipletimestepping.h"
63 #include "gromacs/mdtypes/pull_params.h"
64 #include "gromacs/options/options.h"
65 #include "gromacs/options/treesupport.h"
66 #include "gromacs/pbcutil/pbc.h"
67 #include "gromacs/selection/indexutil.h"
68 #include "gromacs/topology/block.h"
69 #include "gromacs/topology/ifunc.h"
70 #include "gromacs/topology/index.h"
71 #include "gromacs/topology/mtop_util.h"
72 #include "gromacs/topology/symtab.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/exceptions.h"
76 #include "gromacs/utility/fatalerror.h"
77 #include "gromacs/utility/filestream.h"
78 #include "gromacs/utility/gmxassert.h"
79 #include "gromacs/utility/ikeyvaluetreeerror.h"
80 #include "gromacs/utility/keyvaluetree.h"
81 #include "gromacs/utility/keyvaluetreebuilder.h"
82 #include "gromacs/utility/keyvaluetreemdpwriter.h"
83 #include "gromacs/utility/keyvaluetreetransform.h"
84 #include "gromacs/utility/mdmodulenotification.h"
85 #include "gromacs/utility/smalloc.h"
86 #include "gromacs/utility/strconvert.h"
87 #include "gromacs/utility/stringcompare.h"
88 #include "gromacs/utility/stringutil.h"
89 #include "gromacs/utility/textwriter.h"
94 /* Resource parameters
95 * Do not change any of these until you read the instruction
96 * in readinp.h. Some cpp's do not take spaces after the backslash
97 * (like the c-shell), which will give you a very weird compiler
101 struct gmx_inputrec_strings
103 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], acc[STRLEN], accgrps[STRLEN], freeze[STRLEN],
104 frdim[STRLEN], energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN],
105 x_compressed_groups[STRLEN], couple_moltype[STRLEN], orirefitgrp[STRLEN],
106 egptable[STRLEN], egpexcl[STRLEN], wall_atomtype[STRLEN], wall_density[STRLEN],
107 deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
108 char fep_lambda[efptNR][STRLEN];
109 char lambda_weights[STRLEN];
110 std::vector<std::string> pullGroupNames;
111 std::vector<std::string> rotateGroupNames;
112 char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
115 static gmx_inputrec_strings* inputrecStrings = nullptr;
117 void init_inputrec_strings()
122 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
123 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
125 inputrecStrings = new gmx_inputrec_strings();
128 void done_inputrec_strings()
130 delete inputrecStrings;
131 inputrecStrings = nullptr;
137 egrptpALL, /* All particles have to be a member of a group. */
138 egrptpALL_GENREST, /* A rest group with name is generated for particles *
139 * that are not part of any group. */
140 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
141 * for the rest group. */
142 egrptpONE /* Merge all selected groups into one group, *
143 * make a rest group for the remaining particles. */
146 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
147 "h-angles", "all-angles", nullptr };
149 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
151 static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lambdas)
156 for (i = 0; i < ntemps; i++)
158 /* simple linear scaling -- allows more control */
159 if (simtemp->eSimTempScale == esimtempLINEAR)
161 simtemp->temperatures[i] =
163 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
165 else if (simtemp->eSimTempScale
166 == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
168 simtemp->temperatures[i] = simtemp->simtemp_low
169 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
170 static_cast<real>((1.0 * i) / (ntemps - 1)));
172 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
174 simtemp->temperatures[i] = simtemp->simtemp_low
175 + (simtemp->simtemp_high - simtemp->simtemp_low)
176 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
181 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
182 gmx_fatal(FARGS, "%s", errorstr);
188 static void _low_check(bool b, const char* s, warninp_t wi)
192 warning_error(wi, s);
196 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
200 if (*p > 0 && *p % nst != 0)
202 /* Round up to the next multiple of nst */
203 *p = ((*p) / nst + 1) * nst;
204 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
209 static int lcd(int n1, int n2)
214 for (i = 2; (i <= n1 && i <= n2); i++)
216 if (n1 % i == 0 && n2 % i == 0)
225 //! Convert legacy mdp entries to modern ones.
226 static void process_interaction_modifier(int* eintmod)
228 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
230 *eintmod = eintmodPOTSHIFT;
234 static void checkMtsRequirement(const t_inputrec& ir, const char* param, const int nstValue, warninp_t wi)
236 GMX_RELEASE_ASSERT(ir.mtsLevels.size() >= 2, "Need at least two levels for MTS");
237 const int mtsFactor = ir.mtsLevels.back().stepFactor;
238 if (nstValue % mtsFactor != 0)
240 auto message = gmx::formatString(
241 "With MTS, %s = %d should be a multiple of mts-factor = %d", param, nstValue, mtsFactor);
242 warning_error(wi, message.c_str());
246 static void setupMtsLevels(gmx::ArrayRef<gmx::MtsLevel> mtsLevels,
247 const t_inputrec& ir,
248 const t_gromppopts& opts,
251 if (!(ir.eI == eiMD || ir.eI == eiSD1))
253 auto message = gmx::formatString(
254 "Multiple time stepping is only supported with integrators %s and %s",
255 ei_names[eiMD], ei_names[eiSD1]);
256 warning_error(wi, message.c_str());
258 if (opts.numMtsLevels != 2)
260 warning_error(wi, "Only mts-levels = 2 is supported");
264 const std::vector<std::string> inputForceGroups = gmx::splitString(opts.mtsLevel2Forces);
265 auto& forceGroups = mtsLevels[1].forceGroups;
266 for (const auto& inputForceGroup : inputForceGroups)
270 for (const auto& forceGroupName : gmx::mtsForceGroupNames)
272 if (gmx::equalCaseInsensitive(inputForceGroup, forceGroupName))
274 forceGroups.set(nameIndex);
282 gmx::formatString("Unknown MTS force group '%s'", inputForceGroup.c_str());
283 warning_error(wi, message.c_str());
287 if (mtsLevels[1].stepFactor <= 1)
289 gmx_fatal(FARGS, "mts-factor should be larger than 1");
292 // Make the level 0 use the complement of the force groups of group 1
293 mtsLevels[0].forceGroups = ~mtsLevels[1].forceGroups;
294 mtsLevels[0].stepFactor = 1;
296 if ((EEL_FULL(ir.coulombtype) || EVDW_PME(ir.vdwtype))
297 && !mtsLevels[1].forceGroups[static_cast<int>(gmx::MtsForceGroups::LongrangeNonbonded)])
300 "With long-range electrostatics and/or LJ treatment, the long-range part "
301 "has to be part of the mts-level2-forces");
304 if (ir.nstcalcenergy > 0)
306 checkMtsRequirement(ir, "nstcalcenergy", ir.nstcalcenergy, wi);
308 checkMtsRequirement(ir, "nstenergy", ir.nstenergy, wi);
309 checkMtsRequirement(ir, "nstlog", ir.nstlog, wi);
310 if (ir.efep != efepNO)
312 checkMtsRequirement(ir, "nstdhdl", ir.fepvals->nstdhdl, wi);
317 void check_ir(const char* mdparin,
318 const gmx::MdModulesNotifier& mdModulesNotifier,
322 /* Check internal consistency.
323 * NOTE: index groups are not set here yet, don't check things
324 * like temperature coupling group options here, but in triple_check
327 /* Strange macro: first one fills the err_buf, and then one can check
328 * the condition, which will print the message and increase the error
331 #define CHECK(b) _low_check(b, err_buf, wi)
332 char err_buf[256], warn_buf[STRLEN];
335 t_lambda* fep = ir->fepvals;
336 t_expanded* expand = ir->expandedvals;
338 set_warning_line(wi, mdparin, -1);
340 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
342 sprintf(warn_buf, "%s electrostatics is no longer supported", eel_names[eelRF_NEC_UNSUPPORTED]);
343 warning_error(wi, warn_buf);
346 /* BASIC CUT-OFF STUFF */
347 if (ir->rcoulomb < 0)
349 warning_error(wi, "rcoulomb should be >= 0");
353 warning_error(wi, "rvdw should be >= 0");
355 if (ir->rlist < 0 && !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
357 warning_error(wi, "rlist should be >= 0");
360 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
361 "neighbour-list update scheme for efficient buffering for improved energy "
362 "conservation, please use the Verlet cut-off scheme instead.)");
363 CHECK(ir->nstlist < 0);
365 process_interaction_modifier(&ir->coulomb_modifier);
366 process_interaction_modifier(&ir->vdw_modifier);
368 if (ir->cutoff_scheme == ecutsGROUP)
371 "The group cutoff scheme has been removed since GROMACS 2020. "
372 "Please use the Verlet cutoff scheme.");
374 if (ir->cutoff_scheme == ecutsVERLET)
378 /* Normal Verlet type neighbor-list, currently only limited feature support */
379 if (inputrec2nboundeddim(ir) < 3)
381 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
384 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
385 if (ir->rcoulomb != ir->rvdw)
387 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
388 // for PME load balancing, we can support this exception.
389 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == evdwCUT
390 && ir->rcoulomb > ir->rvdw);
391 if (!bUsesPmeTwinRangeKernel)
394 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
395 "rcoulomb>rvdw with PME electrostatics)");
399 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
401 if (ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT)
403 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
406 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
408 evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
409 warning_note(wi, warn_buf);
411 ir->vdwtype = evdwCUT;
415 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
416 evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
417 warning_error(wi, warn_buf);
421 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
424 "With Verlet lists only cut-off and PME LJ interactions are supported");
426 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) || EEL_PME(ir->coulombtype)
427 || ir->coulombtype == eelEWALD))
430 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
431 "electrostatics are supported");
433 if (!(ir->coulomb_modifier == eintmodNONE || ir->coulomb_modifier == eintmodPOTSHIFT))
435 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
436 warning_error(wi, warn_buf);
439 if (EEL_USER(ir->coulombtype))
441 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme",
442 eel_names[ir->coulombtype]);
443 warning_error(wi, warn_buf);
446 if (ir->nstlist <= 0)
448 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
451 if (ir->nstlist < 10)
454 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
455 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
459 rc_max = std::max(ir->rvdw, ir->rcoulomb);
463 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
464 ir->verletbuf_tol = 0;
467 else if (ir->verletbuf_tol <= 0)
469 if (ir->verletbuf_tol == 0)
471 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
474 if (ir->rlist < rc_max)
477 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
480 if (ir->rlist == rc_max && ir->nstlist > 1)
484 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
485 "buffer. The cluster pair list does have a buffering effect, but choosing "
486 "a larger rlist might be necessary for good energy conservation.");
491 if (ir->rlist > rc_max)
494 "You have set rlist larger than the interaction cut-off, but you also "
495 "have verlet-buffer-tolerance > 0. Will set rlist using "
496 "verlet-buffer-tolerance.");
499 if (ir->nstlist == 1)
501 /* No buffer required */
506 if (EI_DYNAMICS(ir->eI))
508 if (inputrec2nboundeddim(ir) < 3)
511 "The box volume is required for calculating rlist from the "
512 "energy drift with verlet-buffer-tolerance > 0. You are "
513 "using at least one unbounded dimension, so no volume can be "
514 "computed. Either use a finite box, or set rlist yourself "
515 "together with verlet-buffer-tolerance = -1.");
517 /* Set rlist temporarily so we can continue processing */
522 /* Set the buffer to 5% of the cut-off */
523 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
529 /* GENERAL INTEGRATOR STUFF */
532 if (ir->etc != etcNO)
534 if (EI_RANDOM(ir->eI))
537 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
538 "implicitly. See the documentation for more information on which "
539 "parameters affect temperature for %s.",
540 etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
545 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
547 etcoupl_names[ir->etc], ei_names[ir->eI]);
549 warning_note(wi, warn_buf);
553 if (ir->eI == eiVVAK)
556 "Integrator method %s is implemented primarily for validation purposes; for "
557 "molecular dynamics, you should probably be using %s or %s",
558 ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
559 warning_note(wi, warn_buf);
561 if (!EI_DYNAMICS(ir->eI))
563 if (ir->epc != epcNO)
566 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
567 epcoupl_names[ir->epc], ei_names[ir->eI]);
568 warning_note(wi, warn_buf);
572 if (EI_DYNAMICS(ir->eI))
574 if (ir->nstcalcenergy < 0)
576 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
577 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
579 /* nstcalcenergy larger than nstener does not make sense.
580 * We ideally want nstcalcenergy=nstener.
584 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
588 ir->nstcalcenergy = ir->nstenergy;
592 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
593 || (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
594 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
597 const char* nsten = "nstenergy";
598 const char* nstdh = "nstdhdl";
599 const char* min_name = nsten;
600 int min_nst = ir->nstenergy;
602 /* find the smallest of ( nstenergy, nstdhdl ) */
603 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
604 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
606 min_nst = ir->fepvals->nstdhdl;
609 /* If the user sets nstenergy small, we should respect that */
610 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy,
612 warning_note(wi, warn_buf);
613 ir->nstcalcenergy = min_nst;
616 if (ir->epc != epcNO)
618 if (ir->nstpcouple < 0)
620 ir->nstpcouple = ir_optimal_nstpcouple(ir);
622 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
625 "With multiple time stepping, nstpcouple should be a mutiple of "
630 if (ir->nstcalcenergy > 0)
632 if (ir->efep != efepNO)
634 /* nstdhdl should be a multiple of nstcalcenergy */
635 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
639 /* nstexpanded should be a multiple of nstcalcenergy */
640 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded",
641 &ir->expandedvals->nstexpanded, wi);
643 /* for storing exact averages nstenergy should be
644 * a multiple of nstcalcenergy
646 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
649 // Inquire all MdModules, if their parameters match with the energy
650 // calculation frequency
651 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
652 mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
654 // Emit all errors from the energy calculation frequency checks
655 for (const std::string& energyFrequencyErrorMessage :
656 energyCalculationFrequencyErrors.errorMessages())
658 warning_error(wi, energyFrequencyErrorMessage);
662 if (ir->nsteps == 0 && !ir->bContinuation)
665 "For a correct single-point energy evaluation with nsteps = 0, use "
666 "continuation = yes to avoid constraining the input coordinates.");
670 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
673 "You are doing a continuation with SD or BD, make sure that ld_seed is "
674 "different from the previous run (using ld_seed=-1 will ensure this)");
680 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
681 CHECK(ir->pbcType != PbcType::Xyz);
682 sprintf(err_buf, "with TPI nstlist should be larger than zero");
683 CHECK(ir->nstlist <= 0);
684 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
685 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
689 if ((opts->nshake > 0) && (opts->bMorse))
691 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
692 warning(wi, warn_buf);
695 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
698 "You are doing a continuation with SD or BD, make sure that ld_seed is "
699 "different from the previous run (using ld_seed=-1 will ensure this)");
701 /* verify simulated tempering options */
705 bool bAllTempZero = TRUE;
706 for (i = 0; i < fep->n_lambda; i++)
708 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
709 efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
710 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
711 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
713 bAllTempZero = FALSE;
716 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
717 CHECK(bAllTempZero == TRUE);
719 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
720 CHECK(ir->eI != eiVV);
722 /* check compatability of the temperature coupling with simulated tempering */
724 if (ir->etc == etcNOSEHOOVER)
727 "Nose-Hoover based temperature control such as [%s] my not be "
728 "entirelyconsistent with simulated tempering",
729 etcoupl_names[ir->etc]);
730 warning_note(wi, warn_buf);
733 /* check that the temperatures make sense */
736 "Higher simulated tempering temperature (%g) must be >= than the simulated "
737 "tempering lower temperature (%g)",
738 ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
739 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
741 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero",
742 ir->simtempvals->simtemp_high);
743 CHECK(ir->simtempvals->simtemp_high <= 0);
745 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero",
746 ir->simtempvals->simtemp_low);
747 CHECK(ir->simtempvals->simtemp_low <= 0);
750 /* verify free energy options */
752 if (ir->efep != efepNO)
755 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
756 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
759 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
761 static_cast<int>(fep->sc_r_power));
762 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
765 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
768 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
770 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
772 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
774 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
775 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
777 sprintf(err_buf, "Free-energy not implemented for Ewald");
778 CHECK(ir->coulombtype == eelEWALD);
780 /* check validty of lambda inputs */
781 if (fep->n_lambda == 0)
783 /* Clear output in case of no states:*/
784 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.",
785 fep->init_fep_state);
786 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
790 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d",
791 fep->init_fep_state, fep->n_lambda - 1);
792 CHECK((fep->init_fep_state >= fep->n_lambda));
796 "Lambda state must be set, either with init-lambda-state or with init-lambda");
797 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
800 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
801 "init-lambda-state or with init-lambda, but not both",
802 fep->init_lambda, fep->init_fep_state);
803 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
806 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
810 for (i = 0; i < efptNR; i++)
812 if (fep->separate_dvdl[i])
817 if (n_lambda_terms > 1)
820 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
821 "use init-lambda to set lambda state (except for slow growth). Use "
822 "init-lambda-state instead.");
823 warning(wi, warn_buf);
826 if (n_lambda_terms < 2 && fep->n_lambda > 0)
829 "init-lambda is deprecated for setting lambda state (except for slow "
830 "growth). Use init-lambda-state instead.");
834 for (j = 0; j < efptNR; j++)
836 for (i = 0; i < fep->n_lambda; i++)
838 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
839 efpt_names[j], fep->all_lambda[j][i]);
840 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
844 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
846 for (i = 0; i < fep->n_lambda; i++)
849 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
850 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
851 "crashes, and is not supported.",
852 i, fep->all_lambda[efptVDW][i], fep->all_lambda[efptCOUL][i]);
853 CHECK((fep->sc_alpha > 0)
854 && (((fep->all_lambda[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL][i] < 1.0))
855 && ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0))));
859 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
861 real sigma, lambda, r_sc;
864 /* Maximum estimate for A and B charges equal with lambda power 1 */
866 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
867 1.0 / fep->sc_r_power);
869 "With PME there is a minor soft core effect present at the cut-off, "
870 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
871 "energy conservation, but usually other effects dominate. With a common sigma "
872 "value of %g nm the fraction of the particle-particle potential at the cut-off "
873 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
874 fep->sc_r_power, sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
875 warning_note(wi, warn_buf);
878 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
879 be treated differently, but that's the next step */
881 for (i = 0; i < efptNR; i++)
883 for (j = 0; j < fep->n_lambda; j++)
885 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
886 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
891 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
895 /* checking equilibration of weights inputs for validity */
898 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
900 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
901 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
904 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
906 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
907 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
910 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
911 expand->equil_steps, elmceq_names[elmceqSTEPS]);
912 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
915 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
916 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
917 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
920 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
921 expand->equil_ratio, elmceq_names[elmceqRATIO]);
922 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
925 "weight-equil-number-all-lambda (%d) must be a positive integer if "
926 "lmc-weights-equil=%s",
927 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
928 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
931 "weight-equil-number-samples (%d) must be a positive integer if "
932 "lmc-weights-equil=%s",
933 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
934 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
937 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
938 expand->equil_steps, elmceq_names[elmceqSTEPS]);
939 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
941 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
942 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
943 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
945 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
946 expand->equil_ratio, elmceq_names[elmceqRATIO]);
947 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
949 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
950 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
951 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
953 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
954 CHECK((expand->lmc_repeats <= 0));
955 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
956 CHECK((expand->minvarmin <= 0));
957 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
958 CHECK((expand->c_range < 0));
960 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
962 fep->init_fep_state, expand->lmc_forced_nstart);
963 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
964 && (expand->elmcmove != elmcmoveNO));
965 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
966 CHECK((expand->lmc_forced_nstart < 0));
967 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
968 fep->init_fep_state);
969 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
971 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
972 CHECK((expand->init_wl_delta < 0));
973 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
974 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
975 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
976 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
978 /* if there is no temperature control, we need to specify an MC temperature */
979 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO)
980 && (expand->mc_temp <= 0.0))
983 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
984 "to a positive number");
985 warning_error(wi, err_buf);
987 if (expand->nstTij > 0)
989 sprintf(err_buf, "nstlog must be non-zero");
990 CHECK(ir->nstlog == 0);
991 // Avoid modulus by zero in the case that already triggered an error exit.
995 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
996 expand->nstTij, ir->nstlog);
997 CHECK((expand->nstTij % ir->nstlog) != 0);
1003 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
1004 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
1007 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
1009 if (ir->pbcType == PbcType::No)
1011 if (ir->epc != epcNO)
1013 warning(wi, "Turning off pressure coupling for vacuum system");
1019 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
1020 c_pbcTypeNames[ir->pbcType].c_str());
1021 CHECK(ir->epc != epcNO);
1023 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
1024 CHECK(EEL_FULL(ir->coulombtype));
1026 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
1027 c_pbcTypeNames[ir->pbcType].c_str());
1028 CHECK(ir->eDispCorr != edispcNO);
1031 if (ir->rlist == 0.0)
1034 "can only have neighborlist cut-off zero (=infinite)\n"
1035 "with coulombtype = %s or coulombtype = %s\n"
1036 "without periodic boundary conditions (pbc = %s) and\n"
1037 "rcoulomb and rvdw set to zero",
1038 eel_names[eelCUT], eel_names[eelUSER], c_pbcTypeNames[PbcType::No].c_str());
1039 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
1040 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1042 if (ir->nstlist > 0)
1045 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1046 "nstype=simple and only one MPI rank");
1051 if (ir->nstcomm == 0)
1053 // TODO Change this behaviour. There should be exactly one way
1054 // to turn off an algorithm.
1055 ir->comm_mode = ecmNO;
1057 if (ir->comm_mode != ecmNO)
1059 if (ir->nstcomm < 0)
1061 // TODO Such input was once valid. Now that we've been
1062 // helpful for a few years, we should reject such input,
1063 // lest we have to support every historical decision
1066 "If you want to remove the rotation around the center of mass, you should set "
1067 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1068 "its absolute value");
1069 ir->nstcomm = abs(ir->nstcomm);
1072 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
1075 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
1076 "nstcomm to nstcalcenergy");
1077 ir->nstcomm = ir->nstcalcenergy;
1080 if (ir->comm_mode == ecmANGULAR)
1083 "Can not remove the rotation around the center of mass with periodic "
1085 CHECK(ir->bPeriodicMols);
1086 if (ir->pbcType != PbcType::No)
1089 "Removing the rotation around the center of mass in a periodic system, "
1090 "this can lead to artifacts. Only use this on a single (cluster of) "
1091 "molecules. This cluster should not cross periodic boundaries.");
1096 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
1099 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1100 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1103 warning_note(wi, warn_buf);
1106 /* TEMPERATURE COUPLING */
1107 if (ir->etc == etcYES)
1109 ir->etc = etcBERENDSEN;
1111 "Old option for temperature coupling given: "
1112 "changing \"yes\" to \"Berendsen\"\n");
1115 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
1117 if (ir->opts.nhchainlength < 1)
1120 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1122 ir->opts.nhchainlength);
1123 ir->opts.nhchainlength = 1;
1124 warning(wi, warn_buf);
1127 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1131 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1132 ir->opts.nhchainlength = 1;
1137 ir->opts.nhchainlength = 0;
1140 if (ir->eI == eiVVAK)
1143 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1146 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1149 if (ETC_ANDERSEN(ir->etc))
1151 sprintf(err_buf, "%s temperature control not supported for integrator %s.",
1152 etcoupl_names[ir->etc], ei_names[ir->eI]);
1153 CHECK(!(EI_VV(ir->eI)));
1155 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
1158 "Center of mass removal not necessary for %s. All velocities of coupled "
1159 "groups are rerandomized periodically, so flying ice cube errors will not "
1161 etcoupl_names[ir->etc]);
1162 warning_note(wi, warn_buf);
1166 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1167 "randomized every time step",
1168 ir->nstcomm, etcoupl_names[ir->etc]);
1169 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
1172 if (ir->etc == etcBERENDSEN)
1175 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1176 "might want to consider using the %s thermostat.",
1177 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
1178 warning_note(wi, warn_buf);
1181 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc)) && ir->epc == epcBERENDSEN)
1184 "Using Berendsen pressure coupling invalidates the "
1185 "true ensemble for the thermostat");
1186 warning(wi, warn_buf);
1189 /* PRESSURE COUPLING */
1190 if (ir->epc == epcISOTROPIC)
1192 ir->epc = epcBERENDSEN;
1194 "Old option for pressure coupling given: "
1195 "changing \"Isotropic\" to \"Berendsen\"\n");
1198 if (ir->epc != epcNO)
1200 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1202 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1203 CHECK(ir->tau_p <= 0);
1205 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1208 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1209 "times larger than nstpcouple*dt (%g)",
1210 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1211 warning(wi, warn_buf);
1215 "compressibility must be > 0 when using pressure"
1217 EPCOUPLTYPE(ir->epc));
1218 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1219 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1220 && ir->compress[ZZ][YY] <= 0));
1222 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1225 "You are generating velocities so I am assuming you "
1226 "are equilibrating a system. You are using "
1227 "%s pressure coupling, but this can be "
1228 "unstable for equilibration. If your system crashes, try "
1229 "equilibrating first with Berendsen pressure coupling. If "
1230 "you are not equilibrating the system, you can probably "
1231 "ignore this warning.",
1232 epcoupl_names[ir->epc]);
1233 warning(wi, warn_buf);
1239 if (ir->epc == epcMTTK)
1241 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1245 /* ELECTROSTATICS */
1246 /* More checks are in triple check (grompp.c) */
1248 if (ir->coulombtype == eelSWITCH)
1251 "coulombtype = %s is only for testing purposes and can lead to serious "
1252 "artifacts, advice: use coulombtype = %s",
1253 eel_names[ir->coulombtype], eel_names[eelRF_ZERO]);
1254 warning(wi, warn_buf);
1257 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1260 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1261 "format and exchanging epsilon-r and epsilon-rf",
1263 warning(wi, warn_buf);
1264 ir->epsilon_rf = ir->epsilon_r;
1265 ir->epsilon_r = 1.0;
1268 if (ir->epsilon_r == 0)
1271 "It is pointless to use long-range electrostatics with infinite relative "
1273 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1275 CHECK(EEL_FULL(ir->coulombtype));
1278 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1280 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1281 CHECK(ir->epsilon_r < 0);
1284 if (EEL_RF(ir->coulombtype))
1286 /* reaction field (at the cut-off) */
1288 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1291 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1292 eel_names[ir->coulombtype]);
1293 warning(wi, warn_buf);
1294 ir->epsilon_rf = 0.0;
1297 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1298 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1299 if (ir->epsilon_rf == ir->epsilon_r)
1301 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1302 eel_names[ir->coulombtype]);
1303 warning(wi, warn_buf);
1306 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1307 * means the interaction is zero outside rcoulomb, but it helps to
1308 * provide accurate energy conservation.
1310 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1312 if (ir_coulomb_switched(ir))
1315 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1316 "potential modifier options!",
1317 eel_names[ir->coulombtype]);
1318 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1322 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1325 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1326 "secondary coulomb-modifier.");
1327 CHECK(ir->coulomb_modifier != eintmodNONE);
1329 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1332 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1333 "secondary vdw-modifier.");
1334 CHECK(ir->vdw_modifier != eintmodNONE);
1337 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->vdwtype == evdwSWITCH
1338 || ir->vdwtype == evdwSHIFT)
1341 "The switch/shift interaction settings are just for compatibility; you will get "
1343 "performance from applying potential modifiers to your interactions!\n");
1344 warning_note(wi, warn_buf);
1347 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1349 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1351 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1353 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1354 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1355 "will be good regardless, since ewald_rtol = %g.",
1356 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1357 warning(wi, warn_buf);
1361 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1363 if (ir->rvdw_switch == 0)
1366 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1367 "potential. This suggests it was not set in the mdp, which can lead to large "
1368 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1369 "switching range.");
1370 warning(wi, warn_buf);
1374 if (EEL_FULL(ir->coulombtype))
1376 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER
1377 || ir->coulombtype == eelPMEUSERSWITCH)
1379 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1380 eel_names[ir->coulombtype]);
1381 CHECK(ir->rcoulomb > ir->rlist);
1385 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1387 // TODO: Move these checks into the ewald module with the options class
1389 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1391 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1393 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d",
1394 eel_names[ir->coulombtype], orderMin, orderMax);
1395 warning_error(wi, warn_buf);
1399 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1401 if (ir->ewald_geometry == eewg3D)
1403 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1404 c_pbcTypeNames[ir->pbcType].c_str(), eewg_names[eewg3DC]);
1405 warning(wi, warn_buf);
1407 /* This check avoids extra pbc coding for exclusion corrections */
1408 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1409 CHECK(ir->wall_ewald_zfac < 2);
1411 if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
1413 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1414 eel_names[ir->coulombtype], eewg_names[eewg3DC], c_pbcTypeNames[PbcType::XY].c_str());
1415 warning(wi, warn_buf);
1417 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1419 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1420 CHECK(ir->bPeriodicMols);
1421 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1422 warning_note(wi, warn_buf);
1424 "With epsilon_surface > 0 you can only use domain decomposition "
1425 "when there are only small molecules with all bonds constrained (mdrun will check "
1427 warning_note(wi, warn_buf);
1430 if (ir_vdw_switched(ir))
1432 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1433 CHECK(ir->rvdw_switch >= ir->rvdw);
1435 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1438 "You are applying a switch function to vdw forces or potentials from %g to %g "
1439 "nm, which is more than half the interaction range, whereas switch functions "
1440 "are intended to act only close to the cut-off.",
1441 ir->rvdw_switch, ir->rvdw);
1442 warning_note(wi, warn_buf);
1446 if (ir->vdwtype == evdwPME)
1448 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1450 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1451 evdw_names[ir->vdwtype], eintmod_names[eintmodPOTSHIFT], eintmod_names[eintmodNONE]);
1452 warning_error(wi, err_buf);
1456 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1459 "You have selected user tables with dispersion correction, the dispersion "
1460 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1461 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1462 "really want dispersion correction to -C6/r^6.");
1465 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT) && ir->rvdw != 0)
1467 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1470 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1472 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1475 /* IMPLICIT SOLVENT */
1476 if (ir->coulombtype == eelGB_NOTUSED)
1478 sprintf(warn_buf, "Invalid option %s for coulombtype", eel_names[ir->coulombtype]);
1479 warning_error(wi, warn_buf);
1484 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1489 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1493 /* interpret a number of doubles from a string and put them in an array,
1494 after allocating space for them.
1495 str = the input string
1496 n = the (pre-allocated) number of doubles read
1497 r = the output array of doubles. */
1498 static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
1500 auto values = gmx::splitString(str);
1504 for (int i = 0; i < *n; i++)
1508 (*r)[i] = gmx::fromString<real>(values[i]);
1510 catch (gmx::GromacsException&)
1512 warning_error(wi, "Invalid value " + values[i]
1513 + " in string in mdp file. Expected a real number.");
1519 static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1522 int i, j, max_n_lambda, nweights, nfep[efptNR];
1523 t_lambda* fep = ir->fepvals;
1524 t_expanded* expand = ir->expandedvals;
1525 real** count_fep_lambdas;
1526 bool bOneLambda = TRUE;
1528 snew(count_fep_lambdas, efptNR);
1530 /* FEP input processing */
1531 /* first, identify the number of lambda values for each type.
1532 All that are nonzero must have the same number */
1534 for (i = 0; i < efptNR; i++)
1536 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1539 /* now, determine the number of components. All must be either zero, or equal. */
1542 for (i = 0; i < efptNR; i++)
1544 if (nfep[i] > max_n_lambda)
1546 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1547 must have the same number if its not zero.*/
1552 for (i = 0; i < efptNR; i++)
1556 ir->fepvals->separate_dvdl[i] = FALSE;
1558 else if (nfep[i] == max_n_lambda)
1560 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute
1561 the derivative with respect to the temperature currently */
1563 ir->fepvals->separate_dvdl[i] = TRUE;
1569 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1571 nfep[i], efpt_names[i], max_n_lambda);
1574 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1575 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1577 /* the number of lambdas is the number we've read in, which is either zero
1578 or the same for all */
1579 fep->n_lambda = max_n_lambda;
1581 /* allocate space for the array of lambda values */
1582 snew(fep->all_lambda, efptNR);
1583 /* if init_lambda is defined, we need to set lambda */
1584 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1586 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1588 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1589 for (i = 0; i < efptNR; i++)
1591 snew(fep->all_lambda[i], fep->n_lambda);
1592 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1595 for (j = 0; j < fep->n_lambda; j++)
1597 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1599 sfree(count_fep_lambdas[i]);
1602 sfree(count_fep_lambdas);
1604 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1605 internal bookkeeping -- for now, init_lambda */
1607 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1609 for (i = 0; i < fep->n_lambda; i++)
1611 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1615 /* check to see if only a single component lambda is defined, and soft core is defined.
1616 In this case, turn on coulomb soft core */
1618 if (max_n_lambda == 0)
1624 for (i = 0; i < efptNR; i++)
1626 if ((nfep[i] != 0) && (i != efptFEP))
1632 if ((bOneLambda) && (fep->sc_alpha > 0))
1634 fep->bScCoul = TRUE;
1637 /* Fill in the others with the efptFEP if they are not explicitly
1638 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1639 they are all zero. */
1641 for (i = 0; i < efptNR; i++)
1643 if ((nfep[i] == 0) && (i != efptFEP))
1645 for (j = 0; j < fep->n_lambda; j++)
1647 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1653 /* now read in the weights */
1654 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1657 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1659 else if (nweights != fep->n_lambda)
1661 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1662 nweights, fep->n_lambda);
1664 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1666 expand->nstexpanded = fep->nstdhdl;
1667 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1672 static void do_simtemp_params(t_inputrec* ir)
1675 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1676 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1679 template<typename T>
1680 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1683 for (const auto& input : inputs)
1687 outputs[i] = gmx::fromStdString<T>(input);
1689 catch (gmx::GromacsException&)
1691 auto message = gmx::formatString(
1692 "Invalid value for mdp option %s. %s should only consist of integers separated "
1695 warning_error(wi, message);
1701 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1704 for (const auto& input : inputs)
1708 outputs[i] = gmx::fromString<real>(input);
1710 catch (gmx::GromacsException&)
1712 auto message = gmx::formatString(
1713 "Invalid value for mdp option %s. %s should only consist of real numbers "
1714 "separated by spaces.",
1716 warning_error(wi, message);
1722 static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
1725 for (const auto& input : inputs)
1729 outputs[i][d] = gmx::fromString<real>(input);
1731 catch (gmx::GromacsException&)
1733 auto message = gmx::formatString(
1734 "Invalid value for mdp option %s. %s should only consist of real numbers "
1735 "separated by spaces.",
1737 warning_error(wi, message);
1748 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1750 opts->wall_atomtype[0] = nullptr;
1751 opts->wall_atomtype[1] = nullptr;
1753 ir->wall_atomtype[0] = -1;
1754 ir->wall_atomtype[1] = -1;
1755 ir->wall_density[0] = 0;
1756 ir->wall_density[1] = 0;
1760 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1761 if (wallAtomTypes.size() != size_t(ir->nwall))
1763 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu", ir->nwall,
1764 wallAtomTypes.size());
1766 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1767 for (int i = 0; i < ir->nwall; i++)
1769 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1772 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1774 auto wallDensity = gmx::splitString(wall_density);
1775 if (wallDensity.size() != size_t(ir->nwall))
1777 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall,
1778 wallDensity.size());
1780 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1781 for (int i = 0; i < ir->nwall; i++)
1783 if (ir->wall_density[i] <= 0)
1785 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1792 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1796 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1797 for (int i = 0; i < nwall; i++)
1799 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1800 grps->emplace_back(groups->groupNames.size() - 1);
1805 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1807 /* read expanded ensemble parameters */
1808 printStringNewline(inp, "expanded ensemble variables");
1809 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1810 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1811 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1812 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1813 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1814 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1815 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1816 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1817 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1818 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1819 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1820 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1821 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1822 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1823 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1824 expand->bSymmetrizedTMatrix =
1825 (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1826 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1827 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1828 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1829 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1830 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1831 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1832 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1835 /*! \brief Return whether an end state with the given coupling-lambda
1836 * value describes fully-interacting VDW.
1838 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1839 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1841 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1843 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1849 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1852 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1854 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1856 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1859 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1860 std::string message = gmx::formatExceptionMessageToString(*ex);
1861 warning_error(wi_, message.c_str());
1866 std::string getOptionName(const gmx::KeyValueTreePath& context)
1868 if (mapping_ != nullptr)
1870 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1871 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1874 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1879 const gmx::IKeyValueTreeBackMapping* mapping_;
1884 void get_ir(const char* mdparin,
1885 const char* mdparout,
1886 gmx::MDModules* mdModules,
1889 WriteMdpHeader writeMdpHeader,
1893 double dumdub[2][6];
1895 char warn_buf[STRLEN];
1896 t_lambda* fep = ir->fepvals;
1897 t_expanded* expand = ir->expandedvals;
1899 const char* no_names[] = { "no", nullptr };
1901 init_inputrec_strings();
1902 gmx::TextInputFile stream(mdparin);
1903 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1905 snew(dumstr[0], STRLEN);
1906 snew(dumstr[1], STRLEN);
1908 /* ignore the following deprecated commands */
1909 replace_inp_entry(inp, "title", nullptr);
1910 replace_inp_entry(inp, "cpp", nullptr);
1911 replace_inp_entry(inp, "domain-decomposition", nullptr);
1912 replace_inp_entry(inp, "andersen-seed", nullptr);
1913 replace_inp_entry(inp, "dihre", nullptr);
1914 replace_inp_entry(inp, "dihre-fc", nullptr);
1915 replace_inp_entry(inp, "dihre-tau", nullptr);
1916 replace_inp_entry(inp, "nstdihreout", nullptr);
1917 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1918 replace_inp_entry(inp, "optimize-fft", nullptr);
1919 replace_inp_entry(inp, "adress_type", nullptr);
1920 replace_inp_entry(inp, "adress_const_wf", nullptr);
1921 replace_inp_entry(inp, "adress_ex_width", nullptr);
1922 replace_inp_entry(inp, "adress_hy_width", nullptr);
1923 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1924 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1925 replace_inp_entry(inp, "adress_site", nullptr);
1926 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1927 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1928 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1929 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1930 replace_inp_entry(inp, "rlistlong", nullptr);
1931 replace_inp_entry(inp, "nstcalclr", nullptr);
1932 replace_inp_entry(inp, "pull-print-com2", nullptr);
1933 replace_inp_entry(inp, "gb-algorithm", nullptr);
1934 replace_inp_entry(inp, "nstgbradii", nullptr);
1935 replace_inp_entry(inp, "rgbradii", nullptr);
1936 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1937 replace_inp_entry(inp, "gb-saltconc", nullptr);
1938 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1939 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1940 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1941 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1942 replace_inp_entry(inp, "sa-algorithm", nullptr);
1943 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1944 replace_inp_entry(inp, "ns-type", nullptr);
1946 /* replace the following commands with the clearer new versions*/
1947 replace_inp_entry(inp, "unconstrained-start", "continuation");
1948 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1949 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1950 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1951 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1952 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1953 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1955 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1956 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1957 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1958 setStringEntry(&inp, "include", opts->include, nullptr);
1959 printStringNoNewline(
1960 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1961 setStringEntry(&inp, "define", opts->define, nullptr);
1963 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1964 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1965 printStringNoNewline(&inp, "Start time and timestep in ps");
1966 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1967 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1968 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1969 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1970 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1971 printStringNoNewline(
1972 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1973 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1974 printStringNoNewline(&inp, "Multiple time-stepping");
1975 ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
1978 opts->numMtsLevels = get_eint(&inp, "mts-levels", 2, wi);
1979 ir->mtsLevels.resize(2);
1980 gmx::MtsLevel& mtsLevel = ir->mtsLevels[1];
1981 opts->mtsLevel2Forces = setStringEntry(&inp, "mts-level2-forces",
1982 "longrange-nonbonded nonbonded pair dihedral");
1983 mtsLevel.stepFactor = get_eint(&inp, "mts-level2-factor", 2, wi);
1985 // We clear after reading without dynamics to not force the user to remove MTS mdp options
1986 if (!EI_DYNAMICS(ir->eI))
1989 ir->mtsLevels.clear();
1992 printStringNoNewline(&inp, "mode for center of mass motion removal");
1993 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
1994 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
1995 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
1996 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
1997 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
1999 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2000 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2001 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2002 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2005 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2006 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2007 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2008 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2009 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2010 ir->niter = get_eint(&inp, "niter", 20, wi);
2011 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2012 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2013 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2014 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2015 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2017 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2018 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2020 /* Output options */
2021 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2022 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2023 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2024 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2025 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2026 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2027 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2028 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2029 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2030 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2031 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2032 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2033 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2034 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2035 printStringNoNewline(&inp, "default, all atoms will be written.");
2036 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2037 printStringNoNewline(&inp, "Selection of energy groups");
2038 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2040 /* Neighbor searching */
2041 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2042 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2043 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
2044 printStringNoNewline(&inp, "nblist update frequency");
2045 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2046 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2047 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2048 std::vector<const char*> pbcTypesNamesChar;
2049 for (const auto& pbcTypeName : c_pbcTypeNames)
2051 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2053 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2054 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
2055 printStringNoNewline(&inp,
2056 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2057 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2058 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2059 printStringNoNewline(&inp, "nblist cut-off");
2060 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2061 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2063 /* Electrostatics */
2064 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2065 printStringNoNewline(&inp, "Method for doing electrostatics");
2066 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
2067 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
2068 printStringNoNewline(&inp, "cut-off lengths");
2069 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2070 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2071 printStringNoNewline(&inp,
2072 "Relative dielectric constant for the medium and the reaction field");
2073 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2074 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2075 printStringNoNewline(&inp, "Method for doing Van der Waals");
2076 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
2077 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
2078 printStringNoNewline(&inp, "cut-off lengths");
2079 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2080 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2081 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2082 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
2083 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2084 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2085 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2086 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2087 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2088 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2089 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2090 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2091 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2092 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2093 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2094 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2095 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2096 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2097 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
2098 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
2099 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2101 /* Implicit solvation is no longer supported, but we need grompp
2102 to be able to refuse old .mdp files that would have built a tpr
2103 to run it. Thus, only "no" is accepted. */
2104 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2106 /* Coupling stuff */
2107 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2108 printStringNoNewline(&inp, "Temperature coupling");
2109 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
2110 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2111 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2112 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
2113 printStringNoNewline(&inp, "Groups to couple separately");
2114 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2115 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2116 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2117 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2118 printStringNoNewline(&inp, "pressure coupling");
2119 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
2120 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
2121 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2122 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2123 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2124 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2125 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2126 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2127 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
2130 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2131 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
2132 printStringNoNewline(&inp, "Groups treated with MiMiC");
2133 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2135 /* Simulated annealing */
2136 printStringNewline(&inp, "SIMULATED ANNEALING");
2137 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2138 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2139 printStringNoNewline(&inp,
2140 "Number of time points to use for specifying annealing in each group");
2141 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2142 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2143 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2144 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2145 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2148 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2149 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
2150 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2151 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2154 printStringNewline(&inp, "OPTIONS FOR BONDS");
2155 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2156 printStringNoNewline(&inp, "Type of constraint algorithm");
2157 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
2158 printStringNoNewline(&inp, "Do not constrain the start configuration");
2159 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
2160 printStringNoNewline(&inp,
2161 "Use successive overrelaxation to reduce the number of shake iterations");
2162 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
2163 printStringNoNewline(&inp, "Relative tolerance of shake");
2164 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2165 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2166 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2167 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2168 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2169 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2170 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2171 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2172 printStringNoNewline(&inp, "rotates over more degrees than");
2173 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2174 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2175 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
2177 /* Energy group exclusions */
2178 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2179 printStringNoNewline(
2180 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2181 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2184 printStringNewline(&inp, "WALLS");
2185 printStringNoNewline(
2186 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2187 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2188 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
2189 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2190 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2191 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2192 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2195 printStringNewline(&inp, "COM PULLING");
2196 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2200 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull, wi);
2204 for (int c = 0; c < ir->pull->ncoord; c++)
2206 if (ir->pull->coord[c].eType == epullCONSTRAINT)
2209 "Constraint COM pulling is not supported in combination with "
2210 "multiple time stepping");
2218 NOTE: needs COM pulling or free energy input */
2219 printStringNewline(&inp, "AWH biasing");
2220 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2223 ir->awhParams = gmx::readAwhParams(&inp, wi);
2226 /* Enforced rotation */
2227 printStringNewline(&inp, "ENFORCED ROTATION");
2228 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2229 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2233 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2236 /* Interactive MD */
2238 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2239 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2240 if (inputrecStrings->imd_grp[0] != '\0')
2247 printStringNewline(&inp, "NMR refinement stuff");
2248 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2249 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2250 printStringNoNewline(
2251 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2252 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2253 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2254 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2255 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2256 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2257 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2258 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2259 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2260 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2261 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2262 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2263 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2264 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2265 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2266 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2268 /* free energy variables */
2269 printStringNewline(&inp, "Free energy variables");
2270 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2271 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2272 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2273 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2274 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2276 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2278 it was not entered */
2279 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2280 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2281 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2282 setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
2283 setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
2284 setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
2285 setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
2286 setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
2287 setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
2288 setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
2289 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2290 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2291 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2292 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2293 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2294 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2295 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2296 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2297 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2298 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2299 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2300 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2301 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2302 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2304 /* Non-equilibrium MD stuff */
2305 printStringNewline(&inp, "Non-equilibrium MD stuff");
2306 setStringEntry(&inp, "acc-grps", inputrecStrings->accgrps, nullptr);
2307 setStringEntry(&inp, "accelerate", inputrecStrings->acc, nullptr);
2308 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2309 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2310 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2311 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2313 /* simulated tempering variables */
2314 printStringNewline(&inp, "simulated tempering variables");
2315 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2316 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2317 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2318 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2320 /* expanded ensemble variables */
2321 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2323 read_expandedparams(&inp, expand, wi);
2326 /* Electric fields */
2328 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2329 gmx::KeyValueTreeTransformer transform;
2330 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2331 mdModules->initMdpTransform(transform.rules());
2332 for (const auto& path : transform.mappedPaths())
2334 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2335 mark_einp_set(inp, path[0].c_str());
2337 MdpErrorHandler errorHandler(wi);
2338 auto result = transform.transform(convertedValues, &errorHandler);
2339 ir->params = new gmx::KeyValueTreeObject(result.object());
2340 mdModules->adjustInputrecBasedOnModules(ir);
2341 errorHandler.setBackMapping(result.backMapping());
2342 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2345 /* Ion/water position swapping ("computational electrophysiology") */
2346 printStringNewline(&inp,
2347 "Ion/water position swapping for computational electrophysiology setups");
2348 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2349 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2350 if (ir->eSwapCoords != eswapNO)
2357 printStringNoNewline(&inp, "Swap attempt frequency");
2358 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2359 printStringNoNewline(&inp, "Number of ion types to be controlled");
2360 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2363 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2365 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2366 snew(ir->swap->grp, ir->swap->ngrp);
2367 for (i = 0; i < ir->swap->ngrp; i++)
2369 snew(ir->swap->grp[i].molname, STRLEN);
2371 printStringNoNewline(&inp,
2372 "Two index groups that contain the compartment-partitioning atoms");
2373 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2374 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2375 printStringNoNewline(&inp,
2376 "Use center of mass of split groups (yes/no), otherwise center of "
2377 "geometry is used");
2378 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2379 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2381 printStringNoNewline(&inp, "Name of solvent molecules");
2382 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2384 printStringNoNewline(&inp,
2385 "Split cylinder: radius, upper and lower extension (nm) (this will "
2386 "define the channels)");
2387 printStringNoNewline(&inp,
2388 "Note that the split cylinder settings do not have an influence on "
2389 "the swapping protocol,");
2390 printStringNoNewline(
2392 "however, if correctly defined, the permeation events are recorded per channel");
2393 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2394 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2395 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2396 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2397 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2398 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2400 printStringNoNewline(
2402 "Average the number of ions per compartment over these many swap attempt steps");
2403 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2405 printStringNoNewline(
2406 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2407 printStringNoNewline(
2408 &inp, "and the requested number of ions of this type in compartments A and B");
2409 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2410 for (i = 0; i < nIonTypes; i++)
2412 int ig = eSwapFixedGrpNR + i;
2414 sprintf(buf, "iontype%d-name", i);
2415 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2416 sprintf(buf, "iontype%d-in-A", i);
2417 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2418 sprintf(buf, "iontype%d-in-B", i);
2419 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2422 printStringNoNewline(
2424 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2425 printStringNoNewline(
2427 "at maximum distance (= bulk concentration) to the split group layers. However,");
2428 printStringNoNewline(&inp,
2429 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2430 "layer from the middle at 0.0");
2431 printStringNoNewline(&inp,
2432 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2433 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2434 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2435 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2436 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2438 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2441 printStringNoNewline(
2442 &inp, "Start to swap ions if threshold difference to requested count is reached");
2443 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2446 /* AdResS is no longer supported, but we need grompp to be able to
2447 refuse to process old .mdp files that used it. */
2448 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2450 /* User defined thingies */
2451 printStringNewline(&inp, "User defined thingies");
2452 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2453 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2454 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2455 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2456 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2457 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2458 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2459 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2460 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2461 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2465 gmx::TextOutputFile stream(mdparout);
2466 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2468 // Transform module data into a flat key-value tree for output.
2469 gmx::KeyValueTreeBuilder builder;
2470 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2471 mdModules->buildMdpOutput(&builderObject);
2473 gmx::TextWriter writer(&stream);
2474 writeKeyValueTreeAsMdp(&writer, builder.build());
2479 /* Process options if necessary */
2480 for (m = 0; m < 2; m++)
2482 for (i = 0; i < 2 * DIM; i++)
2491 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2495 "Pressure coupling incorrect number of values (I need exactly 1)");
2497 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2499 case epctSEMIISOTROPIC:
2500 case epctSURFACETENSION:
2501 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2505 "Pressure coupling incorrect number of values (I need exactly 2)");
2507 dumdub[m][YY] = dumdub[m][XX];
2509 case epctANISOTROPIC:
2510 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf", &(dumdub[m][XX]), &(dumdub[m][YY]),
2511 &(dumdub[m][ZZ]), &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5]))
2516 "Pressure coupling incorrect number of values (I need exactly 6)");
2520 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2521 epcoupltype_names[ir->epct]);
2525 clear_mat(ir->ref_p);
2526 clear_mat(ir->compress);
2527 for (i = 0; i < DIM; i++)
2529 ir->ref_p[i][i] = dumdub[1][i];
2530 ir->compress[i][i] = dumdub[0][i];
2532 if (ir->epct == epctANISOTROPIC)
2534 ir->ref_p[XX][YY] = dumdub[1][3];
2535 ir->ref_p[XX][ZZ] = dumdub[1][4];
2536 ir->ref_p[YY][ZZ] = dumdub[1][5];
2537 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2540 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2541 "apply a threefold shear stress?\n");
2543 ir->compress[XX][YY] = dumdub[0][3];
2544 ir->compress[XX][ZZ] = dumdub[0][4];
2545 ir->compress[YY][ZZ] = dumdub[0][5];
2546 for (i = 0; i < DIM; i++)
2548 for (m = 0; m < i; m++)
2550 ir->ref_p[i][m] = ir->ref_p[m][i];
2551 ir->compress[i][m] = ir->compress[m][i];
2556 if (ir->comm_mode == ecmNO)
2561 opts->couple_moltype = nullptr;
2562 if (strlen(inputrecStrings->couple_moltype) > 0)
2564 if (ir->efep != efepNO)
2566 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2567 if (opts->couple_lam0 == opts->couple_lam1)
2569 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2571 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2575 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2582 "Free energy is turned off, so we will not decouple the molecule listed "
2586 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2587 if (ir->efep != efepNO)
2589 if (fep->delta_lambda != 0)
2591 ir->efep = efepSLOWGROWTH;
2595 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2597 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2599 "Old option for dhdl-print-energy given: "
2600 "changing \"yes\" to \"total\"\n");
2603 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2605 /* always print out the energy to dhdl if we are doing
2606 expanded ensemble, since we need the total energy for
2607 analysis if the temperature is changing. In some
2608 conditions one may only want the potential energy, so
2609 we will allow that if the appropriate mdp setting has
2610 been enabled. Otherwise, total it is:
2612 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2615 if ((ir->efep != efepNO) || ir->bSimTemp)
2617 ir->bExpanded = FALSE;
2618 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2620 ir->bExpanded = TRUE;
2622 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2623 if (ir->bSimTemp) /* done after fep params */
2625 do_simtemp_params(ir);
2628 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2629 * setup and not on the old way of specifying the free-energy setup,
2630 * we should check for using soft-core when not needed, since that
2631 * can complicate the sampling significantly.
2632 * Note that we only check for the automated coupling setup.
2633 * If the (advanced) user does FEP through manual topology changes,
2634 * this check will not be triggered.
2636 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 && ir->fepvals->sc_alpha != 0
2637 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2640 "You are using soft-core interactions while the Van der Waals interactions are "
2641 "not decoupled (note that the sc-coul option is only active when using lambda "
2642 "states). Although this will not lead to errors, you will need much more "
2643 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2648 ir->fepvals->n_lambda = 0;
2651 /* WALL PARAMETERS */
2653 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2655 /* ORIENTATION RESTRAINT PARAMETERS */
2657 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2659 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2662 /* DEFORMATION PARAMETERS */
2664 clear_mat(ir->deform);
2665 for (i = 0; i < 6; i++)
2670 double gmx_unused canary;
2671 int ndeform = sscanf(inputrecStrings->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]),
2672 &(dumdub[0][1]), &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]),
2673 &(dumdub[0][5]), &canary);
2675 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2679 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2680 inputrecStrings->deform)
2683 for (i = 0; i < 3; i++)
2685 ir->deform[i][i] = dumdub[0][i];
2687 ir->deform[YY][XX] = dumdub[0][3];
2688 ir->deform[ZZ][XX] = dumdub[0][4];
2689 ir->deform[ZZ][YY] = dumdub[0][5];
2690 if (ir->epc != epcNO)
2692 for (i = 0; i < 3; i++)
2694 for (j = 0; j <= i; j++)
2696 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2698 warning_error(wi, "A box element has deform set and compressibility > 0");
2702 for (i = 0; i < 3; i++)
2704 for (j = 0; j < i; j++)
2706 if (ir->deform[i][j] != 0)
2708 for (m = j; m < DIM; m++)
2710 if (ir->compress[m][j] != 0)
2713 "An off-diagonal box element has deform set while "
2714 "compressibility > 0 for the same component of another box "
2715 "vector, this might lead to spurious periodicity effects.");
2716 warning(wi, warn_buf);
2724 /* Ion/water position swapping checks */
2725 if (ir->eSwapCoords != eswapNO)
2727 if (ir->swap->nstswap < 1)
2729 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2731 if (ir->swap->nAverage < 1)
2733 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2735 if (ir->swap->threshold < 1.0)
2737 warning_error(wi, "Ion count threshold must be at least 1.\n");
2741 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2744 setupMtsLevels(ir->mtsLevels, *ir, *opts, wi);
2749 gmx::checkAwhParams(ir->awhParams, ir, wi);
2756 /* We would like gn to be const as well, but C doesn't allow this */
2757 /* TODO this is utility functionality (search for the index of a
2758 string in a collection), so should be refactored and located more
2760 int search_string(const char* s, int ng, char* gn[])
2764 for (i = 0; (i < ng); i++)
2766 if (gmx_strcasecmp(s, gn[i]) == 0)
2773 "Group %s referenced in the .mdp file was not found in the index file.\n"
2774 "Group names must match either [moleculetype] names or custom index group\n"
2775 "names, in which case you must supply an index file to the '-n' option\n"
2780 static void do_numbering(int natoms,
2781 SimulationGroups* groups,
2782 gmx::ArrayRef<std::string> groupsFromMdpFile,
2785 SimulationAtomGroupType gtype,
2791 unsigned short* cbuf;
2792 AtomGroupIndices* grps = &(groups->groups[gtype]);
2793 int j, gid, aj, ognr, ntot = 0;
2795 char warn_buf[STRLEN];
2797 title = shortName(gtype);
2800 /* Mark all id's as not set */
2801 for (int i = 0; (i < natoms); i++)
2806 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2808 /* Lookup the group name in the block structure */
2809 gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2810 if ((grptp != egrptpONE) || (i == 0))
2812 grps->emplace_back(gid);
2815 /* Now go over the atoms in the group */
2816 for (j = block->index[gid]; (j < block->index[gid + 1]); j++)
2821 /* Range checking */
2822 if ((aj < 0) || (aj >= natoms))
2824 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2826 /* Lookup up the old group number */
2830 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
2835 /* Store the group number in buffer */
2836 if (grptp == egrptpONE)
2849 /* Now check whether we have done all atoms */
2852 if (grptp == egrptpALL)
2854 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2856 else if (grptp == egrptpPART)
2858 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2859 warning_note(wi, warn_buf);
2861 /* Assign all atoms currently unassigned to a rest group */
2862 for (j = 0; (j < natoms); j++)
2864 if (cbuf[j] == NOGID)
2866 cbuf[j] = grps->size();
2869 if (grptp != egrptpPART)
2873 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title,
2876 /* Add group name "rest" */
2877 grps->emplace_back(restnm);
2879 /* Assign the rest name to all atoms not currently assigned to a group */
2880 for (j = 0; (j < natoms); j++)
2882 if (cbuf[j] == NOGID)
2884 // group size was not updated before this here, so need to use -1.
2885 cbuf[j] = grps->size() - 1;
2891 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2893 /* All atoms are part of one (or no) group, no index required */
2894 groups->groupNumbers[gtype].clear();
2898 for (int j = 0; (j < natoms); j++)
2900 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2907 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2910 pull_params_t* pull;
2911 int natoms, imin, jmin;
2912 int * nrdf2, *na_vcm, na_tot;
2913 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2918 * First calc 3xnr-atoms for each group
2919 * then subtract half a degree of freedom for each constraint
2921 * Only atoms and nuclei contribute to the degrees of freedom...
2926 const SimulationGroups& groups = mtop->groups;
2927 natoms = mtop->natoms;
2929 /* Allocate one more for a possible rest group */
2930 /* We need to sum degrees of freedom into doubles,
2931 * since floats give too low nrdf's above 3 million atoms.
2933 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
2934 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2935 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2936 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2937 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2939 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2943 for (gmx::index i = 0;
2944 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; i++)
2947 clear_ivec(dof_vcm[i]);
2949 nrdf_vcm_sub[i] = 0;
2951 snew(nrdf2, natoms);
2952 for (const AtomProxy atomP : AtomRange(*mtop))
2954 const t_atom& local = atomP.atom();
2955 int i = atomP.globalAtomNumber();
2957 if (local.ptype == eptAtom || local.ptype == eptNucleus)
2959 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
2960 for (int d = 0; d < DIM; d++)
2962 if (opts->nFreeze[g][d] == 0)
2964 /* Add one DOF for particle i (counted as 2*1) */
2966 /* VCM group i has dim d as a DOF */
2967 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
2971 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
2973 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
2979 for (const gmx_molblock_t& molb : mtop->molblock)
2981 const gmx_moltype_t& molt = mtop->moltype[molb.type];
2982 const t_atom* atom = molt.atoms.atom;
2983 for (int mol = 0; mol < molb.nmol; mol++)
2985 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2987 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
2988 for (int i = 0; i < molt.ilist[ftype].size();)
2990 /* Subtract degrees of freedom for the constraints,
2991 * if the particles still have degrees of freedom left.
2992 * If one of the particles is a vsite or a shell, then all
2993 * constraint motion will go there, but since they do not
2994 * contribute to the constraints the degrees of freedom do not
2997 int ai = as + ia[i + 1];
2998 int aj = as + ia[i + 2];
2999 if (((atom[ia[i + 1]].ptype == eptNucleus) || (atom[ia[i + 1]].ptype == eptAtom))
3000 && ((atom[ia[i + 2]].ptype == eptNucleus) || (atom[ia[i + 2]].ptype == eptAtom)))
3018 imin = std::min(imin, nrdf2[ai]);
3019 jmin = std::min(jmin, nrdf2[aj]);
3022 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3024 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3026 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3028 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3031 i += interaction_function[ftype].nratoms + 1;
3034 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3035 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3037 /* Subtract 1 dof from every atom in the SETTLE */
3038 for (int j = 0; j < 3; j++)
3040 int ai = as + ia[i + 1 + j];
3041 imin = std::min(2, nrdf2[ai]);
3043 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3045 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3050 as += molt.atoms.nr;
3056 /* Correct nrdf for the COM constraints.
3057 * We correct using the TC and VCM group of the first atom
3058 * in the reference and pull group. If atoms in one pull group
3059 * belong to different TC or VCM groups it is anyhow difficult
3060 * to determine the optimal nrdf assignment.
3064 for (int i = 0; i < pull->ncoord; i++)
3066 if (pull->coord[i].eType != epullCONSTRAINT)
3073 for (int j = 0; j < 2; j++)
3075 const t_pull_group* pgrp;
3077 pgrp = &pull->group[pull->coord[i].group[j]];
3081 /* Subtract 1/2 dof from each group */
3082 int ai = pgrp->ind[0];
3083 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3085 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3087 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3090 "Center of mass pulling constraints caused the number of degrees "
3091 "of freedom for temperature coupling group %s to be negative",
3092 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3093 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3098 /* We need to subtract the whole DOF from group j=1 */
3105 if (ir->nstcomm != 0)
3107 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3108 "Expect at least one group when removing COM motion");
3110 /* We remove COM motion up to dim ndof_com() */
3111 const int ndim_rm_vcm = ndof_com(ir);
3113 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3114 * the number of degrees of freedom in each vcm group when COM
3115 * translation is removed and 6 when rotation is removed as well.
3116 * Note that we do not and should not include the rest group here.
3118 for (gmx::index j = 0;
3119 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]); j++)
3121 switch (ir->comm_mode)
3124 case ecmLINEAR_ACCELERATION_CORRECTION:
3125 nrdf_vcm_sub[j] = 0;
3126 for (int d = 0; d < ndim_rm_vcm; d++)
3134 case ecmANGULAR: nrdf_vcm_sub[j] = 6; break;
3135 default: gmx_incons("Checking comm_mode");
3139 for (gmx::index i = 0;
3140 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
3142 /* Count the number of atoms of TC group i for every VCM group */
3143 for (gmx::index j = 0;
3144 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3149 for (int ai = 0; ai < natoms; ai++)
3151 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3153 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3157 /* Correct for VCM removal according to the fraction of each VCM
3158 * group present in this TC group.
3160 nrdf_uc = nrdf_tc[i];
3162 for (gmx::index j = 0;
3163 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3165 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3167 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3168 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3173 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3175 opts->nrdf[i] = nrdf_tc[i];
3176 if (opts->nrdf[i] < 0)
3180 fprintf(stderr, "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3181 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
3189 sfree(nrdf_vcm_sub);
3192 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3194 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3195 * But since this is much larger than STRLEN, such a line can not be parsed.
3196 * The real maximum is the number of names that fit in a string: STRLEN/2.
3198 #define EGP_MAX (STRLEN / 2)
3202 auto names = gmx::splitString(val);
3203 if (names.size() % 2 != 0)
3205 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3207 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3209 for (size_t i = 0; i < names.size() / 2; i++)
3211 // TODO this needs to be replaced by a solution using std::find_if
3215 names[2 * i].c_str(),
3216 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3222 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3227 names[2 * i + 1].c_str(),
3228 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3234 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3236 if ((j < nr) && (k < nr))
3238 ir->opts.egp_flags[nr * j + k] |= flag;
3239 ir->opts.egp_flags[nr * k + j] |= flag;
3248 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3250 int ig = -1, i = 0, gind;
3254 /* Just a quick check here, more thorough checks are in mdrun */
3255 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3257 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3260 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3261 for (ig = 0; ig < swap->ngrp; ig++)
3263 swapg = &swap->grp[ig];
3264 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3265 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3269 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3270 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap", swap->grp[ig].molname, swapg->nat);
3271 snew(swapg->ind, swapg->nat);
3272 for (i = 0; i < swapg->nat; i++)
3274 swapg->ind[i] = grps->a[grps->index[gind] + i];
3279 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3285 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3290 ig = search_string(IMDgname, grps->nr, gnames);
3291 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3293 if (IMDgroup->nat > 0)
3296 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3298 IMDgname, IMDgroup->nat);
3299 snew(IMDgroup->ind, IMDgroup->nat);
3300 for (i = 0; i < IMDgroup->nat; i++)
3302 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3307 /* Checks whether atoms are both part of a COM removal group and frozen.
3308 * If a fully frozen atom is part of a COM removal group, it is removed
3309 * from the COM removal group. A note is issued if such atoms are present.
3310 * A warning is issued for atom with one or two dimensions frozen that
3311 * are part of a COM removal group (mdrun would need to compute COM mass
3312 * per dimension to handle this correctly).
3313 * Also issues a warning when non-frozen atoms are not part of a COM
3314 * removal group while COM removal is active.
3316 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3318 const t_grpopts& opts,
3321 const int vcmRestGroup =
3322 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3324 int numFullyFrozenVcmAtoms = 0;
3325 int numPartiallyFrozenVcmAtoms = 0;
3326 int numNonVcmAtoms = 0;
3327 for (int a = 0; a < numAtoms; a++)
3329 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3330 int numFrozenDims = 0;
3331 for (int d = 0; d < DIM; d++)
3333 numFrozenDims += opts.nFreeze[freezeGroup][d];
3336 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3337 if (vcmGroup < vcmRestGroup)
3339 if (numFrozenDims == DIM)
3341 /* Do not remove COM motion for this fully frozen atom */
3342 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3344 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3347 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3348 numFullyFrozenVcmAtoms++;
3350 else if (numFrozenDims > 0)
3352 numPartiallyFrozenVcmAtoms++;
3355 else if (numFrozenDims < DIM)
3361 if (numFullyFrozenVcmAtoms > 0)
3363 std::string warningText = gmx::formatString(
3364 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3365 "removing these atoms from the COMM removal group(s)",
3366 numFullyFrozenVcmAtoms);
3367 warning_note(wi, warningText.c_str());
3369 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3371 std::string warningText = gmx::formatString(
3372 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3373 "removal group(s), due to limitations in the code these still contribute to the "
3374 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3376 numPartiallyFrozenVcmAtoms, DIM);
3377 warning(wi, warningText.c_str());
3379 if (numNonVcmAtoms > 0)
3381 std::string warningText = gmx::formatString(
3382 "%d atoms are not part of any center of mass motion removal group.\n"
3383 "This may lead to artifacts.\n"
3384 "In most cases one should use one group for the whole system.",
3386 warning(wi, warningText.c_str());
3390 void do_index(const char* mdparin,
3394 const gmx::MdModulesNotifier& notifier,
3398 t_blocka* defaultIndexGroups;
3406 int i, j, k, restnm;
3407 bool bExcl, bTable, bAnneal;
3408 char warn_buf[STRLEN];
3412 fprintf(stderr, "processing index file...\n");
3416 snew(defaultIndexGroups, 1);
3417 snew(defaultIndexGroups->index, 1);
3419 atoms_all = gmx_mtop_global_atoms(mtop);
3420 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3421 done_atom(&atoms_all);
3425 defaultIndexGroups = init_index(ndx, &gnames);
3428 SimulationGroups* groups = &mtop->groups;
3429 natoms = mtop->natoms;
3430 symtab = &mtop->symtab;
3432 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3434 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3436 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3437 restnm = groups->groupNames.size() - 1;
3438 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3439 srenew(gnames, defaultIndexGroups->nr + 1);
3440 gnames[restnm] = *(groups->groupNames.back());
3442 set_warning_line(wi, mdparin, -1);
3444 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3445 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3446 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3447 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3448 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3451 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3453 temperatureCouplingGroupNames.size(), temperatureCouplingReferenceValues.size(),
3454 temperatureCouplingTauValues.size());
3457 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3458 do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
3459 SimulationAtomGroupType::TemperatureCoupling, restnm,
3460 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3461 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3463 snew(ir->opts.nrdf, nr);
3464 snew(ir->opts.tau_t, nr);
3465 snew(ir->opts.ref_t, nr);
3466 if (ir->eI == eiBD && ir->bd_fric == 0)
3468 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3471 if (useReferenceTemperature)
3473 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3475 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3479 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3480 for (i = 0; (i < nr); i++)
3482 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3484 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3485 warning_error(wi, warn_buf);
3488 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3492 "tau-t = -1 is the value to signal that a group should not have "
3493 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3496 if (ir->opts.tau_t[i] >= 0)
3498 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3501 if (ir->etc != etcNO && ir->nsttcouple == -1)
3503 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3508 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3511 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3512 "md-vv; use either vrescale temperature with berendsen pressure or "
3513 "Nose-Hoover temperature with MTTK pressure");
3515 if (ir->epc == epcMTTK)
3517 if (ir->etc != etcNOSEHOOVER)
3520 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3525 if (ir->nstpcouple != ir->nsttcouple)
3527 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3528 ir->nstpcouple = ir->nsttcouple = mincouple;
3530 "for current Trotter decomposition methods with vv, nsttcouple and "
3531 "nstpcouple must be equal. Both have been reset to "
3532 "min(nsttcouple,nstpcouple) = %d",
3534 warning_note(wi, warn_buf);
3539 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3540 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3542 if (ETC_ANDERSEN(ir->etc))
3544 if (ir->nsttcouple != 1)
3548 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3549 "need for larger nsttcouple > 1, since no global parameters are computed. "
3550 "nsttcouple has been reset to 1");
3551 warning_note(wi, warn_buf);
3554 nstcmin = tcouple_min_integration_steps(ir->etc);
3557 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3560 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3561 "least %d times larger than nsttcouple*dt (%g)",
3562 ETCOUPLTYPE(ir->etc), tau_min, nstcmin, ir->nsttcouple * ir->delta_t);
3563 warning(wi, warn_buf);
3566 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3567 for (i = 0; (i < nr); i++)
3569 if (ir->opts.ref_t[i] < 0)
3571 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3574 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3575 if we are in this conditional) if mc_temp is negative */
3576 if (ir->expandedvals->mc_temp < 0)
3578 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3582 /* Simulated annealing for each group. There are nr groups */
3583 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3584 if (simulatedAnnealingGroupNames.size() == 1
3585 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3587 simulatedAnnealingGroupNames.resize(0);
3589 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3591 gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
3592 simulatedAnnealingGroupNames.size(), nr);
3596 snew(ir->opts.annealing, nr);
3597 snew(ir->opts.anneal_npoints, nr);
3598 snew(ir->opts.anneal_time, nr);
3599 snew(ir->opts.anneal_temp, nr);
3600 for (i = 0; i < nr; i++)
3602 ir->opts.annealing[i] = eannNO;
3603 ir->opts.anneal_npoints[i] = 0;
3604 ir->opts.anneal_time[i] = nullptr;
3605 ir->opts.anneal_temp[i] = nullptr;
3607 if (!simulatedAnnealingGroupNames.empty())
3610 for (i = 0; i < nr; i++)
3612 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3614 ir->opts.annealing[i] = eannNO;
3616 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3618 ir->opts.annealing[i] = eannSINGLE;
3621 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3623 ir->opts.annealing[i] = eannPERIODIC;
3629 /* Read the other fields too */
3630 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3631 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3633 gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
3634 simulatedAnnealingPoints.size(), simulatedAnnealingGroupNames.size());
3636 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3637 size_t numSimulatedAnnealingFields = 0;
3638 for (i = 0; i < nr; i++)
3640 if (ir->opts.anneal_npoints[i] == 1)
3644 "Please specify at least a start and an end point for annealing\n");
3646 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3647 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3648 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3651 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3653 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3655 gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
3656 simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
3658 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3659 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3661 gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
3662 simulatedAnnealingTemperatures.size(), numSimulatedAnnealingFields);
3665 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3666 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3667 convertReals(wi, simulatedAnnealingTimes, "anneal-time",
3668 allSimulatedAnnealingTimes.data());
3669 convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp",
3670 allSimulatedAnnealingTemperatures.data());
3671 for (i = 0, k = 0; i < nr; i++)
3673 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3675 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3676 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3679 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3681 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3687 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3690 "Annealing timepoints out of order: t=%f comes after "
3692 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j - 1]);
3695 if (ir->opts.anneal_temp[i][j] < 0)
3697 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n",
3698 ir->opts.anneal_temp[i][j]);
3703 /* Print out some summary information, to make sure we got it right */
3704 for (i = 0; i < nr; i++)
3706 if (ir->opts.annealing[i] != eannNO)
3708 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3709 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3710 *(groups->groupNames[j]), eann_names[ir->opts.annealing[i]],
3711 ir->opts.anneal_npoints[i]);
3712 fprintf(stderr, "Time (ps) Temperature (K)\n");
3713 /* All terms except the last one */
3714 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3716 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3717 ir->opts.anneal_temp[i][j]);
3720 /* Finally the last one */
3721 j = ir->opts.anneal_npoints[i] - 1;
3722 if (ir->opts.annealing[i] == eannSINGLE)
3724 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j],
3725 ir->opts.anneal_temp[i][j]);
3729 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3730 ir->opts.anneal_temp[i][j]);
3731 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3734 "There is a temperature jump when your annealing "
3746 make_pull_groups(ir->pull, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3748 make_pull_coords(ir->pull);
3753 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3756 if (ir->eSwapCoords != eswapNO)
3758 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3761 /* Make indices for IMD session */
3764 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3767 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3768 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3769 notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
3771 auto accelerations = gmx::splitString(inputrecStrings->acc);
3772 auto accelerationGroupNames = gmx::splitString(inputrecStrings->accgrps);
3773 if (accelerationGroupNames.size() * DIM != accelerations.size())
3775 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3776 accelerationGroupNames.size(), accelerations.size());
3778 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3779 SimulationAtomGroupType::Acceleration, restnm, egrptpALL_GENREST, bVerbose, wi);
3780 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3781 snew(ir->opts.acc, nr);
3782 ir->opts.ngacc = nr;
3784 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3786 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3787 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3788 if (freezeDims.size() != DIM * freezeGroupNames.size())
3790 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3791 freezeGroupNames.size(), freezeDims.size());
3793 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3794 SimulationAtomGroupType::Freeze, restnm, egrptpALL_GENREST, bVerbose, wi);
3795 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3796 ir->opts.ngfrz = nr;
3797 snew(ir->opts.nFreeze, nr);
3798 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3800 for (j = 0; (j < DIM); j++, k++)
3802 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3803 if (!ir->opts.nFreeze[i][j])
3805 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3808 "Please use Y(ES) or N(O) for freezedim only "
3810 freezeDims[k].c_str());
3811 warning(wi, warn_buf);
3816 for (; (i < nr); i++)
3818 for (j = 0; (j < DIM); j++)
3820 ir->opts.nFreeze[i][j] = 0;
3824 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3825 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3826 SimulationAtomGroupType::EnergyOutput, restnm, egrptpALL_GENREST, bVerbose, wi);
3827 add_wall_energrps(groups, ir->nwall, symtab);
3828 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3829 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3830 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3831 SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
3832 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3834 if (ir->comm_mode != ecmNO)
3836 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3839 /* Now we have filled the freeze struct, so we can calculate NRDF */
3840 calc_nrdf(mtop, ir, gnames);
3842 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3843 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3844 SimulationAtomGroupType::User1, restnm, egrptpALL_GENREST, bVerbose, wi);
3845 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3846 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3847 SimulationAtomGroupType::User2, restnm, egrptpALL_GENREST, bVerbose, wi);
3848 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3849 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3850 SimulationAtomGroupType::CompressedPositionOutput, restnm, egrptpONE, bVerbose, wi);
3851 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3852 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3853 SimulationAtomGroupType::OrientationRestraintsFit, restnm, egrptpALL_GENREST,
3856 /* MiMiC QMMM input processing */
3857 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
3858 if (qmGroupNames.size() > 1)
3860 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3862 /* group rest, if any, is always MM! */
3863 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3864 SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
3865 ir->opts.ngQM = qmGroupNames.size();
3867 /* end of MiMiC QMMM input */
3871 for (auto group : gmx::keysOf(groups->groups))
3873 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3874 for (const auto& entry : groups->groups[group])
3876 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3878 fprintf(stderr, "\n");
3882 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3883 snew(ir->opts.egp_flags, nr * nr);
3885 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
3886 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3888 warning_error(wi, "Energy group exclusions are currently not supported");
3890 if (bExcl && EEL_FULL(ir->coulombtype))
3892 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3895 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
3896 if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
3897 && !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
3900 "Can only have energy group pair tables in combination with user tables for VdW "
3904 /* final check before going out of scope if simulated tempering variables
3905 * need to be set to default values.
3907 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3909 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
3910 warning(wi, gmx::formatString(
3911 "the value for nstexpanded was not specified for "
3912 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3913 "by default, but it is recommended to set it to an explicit value!",
3914 ir->expandedvals->nstexpanded));
3916 for (i = 0; (i < defaultIndexGroups->nr); i++)
3921 done_blocka(defaultIndexGroups);
3922 sfree(defaultIndexGroups);
3926 static void check_disre(const gmx_mtop_t* mtop)
3928 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3930 const gmx_ffparams_t& ffparams = mtop->ffparams;
3933 for (int i = 0; i < ffparams.numTypes(); i++)
3935 int ftype = ffparams.functype[i];
3936 if (ftype == F_DISRES)
3938 int label = ffparams.iparams[i].disres.label;
3939 if (label == old_label)
3941 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3950 "Found %d double distance restraint indices,\n"
3951 "probably the parameters for multiple pairs in one restraint "
3952 "are not identical\n",
3958 static bool absolute_reference(const t_inputrec* ir, const gmx_mtop_t* sys, const bool posres_only, ivec AbsRef)
3961 gmx_mtop_ilistloop_t iloop;
3963 const t_iparams* pr;
3970 for (d = 0; d < DIM; d++)
3972 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3974 /* Check for freeze groups */
3975 for (g = 0; g < ir->opts.ngfrz; g++)
3977 for (d = 0; d < DIM; d++)
3979 if (ir->opts.nFreeze[g][d] != 0)
3987 /* Check for position restraints */
3988 iloop = gmx_mtop_ilistloop_init(sys);
3989 while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
3991 if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3993 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
3995 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
3996 for (d = 0; d < DIM; d++)
3998 if (pr->posres.fcA[d] != 0)
4004 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
4006 /* Check for flat-bottom posres */
4007 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
4008 if (pr->fbposres.k != 0)
4010 switch (pr->fbposres.geom)
4012 case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
4013 case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
4014 case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
4015 case efbposresCYLINDER:
4016 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4017 case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
4018 case efbposresX: /* d=XX */
4019 case efbposresY: /* d=YY */
4020 case efbposresZ: /* d=ZZ */
4021 d = pr->fbposres.geom - efbposresX;
4026 " Invalid geometry for flat-bottom position restraint.\n"
4027 "Expected nr between 1 and %d. Found %d\n",
4028 efbposresNR - 1, pr->fbposres.geom);
4035 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
4038 static void check_combination_rule_differences(const gmx_mtop_t* mtop,
4040 bool* bC6ParametersWorkWithGeometricRules,
4041 bool* bC6ParametersWorkWithLBRules,
4042 bool* bLBRulesPossible)
4044 int ntypes, tpi, tpj;
4047 double c6i, c6j, c12i, c12j;
4048 double c6, c6_geometric, c6_LB;
4049 double sigmai, sigmaj, epsi, epsj;
4050 bool bCanDoLBRules, bCanDoGeometricRules;
4053 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4054 * force-field floating point parameters.
4057 ptr = getenv("GMX_LJCOMB_TOL");
4061 double gmx_unused canary;
4063 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4066 "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4071 *bC6ParametersWorkWithLBRules = TRUE;
4072 *bC6ParametersWorkWithGeometricRules = TRUE;
4073 bCanDoLBRules = TRUE;
4074 ntypes = mtop->ffparams.atnr;
4075 snew(typecount, ntypes);
4076 gmx_mtop_count_atomtypes(mtop, state, typecount);
4077 *bLBRulesPossible = TRUE;
4078 for (tpi = 0; tpi < ntypes; ++tpi)
4080 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4081 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4082 for (tpj = tpi; tpj < ntypes; ++tpj)
4084 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4085 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4086 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4087 c6_geometric = std::sqrt(c6i * c6j);
4088 if (!gmx_numzero(c6_geometric))
4090 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4092 sigmai = gmx::sixthroot(c12i / c6i);
4093 sigmaj = gmx::sixthroot(c12j / c6j);
4094 epsi = c6i * c6i / (4.0 * c12i);
4095 epsj = c6j * c6j / (4.0 * c12j);
4096 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4100 *bLBRulesPossible = FALSE;
4101 c6_LB = c6_geometric;
4103 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4108 *bC6ParametersWorkWithLBRules = FALSE;
4111 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4113 if (!bCanDoGeometricRules)
4115 *bC6ParametersWorkWithGeometricRules = FALSE;
4122 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop, warninp_t wi)
4124 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4126 check_combination_rule_differences(mtop, 0, &bC6ParametersWorkWithGeometricRules,
4127 &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4128 if (ir->ljpme_combination_rule == eljpmeLB)
4130 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4133 "You are using arithmetic-geometric combination rules "
4134 "in LJ-PME, but your non-bonded C6 parameters do not "
4135 "follow these rules.");
4140 if (!bC6ParametersWorkWithGeometricRules)
4142 if (ir->eDispCorr != edispcNO)
4145 "You are using geometric combination rules in "
4146 "LJ-PME, but your non-bonded C6 parameters do "
4147 "not follow these rules. "
4148 "This will introduce very small errors in the forces and energies in "
4149 "your simulations. Dispersion correction will correct total energy "
4150 "and/or pressure for isotropic systems, but not forces or surface "
4156 "You are using geometric combination rules in "
4157 "LJ-PME, but your non-bonded C6 parameters do "
4158 "not follow these rules. "
4159 "This will introduce very small errors in the forces and energies in "
4160 "your simulations. If your system is homogeneous, consider using "
4161 "dispersion correction "
4162 "for the total energy and pressure.");
4168 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4170 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4171 gmx::assertMtsRequirements(*ir);
4173 char err_buf[STRLEN];
4178 gmx_mtop_atomloop_block_t aloopb;
4180 char warn_buf[STRLEN];
4182 set_warning_line(wi, mdparin, -1);
4184 if (absolute_reference(ir, sys, false, AbsRef))
4187 "Removing center of mass motion in the presence of position restraints might "
4188 "cause artifacts. When you are using position restraints to equilibrate a "
4189 "macro-molecule, the artifacts are usually negligible.");
4192 if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 && ir->nstlist > 1
4193 && ((EI_MD(ir->eI) || EI_SD(ir->eI)) && (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4195 /* Check if a too small Verlet buffer might potentially
4196 * cause more drift than the thermostat can couple off.
4198 /* Temperature error fraction for warning and suggestion */
4199 const real T_error_warn = 0.002;
4200 const real T_error_suggest = 0.001;
4201 /* For safety: 2 DOF per atom (typical with constraints) */
4202 const real nrdf_at = 2;
4203 real T, tau, max_T_error;
4208 for (i = 0; i < ir->opts.ngtc; i++)
4210 T = std::max(T, ir->opts.ref_t[i]);
4211 tau = std::max(tau, ir->opts.tau_t[i]);
4215 /* This is a worst case estimate of the temperature error,
4216 * assuming perfect buffer estimation and no cancelation
4217 * of errors. The factor 0.5 is because energy distributes
4218 * equally over Ekin and Epot.
4220 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * BOLTZ * T);
4221 if (max_T_error > T_error_warn)
4224 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4225 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4226 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4227 "%.0e or decrease tau_t.",
4228 ir->verletbuf_tol, T, tau, 100 * max_T_error, 100 * T_error_suggest,
4229 ir->verletbuf_tol * T_error_suggest / max_T_error);
4230 warning(wi, warn_buf);
4235 if (ETC_ANDERSEN(ir->etc))
4239 for (i = 0; i < ir->opts.ngtc; i++)
4242 "all tau_t must currently be equal using Andersen temperature control, "
4243 "violated for group %d",
4245 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4247 "all tau_t must be positive using Andersen temperature control, "
4249 i, ir->opts.tau_t[i]);
4250 CHECK(ir->opts.tau_t[i] < 0);
4253 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
4255 for (i = 0; i < ir->opts.ngtc; i++)
4257 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4259 "tau_t/delta_t for group %d for temperature control method %s must be a "
4260 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4261 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4263 i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
4264 CHECK(nsteps % ir->nstcomm != 0);
4269 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD && ir->comm_mode == ecmNO
4270 && !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
4273 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4274 "rounding errors can lead to build up of kinetic energy of the center of mass");
4277 if (ir->epc == epcPARRINELLORAHMAN && ir->etc == etcNOSEHOOVER)
4280 for (int g = 0; g < ir->opts.ngtc; g++)
4282 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4284 if (ir->tau_p < 1.9 * tau_t_max)
4286 std::string message = gmx::formatString(
4287 "With %s T-coupling and %s p-coupling, "
4288 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4289 etcoupl_names[ir->etc], epcoupl_names[ir->epc], "tau-p", ir->tau_p, "tau-t",
4291 warning(wi, message.c_str());
4295 /* Check for pressure coupling with absolute position restraints */
4296 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4298 absolute_reference(ir, sys, TRUE, AbsRef);
4300 for (m = 0; m < DIM; m++)
4302 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4305 "You are using pressure coupling with absolute position restraints, "
4306 "this will give artifacts. Use the refcoord_scaling option.");
4314 aloopb = gmx_mtop_atomloop_block_init(sys);
4316 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4318 if (atom->q != 0 || atom->qB != 0)
4326 if (EEL_FULL(ir->coulombtype))
4329 "You are using full electrostatics treatment %s for a system without charges.\n"
4330 "This costs a lot of performance for just processing zeros, consider using %s "
4332 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4333 warning(wi, err_buf);
4338 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4341 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4342 "You might want to consider using %s electrostatics.\n",
4344 warning_note(wi, err_buf);
4348 /* Check if combination rules used in LJ-PME are the same as in the force field */
4349 if (EVDW_PME(ir->vdwtype))
4351 check_combination_rules(ir, sys, wi);
4354 /* Generalized reaction field */
4355 if (ir->coulombtype == eelGRF_NOTUSED)
4358 "Generalized reaction-field electrostatics is no longer supported. "
4359 "You can use normal reaction-field instead and compute the reaction-field "
4360 "constant by hand.");
4364 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4366 for (m = 0; (m < DIM); m++)
4368 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4377 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4378 for (const AtomProxy atomP : AtomRange(*sys))
4380 const t_atom& local = atomP.atom();
4381 int i = atomP.globalAtomNumber();
4382 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4385 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4387 for (m = 0; (m < DIM); m++)
4389 acc[m] += ir->opts.acc[i][m] * mgrp[i];
4393 for (m = 0; (m < DIM); m++)
4395 if (fabs(acc[m]) > 1e-6)
4397 const char* dim[DIM] = { "X", "Y", "Z" };
4398 fprintf(stderr, "Net Acceleration in %s direction, will %s be corrected\n", dim[m],
4399 ir->nstcomm != 0 ? "" : "not");
4400 if (ir->nstcomm != 0 && m < ndof_com(ir))
4404 (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4406 ir->opts.acc[i][m] -= acc[m];
4414 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0
4415 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4417 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4425 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4427 if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
4429 absolute_reference(ir, sys, FALSE, AbsRef);
4430 for (m = 0; m < DIM; m++)
4432 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4435 "You are using an absolute reference for pulling, but the rest of "
4436 "the system does not have an absolute reference. This will lead to "
4445 for (i = 0; i < 3; i++)
4447 for (m = 0; m <= i; m++)
4449 if ((ir->epc != epcNO && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4451 for (c = 0; c < ir->pull->ncoord; c++)
4453 if (ir->pull->coord[c].eGeom == epullgDIRPBC && ir->pull->coord[c].vec[m] != 0)
4456 "Can not have dynamic box while using pull geometry '%s' "
4458 EPULLGEOM(ir->pull->coord[c].eGeom), 'x' + m);
4469 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4471 char warn_buf[STRLEN];
4474 ptr = check_box(ir->pbcType, box);
4477 warning_error(wi, ptr);
4480 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4482 if (ir->shake_tol <= 0.0)
4484 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4485 warning_error(wi, warn_buf);
4489 if ((ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4491 /* If we have Lincs constraints: */
4492 if (ir->eI == eiMD && ir->etc == etcNO && ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4495 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4496 warning_note(wi, warn_buf);
4499 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4502 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4504 warning_note(wi, warn_buf);
4506 if (ir->epc == epcMTTK)
4508 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4512 if (bHasAnyConstraints && ir->epc == epcMTTK)
4514 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4517 if (ir->LincsWarnAngle > 90.0)
4519 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4520 warning(wi, warn_buf);
4521 ir->LincsWarnAngle = 90.0;
4524 if (ir->pbcType != PbcType::No)
4526 if (ir->nstlist == 0)
4529 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4530 "atoms might cause the simulation to crash.");
4532 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4535 "ERROR: The cut-off length is longer than half the shortest box vector or "
4536 "longer than the smallest box diagonal element. Increase the box size or "
4537 "decrease rlist.\n");
4538 warning_error(wi, warn_buf);