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.
51 #include "gromacs/applied_forces/awh/read_params.h"
52 #include "gromacs/fileio/readinp.h"
53 #include "gromacs/fileio/warninp.h"
54 #include "gromacs/gmxlib/network.h"
55 #include "gromacs/gmxpreprocess/toputil.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/calc_verletbuf.h"
60 #include "gromacs/mdrun/mdmodules.h"
61 #include "gromacs/mdtypes/inputrec.h"
62 #include "gromacs/mdtypes/md_enums.h"
63 #include "gromacs/mdtypes/multipletimestepping.h"
64 #include "gromacs/mdtypes/pull_params.h"
65 #include "gromacs/options/options.h"
66 #include "gromacs/options/treesupport.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/selection/indexutil.h"
69 #include "gromacs/topology/block.h"
70 #include "gromacs/topology/ifunc.h"
71 #include "gromacs/topology/index.h"
72 #include "gromacs/topology/mtop_util.h"
73 #include "gromacs/topology/symtab.h"
74 #include "gromacs/topology/topology.h"
75 #include "gromacs/utility/cstringutil.h"
76 #include "gromacs/utility/exceptions.h"
77 #include "gromacs/utility/fatalerror.h"
78 #include "gromacs/utility/filestream.h"
79 #include "gromacs/utility/gmxassert.h"
80 #include "gromacs/utility/ikeyvaluetreeerror.h"
81 #include "gromacs/utility/keyvaluetree.h"
82 #include "gromacs/utility/keyvaluetreebuilder.h"
83 #include "gromacs/utility/keyvaluetreemdpwriter.h"
84 #include "gromacs/utility/keyvaluetreetransform.h"
85 #include "gromacs/utility/mdmodulenotification.h"
86 #include "gromacs/utility/smalloc.h"
87 #include "gromacs/utility/strconvert.h"
88 #include "gromacs/utility/stringcompare.h"
89 #include "gromacs/utility/stringutil.h"
90 #include "gromacs/utility/textwriter.h"
95 /* Resource parameters
96 * Do not change any of these until you read the instruction
97 * in readinp.h. Some cpp's do not take spaces after the backslash
98 * (like the c-shell), which will give you a very weird compiler
102 struct gmx_inputrec_strings
104 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], acc[STRLEN], accgrps[STRLEN], freeze[STRLEN],
105 frdim[STRLEN], energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN],
106 x_compressed_groups[STRLEN], couple_moltype[STRLEN], orirefitgrp[STRLEN],
107 egptable[STRLEN], egpexcl[STRLEN], wall_atomtype[STRLEN], wall_density[STRLEN],
108 deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
109 char fep_lambda[efptNR][STRLEN];
110 char lambda_weights[STRLEN];
111 std::vector<std::string> pullGroupNames;
112 std::vector<std::string> rotateGroupNames;
113 char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
116 static gmx_inputrec_strings* inputrecStrings = nullptr;
118 void init_inputrec_strings()
123 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
124 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
126 inputrecStrings = new gmx_inputrec_strings();
129 void done_inputrec_strings()
131 delete inputrecStrings;
132 inputrecStrings = nullptr;
138 egrptpALL, /* All particles have to be a member of a group. */
139 egrptpALL_GENREST, /* A rest group with name is generated for particles *
140 * that are not part of any group. */
141 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
142 * for the rest group. */
143 egrptpONE /* Merge all selected groups into one group, *
144 * make a rest group for the remaining particles. */
147 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
148 "h-angles", "all-angles", nullptr };
150 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
152 static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lambdas)
157 for (i = 0; i < ntemps; i++)
159 /* simple linear scaling -- allows more control */
160 if (simtemp->eSimTempScale == esimtempLINEAR)
162 simtemp->temperatures[i] =
164 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
166 else if (simtemp->eSimTempScale
167 == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
169 simtemp->temperatures[i] = simtemp->simtemp_low
170 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
171 static_cast<real>((1.0 * i) / (ntemps - 1)));
173 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
175 simtemp->temperatures[i] = simtemp->simtemp_low
176 + (simtemp->simtemp_high - simtemp->simtemp_low)
177 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
182 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
183 gmx_fatal(FARGS, "%s", errorstr);
189 static void _low_check(bool b, const char* s, warninp_t wi)
193 warning_error(wi, s);
197 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
201 if (*p > 0 && *p % nst != 0)
203 /* Round up to the next multiple of nst */
204 *p = ((*p) / nst + 1) * nst;
205 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
210 static int lcd(int n1, int n2)
215 for (i = 2; (i <= n1 && i <= n2); i++)
217 if (n1 % i == 0 && n2 % i == 0)
226 //! Convert legacy mdp entries to modern ones.
227 static void process_interaction_modifier(int* eintmod)
229 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
231 *eintmod = eintmodPOTSHIFT;
235 static void checkMtsRequirement(const t_inputrec& ir, const char* param, const int nstValue, warninp_t wi)
237 GMX_RELEASE_ASSERT(ir.mtsLevels.size() >= 2, "Need at least two levels for MTS");
238 const int mtsFactor = ir.mtsLevels.back().stepFactor;
239 if (nstValue % mtsFactor != 0)
241 auto message = gmx::formatString(
242 "With MTS, %s = %d should be a multiple of mts-factor = %d", param, nstValue, mtsFactor);
243 warning_error(wi, message.c_str());
247 static void setupMtsLevels(gmx::ArrayRef<gmx::MtsLevel> mtsLevels,
248 const t_inputrec& ir,
249 const t_gromppopts& opts,
252 if (!(ir.eI == eiMD || ir.eI == eiSD1))
254 auto message = gmx::formatString(
255 "Multiple time stepping is only supported with integrators %s and %s",
256 ei_names[eiMD], ei_names[eiSD1]);
257 warning_error(wi, message.c_str());
259 if (opts.numMtsLevels != 2)
261 warning_error(wi, "Only mts-levels = 2 is supported");
265 const std::vector<std::string> inputForceGroups = gmx::splitString(opts.mtsLevel2Forces);
266 auto& forceGroups = mtsLevels[1].forceGroups;
267 for (const auto& inputForceGroup : inputForceGroups)
271 for (const auto& forceGroupName : gmx::mtsForceGroupNames)
273 if (gmx::equalCaseInsensitive(inputForceGroup, forceGroupName))
275 forceGroups.set(nameIndex);
283 gmx::formatString("Unknown MTS force group '%s'", inputForceGroup.c_str());
284 warning_error(wi, message.c_str());
288 if (mtsLevels[1].stepFactor <= 1)
290 gmx_fatal(FARGS, "mts-factor should be larger than 1");
293 // Make the level 0 use the complement of the force groups of group 1
294 mtsLevels[0].forceGroups = ~mtsLevels[1].forceGroups;
295 mtsLevels[0].stepFactor = 1;
297 if ((EEL_FULL(ir.coulombtype) || EVDW_PME(ir.vdwtype))
298 && !mtsLevels[1].forceGroups[static_cast<int>(gmx::MtsForceGroups::LongrangeNonbonded)])
301 "With long-range electrostatics and/or LJ treatment, the long-range part "
302 "has to be part of the mts-level2-forces");
305 if (ir.nstcalcenergy > 0)
307 checkMtsRequirement(ir, "nstcalcenergy", ir.nstcalcenergy, wi);
309 checkMtsRequirement(ir, "nstenergy", ir.nstenergy, wi);
310 checkMtsRequirement(ir, "nstlog", ir.nstlog, wi);
311 if (ir.efep != efepNO)
313 checkMtsRequirement(ir, "nstdhdl", ir.fepvals->nstdhdl, wi);
318 const int pullMtsLevel = gmx::forceGroupMtsLevel(ir.mtsLevels, gmx::MtsForceGroups::Pull);
319 if (ir.pull->nstxout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
321 warning_error(wi, "pull-nstxout should be a multiple of mts-factor");
323 if (ir.pull->nstfout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
325 warning_error(wi, "pull-nstfout should be a multiple of mts-factor");
331 void check_ir(const char* mdparin,
332 const gmx::MdModulesNotifier& mdModulesNotifier,
336 /* Check internal consistency.
337 * NOTE: index groups are not set here yet, don't check things
338 * like temperature coupling group options here, but in triple_check
341 /* Strange macro: first one fills the err_buf, and then one can check
342 * the condition, which will print the message and increase the error
345 #define CHECK(b) _low_check(b, err_buf, wi)
346 char err_buf[256], warn_buf[STRLEN];
349 t_lambda* fep = ir->fepvals;
350 t_expanded* expand = ir->expandedvals;
352 set_warning_line(wi, mdparin, -1);
354 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
356 sprintf(warn_buf, "%s electrostatics is no longer supported", eel_names[eelRF_NEC_UNSUPPORTED]);
357 warning_error(wi, warn_buf);
360 /* BASIC CUT-OFF STUFF */
361 if (ir->rcoulomb < 0)
363 warning_error(wi, "rcoulomb should be >= 0");
367 warning_error(wi, "rvdw should be >= 0");
369 if (ir->rlist < 0 && !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
371 warning_error(wi, "rlist should be >= 0");
374 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
375 "neighbour-list update scheme for efficient buffering for improved energy "
376 "conservation, please use the Verlet cut-off scheme instead.)");
377 CHECK(ir->nstlist < 0);
379 process_interaction_modifier(&ir->coulomb_modifier);
380 process_interaction_modifier(&ir->vdw_modifier);
382 if (ir->cutoff_scheme == ecutsGROUP)
385 "The group cutoff scheme has been removed since GROMACS 2020. "
386 "Please use the Verlet cutoff scheme.");
388 if (ir->cutoff_scheme == ecutsVERLET)
392 /* Normal Verlet type neighbor-list, currently only limited feature support */
393 if (inputrec2nboundeddim(ir) < 3)
395 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
398 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
399 if (ir->rcoulomb != ir->rvdw)
401 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
402 // for PME load balancing, we can support this exception.
403 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == evdwCUT
404 && ir->rcoulomb > ir->rvdw);
405 if (!bUsesPmeTwinRangeKernel)
408 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
409 "rcoulomb>rvdw with PME electrostatics)");
413 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
415 if (ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT)
417 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
420 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
422 evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
423 warning_note(wi, warn_buf);
425 ir->vdwtype = evdwCUT;
429 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
430 evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
431 warning_error(wi, warn_buf);
435 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
438 "With Verlet lists only cut-off and PME LJ interactions are supported");
440 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) || EEL_PME(ir->coulombtype)
441 || ir->coulombtype == eelEWALD))
444 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
445 "electrostatics are supported");
447 if (!(ir->coulomb_modifier == eintmodNONE || ir->coulomb_modifier == eintmodPOTSHIFT))
449 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
450 warning_error(wi, warn_buf);
453 if (EEL_USER(ir->coulombtype))
455 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme",
456 eel_names[ir->coulombtype]);
457 warning_error(wi, warn_buf);
460 if (ir->nstlist <= 0)
462 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
465 if (ir->nstlist < 10)
468 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
469 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
473 rc_max = std::max(ir->rvdw, ir->rcoulomb);
477 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
478 ir->verletbuf_tol = 0;
481 else if (ir->verletbuf_tol <= 0)
483 if (ir->verletbuf_tol == 0)
485 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
488 if (ir->rlist < rc_max)
491 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
494 if (ir->rlist == rc_max && ir->nstlist > 1)
498 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
499 "buffer. The cluster pair list does have a buffering effect, but choosing "
500 "a larger rlist might be necessary for good energy conservation.");
505 if (ir->rlist > rc_max)
508 "You have set rlist larger than the interaction cut-off, but you also "
509 "have verlet-buffer-tolerance > 0. Will set rlist using "
510 "verlet-buffer-tolerance.");
513 if (ir->nstlist == 1)
515 /* No buffer required */
520 if (EI_DYNAMICS(ir->eI))
522 if (inputrec2nboundeddim(ir) < 3)
525 "The box volume is required for calculating rlist from the "
526 "energy drift with verlet-buffer-tolerance > 0. You are "
527 "using at least one unbounded dimension, so no volume can be "
528 "computed. Either use a finite box, or set rlist yourself "
529 "together with verlet-buffer-tolerance = -1.");
531 /* Set rlist temporarily so we can continue processing */
536 /* Set the buffer to 5% of the cut-off */
537 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
543 /* GENERAL INTEGRATOR STUFF */
546 if (ir->etc != etcNO)
548 if (EI_RANDOM(ir->eI))
551 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
552 "implicitly. See the documentation for more information on which "
553 "parameters affect temperature for %s.",
554 etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
559 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
561 etcoupl_names[ir->etc], ei_names[ir->eI]);
563 warning_note(wi, warn_buf);
567 if (ir->eI == eiVVAK)
570 "Integrator method %s is implemented primarily for validation purposes; for "
571 "molecular dynamics, you should probably be using %s or %s",
572 ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
573 warning_note(wi, warn_buf);
575 if (!EI_DYNAMICS(ir->eI))
577 if (ir->epc != epcNO)
580 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
581 epcoupl_names[ir->epc], ei_names[ir->eI]);
582 warning_note(wi, warn_buf);
586 if (EI_DYNAMICS(ir->eI))
588 if (ir->nstcalcenergy < 0)
590 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
591 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
593 /* nstcalcenergy larger than nstener does not make sense.
594 * We ideally want nstcalcenergy=nstener.
598 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
602 ir->nstcalcenergy = ir->nstenergy;
606 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
607 || (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
608 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
611 const char* nsten = "nstenergy";
612 const char* nstdh = "nstdhdl";
613 const char* min_name = nsten;
614 int min_nst = ir->nstenergy;
616 /* find the smallest of ( nstenergy, nstdhdl ) */
617 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
618 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
620 min_nst = ir->fepvals->nstdhdl;
623 /* If the user sets nstenergy small, we should respect that */
624 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy,
626 warning_note(wi, warn_buf);
627 ir->nstcalcenergy = min_nst;
630 if (ir->epc != epcNO)
632 if (ir->nstpcouple < 0)
634 ir->nstpcouple = ir_optimal_nstpcouple(ir);
636 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
639 "With multiple time stepping, nstpcouple should be a mutiple of "
644 if (ir->nstcalcenergy > 0)
646 if (ir->efep != efepNO)
648 /* nstdhdl should be a multiple of nstcalcenergy */
649 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
653 /* nstexpanded should be a multiple of nstcalcenergy */
654 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded",
655 &ir->expandedvals->nstexpanded, wi);
657 /* for storing exact averages nstenergy should be
658 * a multiple of nstcalcenergy
660 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
663 // Inquire all MdModules, if their parameters match with the energy
664 // calculation frequency
665 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
666 mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
668 // Emit all errors from the energy calculation frequency checks
669 for (const std::string& energyFrequencyErrorMessage :
670 energyCalculationFrequencyErrors.errorMessages())
672 warning_error(wi, energyFrequencyErrorMessage);
676 if (ir->nsteps == 0 && !ir->bContinuation)
679 "For a correct single-point energy evaluation with nsteps = 0, use "
680 "continuation = yes to avoid constraining the input coordinates.");
684 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
687 "You are doing a continuation with SD or BD, make sure that ld_seed is "
688 "different from the previous run (using ld_seed=-1 will ensure this)");
694 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
695 CHECK(ir->pbcType != PbcType::Xyz);
696 sprintf(err_buf, "with TPI nstlist should be larger than zero");
697 CHECK(ir->nstlist <= 0);
698 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
699 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
703 if ((opts->nshake > 0) && (opts->bMorse))
705 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
706 warning(wi, warn_buf);
709 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
712 "You are doing a continuation with SD or BD, make sure that ld_seed is "
713 "different from the previous run (using ld_seed=-1 will ensure this)");
715 /* verify simulated tempering options */
719 bool bAllTempZero = TRUE;
720 for (i = 0; i < fep->n_lambda; i++)
722 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
723 efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
724 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
725 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
727 bAllTempZero = FALSE;
730 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
731 CHECK(bAllTempZero == TRUE);
733 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
734 CHECK(ir->eI != eiVV);
736 /* check compatability of the temperature coupling with simulated tempering */
738 if (ir->etc == etcNOSEHOOVER)
741 "Nose-Hoover based temperature control such as [%s] my not be "
742 "entirelyconsistent with simulated tempering",
743 etcoupl_names[ir->etc]);
744 warning_note(wi, warn_buf);
747 /* check that the temperatures make sense */
750 "Higher simulated tempering temperature (%g) must be >= than the simulated "
751 "tempering lower temperature (%g)",
752 ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
753 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
755 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero",
756 ir->simtempvals->simtemp_high);
757 CHECK(ir->simtempvals->simtemp_high <= 0);
759 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero",
760 ir->simtempvals->simtemp_low);
761 CHECK(ir->simtempvals->simtemp_low <= 0);
764 /* verify free energy options */
766 if (ir->efep != efepNO)
769 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
770 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
773 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
775 static_cast<int>(fep->sc_r_power));
776 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
779 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
782 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
784 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
786 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
788 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
789 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
791 sprintf(err_buf, "Free-energy not implemented for Ewald");
792 CHECK(ir->coulombtype == eelEWALD);
794 /* check validty of lambda inputs */
795 if (fep->n_lambda == 0)
797 /* Clear output in case of no states:*/
798 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.",
799 fep->init_fep_state);
800 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
804 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d",
805 fep->init_fep_state, fep->n_lambda - 1);
806 CHECK((fep->init_fep_state >= fep->n_lambda));
810 "Lambda state must be set, either with init-lambda-state or with init-lambda");
811 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
814 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
815 "init-lambda-state or with init-lambda, but not both",
816 fep->init_lambda, fep->init_fep_state);
817 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
820 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
824 for (i = 0; i < efptNR; i++)
826 if (fep->separate_dvdl[i])
831 if (n_lambda_terms > 1)
834 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
835 "use init-lambda to set lambda state (except for slow growth). Use "
836 "init-lambda-state instead.");
837 warning(wi, warn_buf);
840 if (n_lambda_terms < 2 && fep->n_lambda > 0)
843 "init-lambda is deprecated for setting lambda state (except for slow "
844 "growth). Use init-lambda-state instead.");
848 for (j = 0; j < efptNR; j++)
850 for (i = 0; i < fep->n_lambda; i++)
852 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
853 efpt_names[j], fep->all_lambda[j][i]);
854 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
858 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
860 for (i = 0; i < fep->n_lambda; i++)
863 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
864 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
865 "crashes, and is not supported.",
866 i, fep->all_lambda[efptVDW][i], fep->all_lambda[efptCOUL][i]);
867 CHECK((fep->sc_alpha > 0)
868 && (((fep->all_lambda[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL][i] < 1.0))
869 && ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0))));
873 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
875 real sigma, lambda, r_sc;
878 /* Maximum estimate for A and B charges equal with lambda power 1 */
880 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
881 1.0 / fep->sc_r_power);
883 "With PME there is a minor soft core effect present at the cut-off, "
884 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
885 "energy conservation, but usually other effects dominate. With a common sigma "
886 "value of %g nm the fraction of the particle-particle potential at the cut-off "
887 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
888 fep->sc_r_power, sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
889 warning_note(wi, warn_buf);
892 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
893 be treated differently, but that's the next step */
895 for (i = 0; i < efptNR; i++)
897 for (j = 0; j < fep->n_lambda; j++)
899 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
900 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
905 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
909 /* checking equilibration of weights inputs for validity */
912 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
914 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
915 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
918 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
920 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
921 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
924 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
925 expand->equil_steps, elmceq_names[elmceqSTEPS]);
926 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
929 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
930 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
931 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
934 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
935 expand->equil_ratio, elmceq_names[elmceqRATIO]);
936 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
939 "weight-equil-number-all-lambda (%d) must be a positive integer if "
940 "lmc-weights-equil=%s",
941 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
942 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
945 "weight-equil-number-samples (%d) must be a positive integer if "
946 "lmc-weights-equil=%s",
947 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
948 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
951 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
952 expand->equil_steps, elmceq_names[elmceqSTEPS]);
953 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
955 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
956 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
957 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
959 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
960 expand->equil_ratio, elmceq_names[elmceqRATIO]);
961 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
963 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
964 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
965 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
967 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
968 CHECK((expand->lmc_repeats <= 0));
969 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
970 CHECK((expand->minvarmin <= 0));
971 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
972 CHECK((expand->c_range < 0));
974 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
976 fep->init_fep_state, expand->lmc_forced_nstart);
977 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
978 && (expand->elmcmove != elmcmoveNO));
979 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
980 CHECK((expand->lmc_forced_nstart < 0));
981 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
982 fep->init_fep_state);
983 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
985 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
986 CHECK((expand->init_wl_delta < 0));
987 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
988 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
989 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
990 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
992 /* if there is no temperature control, we need to specify an MC temperature */
993 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO)
994 && (expand->mc_temp <= 0.0))
997 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
998 "to a positive number");
999 warning_error(wi, err_buf);
1001 if (expand->nstTij > 0)
1003 sprintf(err_buf, "nstlog must be non-zero");
1004 CHECK(ir->nstlog == 0);
1005 // Avoid modulus by zero in the case that already triggered an error exit.
1006 if (ir->nstlog != 0)
1009 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
1010 expand->nstTij, ir->nstlog);
1011 CHECK((expand->nstTij % ir->nstlog) != 0);
1017 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
1018 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
1021 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
1023 if (ir->pbcType == PbcType::No)
1025 if (ir->epc != epcNO)
1027 warning(wi, "Turning off pressure coupling for vacuum system");
1033 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
1034 c_pbcTypeNames[ir->pbcType].c_str());
1035 CHECK(ir->epc != epcNO);
1037 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
1038 CHECK(EEL_FULL(ir->coulombtype));
1040 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
1041 c_pbcTypeNames[ir->pbcType].c_str());
1042 CHECK(ir->eDispCorr != edispcNO);
1045 if (ir->rlist == 0.0)
1048 "can only have neighborlist cut-off zero (=infinite)\n"
1049 "with coulombtype = %s or coulombtype = %s\n"
1050 "without periodic boundary conditions (pbc = %s) and\n"
1051 "rcoulomb and rvdw set to zero",
1052 eel_names[eelCUT], eel_names[eelUSER], c_pbcTypeNames[PbcType::No].c_str());
1053 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
1054 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1056 if (ir->nstlist > 0)
1059 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1060 "nstype=simple and only one MPI rank");
1065 if (ir->nstcomm == 0)
1067 // TODO Change this behaviour. There should be exactly one way
1068 // to turn off an algorithm.
1069 ir->comm_mode = ecmNO;
1071 if (ir->comm_mode != ecmNO)
1073 if (ir->nstcomm < 0)
1075 // TODO Such input was once valid. Now that we've been
1076 // helpful for a few years, we should reject such input,
1077 // lest we have to support every historical decision
1080 "If you want to remove the rotation around the center of mass, you should set "
1081 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1082 "its absolute value");
1083 ir->nstcomm = abs(ir->nstcomm);
1086 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
1089 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
1090 "nstcomm to nstcalcenergy");
1091 ir->nstcomm = ir->nstcalcenergy;
1094 if (ir->comm_mode == ecmANGULAR)
1097 "Can not remove the rotation around the center of mass with periodic "
1099 CHECK(ir->bPeriodicMols);
1100 if (ir->pbcType != PbcType::No)
1103 "Removing the rotation around the center of mass in a periodic system, "
1104 "this can lead to artifacts. Only use this on a single (cluster of) "
1105 "molecules. This cluster should not cross periodic boundaries.");
1110 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
1113 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1114 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1117 warning_note(wi, warn_buf);
1120 /* TEMPERATURE COUPLING */
1121 if (ir->etc == etcYES)
1123 ir->etc = etcBERENDSEN;
1125 "Old option for temperature coupling given: "
1126 "changing \"yes\" to \"Berendsen\"\n");
1129 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
1131 if (ir->opts.nhchainlength < 1)
1134 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1136 ir->opts.nhchainlength);
1137 ir->opts.nhchainlength = 1;
1138 warning(wi, warn_buf);
1141 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1145 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1146 ir->opts.nhchainlength = 1;
1151 ir->opts.nhchainlength = 0;
1154 if (ir->eI == eiVVAK)
1157 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1160 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1163 if (ETC_ANDERSEN(ir->etc))
1165 sprintf(err_buf, "%s temperature control not supported for integrator %s.",
1166 etcoupl_names[ir->etc], ei_names[ir->eI]);
1167 CHECK(!(EI_VV(ir->eI)));
1169 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
1172 "Center of mass removal not necessary for %s. All velocities of coupled "
1173 "groups are rerandomized periodically, so flying ice cube errors will not "
1175 etcoupl_names[ir->etc]);
1176 warning_note(wi, warn_buf);
1180 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1181 "randomized every time step",
1182 ir->nstcomm, etcoupl_names[ir->etc]);
1183 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
1186 if (ir->etc == etcBERENDSEN)
1189 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1190 "might want to consider using the %s thermostat.",
1191 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
1192 warning_note(wi, warn_buf);
1195 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc)) && ir->epc == epcBERENDSEN)
1198 "Using Berendsen pressure coupling invalidates the "
1199 "true ensemble for the thermostat");
1200 warning(wi, warn_buf);
1203 /* PRESSURE COUPLING */
1204 if (ir->epc == epcISOTROPIC)
1206 ir->epc = epcBERENDSEN;
1208 "Old option for pressure coupling given: "
1209 "changing \"Isotropic\" to \"Berendsen\"\n");
1212 if (ir->epc != epcNO)
1214 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1216 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1217 CHECK(ir->tau_p <= 0);
1219 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1222 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1223 "times larger than nstpcouple*dt (%g)",
1224 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1225 warning(wi, warn_buf);
1229 "compressibility must be > 0 when using pressure"
1231 EPCOUPLTYPE(ir->epc));
1232 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1233 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1234 && ir->compress[ZZ][YY] <= 0));
1236 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1239 "You are generating velocities so I am assuming you "
1240 "are equilibrating a system. You are using "
1241 "%s pressure coupling, but this can be "
1242 "unstable for equilibration. If your system crashes, try "
1243 "equilibrating first with Berendsen pressure coupling. If "
1244 "you are not equilibrating the system, you can probably "
1245 "ignore this warning.",
1246 epcoupl_names[ir->epc]);
1247 warning(wi, warn_buf);
1253 if (ir->epc == epcMTTK)
1255 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1259 /* ELECTROSTATICS */
1260 /* More checks are in triple check (grompp.c) */
1262 if (ir->coulombtype == eelSWITCH)
1265 "coulombtype = %s is only for testing purposes and can lead to serious "
1266 "artifacts, advice: use coulombtype = %s",
1267 eel_names[ir->coulombtype], eel_names[eelRF_ZERO]);
1268 warning(wi, warn_buf);
1271 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1274 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1275 "format and exchanging epsilon-r and epsilon-rf",
1277 warning(wi, warn_buf);
1278 ir->epsilon_rf = ir->epsilon_r;
1279 ir->epsilon_r = 1.0;
1282 if (ir->epsilon_r == 0)
1285 "It is pointless to use long-range electrostatics with infinite relative "
1287 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1289 CHECK(EEL_FULL(ir->coulombtype));
1292 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1294 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1295 CHECK(ir->epsilon_r < 0);
1298 if (EEL_RF(ir->coulombtype))
1300 /* reaction field (at the cut-off) */
1302 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1305 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1306 eel_names[ir->coulombtype]);
1307 warning(wi, warn_buf);
1308 ir->epsilon_rf = 0.0;
1311 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1312 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1313 if (ir->epsilon_rf == ir->epsilon_r)
1315 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1316 eel_names[ir->coulombtype]);
1317 warning(wi, warn_buf);
1320 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1321 * means the interaction is zero outside rcoulomb, but it helps to
1322 * provide accurate energy conservation.
1324 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1326 if (ir_coulomb_switched(ir))
1329 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1330 "potential modifier options!",
1331 eel_names[ir->coulombtype]);
1332 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1336 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1339 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1340 "secondary coulomb-modifier.");
1341 CHECK(ir->coulomb_modifier != eintmodNONE);
1343 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1346 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1347 "secondary vdw-modifier.");
1348 CHECK(ir->vdw_modifier != eintmodNONE);
1351 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->vdwtype == evdwSWITCH
1352 || ir->vdwtype == evdwSHIFT)
1355 "The switch/shift interaction settings are just for compatibility; you will get "
1357 "performance from applying potential modifiers to your interactions!\n");
1358 warning_note(wi, warn_buf);
1361 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1363 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1365 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1367 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1368 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1369 "will be good regardless, since ewald_rtol = %g.",
1370 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1371 warning(wi, warn_buf);
1375 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1377 if (ir->rvdw_switch == 0)
1380 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1381 "potential. This suggests it was not set in the mdp, which can lead to large "
1382 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1383 "switching range.");
1384 warning(wi, warn_buf);
1388 if (EEL_FULL(ir->coulombtype))
1390 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER
1391 || ir->coulombtype == eelPMEUSERSWITCH)
1393 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1394 eel_names[ir->coulombtype]);
1395 CHECK(ir->rcoulomb > ir->rlist);
1399 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1401 // TODO: Move these checks into the ewald module with the options class
1403 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1405 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1407 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d",
1408 eel_names[ir->coulombtype], orderMin, orderMax);
1409 warning_error(wi, warn_buf);
1413 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1415 if (ir->ewald_geometry == eewg3D)
1417 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1418 c_pbcTypeNames[ir->pbcType].c_str(), eewg_names[eewg3DC]);
1419 warning(wi, warn_buf);
1421 /* This check avoids extra pbc coding for exclusion corrections */
1422 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1423 CHECK(ir->wall_ewald_zfac < 2);
1425 if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
1427 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1428 eel_names[ir->coulombtype], eewg_names[eewg3DC], c_pbcTypeNames[PbcType::XY].c_str());
1429 warning(wi, warn_buf);
1431 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1433 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1434 CHECK(ir->bPeriodicMols);
1435 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1436 warning_note(wi, warn_buf);
1438 "With epsilon_surface > 0 you can only use domain decomposition "
1439 "when there are only small molecules with all bonds constrained (mdrun will check "
1441 warning_note(wi, warn_buf);
1444 if (ir_vdw_switched(ir))
1446 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1447 CHECK(ir->rvdw_switch >= ir->rvdw);
1449 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1452 "You are applying a switch function to vdw forces or potentials from %g to %g "
1453 "nm, which is more than half the interaction range, whereas switch functions "
1454 "are intended to act only close to the cut-off.",
1455 ir->rvdw_switch, ir->rvdw);
1456 warning_note(wi, warn_buf);
1460 if (ir->vdwtype == evdwPME)
1462 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1464 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1465 evdw_names[ir->vdwtype], eintmod_names[eintmodPOTSHIFT], eintmod_names[eintmodNONE]);
1466 warning_error(wi, err_buf);
1470 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1473 "You have selected user tables with dispersion correction, the dispersion "
1474 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1475 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1476 "really want dispersion correction to -C6/r^6.");
1479 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT) && ir->rvdw != 0)
1481 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1484 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1486 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1489 /* IMPLICIT SOLVENT */
1490 if (ir->coulombtype == eelGB_NOTUSED)
1492 sprintf(warn_buf, "Invalid option %s for coulombtype", eel_names[ir->coulombtype]);
1493 warning_error(wi, warn_buf);
1498 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1503 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1507 /* interpret a number of doubles from a string and put them in an array,
1508 after allocating space for them.
1509 str = the input string
1510 n = the (pre-allocated) number of doubles read
1511 r = the output array of doubles. */
1512 static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
1514 auto values = gmx::splitString(str);
1518 for (int i = 0; i < *n; i++)
1522 (*r)[i] = gmx::fromString<real>(values[i]);
1524 catch (gmx::GromacsException&)
1526 warning_error(wi, "Invalid value " + values[i]
1527 + " in string in mdp file. Expected a real number.");
1533 static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1536 int i, j, max_n_lambda, nweights, nfep[efptNR];
1537 t_lambda* fep = ir->fepvals;
1538 t_expanded* expand = ir->expandedvals;
1539 real** count_fep_lambdas;
1540 bool bOneLambda = TRUE;
1542 snew(count_fep_lambdas, efptNR);
1544 /* FEP input processing */
1545 /* first, identify the number of lambda values for each type.
1546 All that are nonzero must have the same number */
1548 for (i = 0; i < efptNR; i++)
1550 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1553 /* now, determine the number of components. All must be either zero, or equal. */
1556 for (i = 0; i < efptNR; i++)
1558 if (nfep[i] > max_n_lambda)
1560 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1561 must have the same number if its not zero.*/
1566 for (i = 0; i < efptNR; i++)
1570 ir->fepvals->separate_dvdl[i] = FALSE;
1572 else if (nfep[i] == max_n_lambda)
1574 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute
1575 the derivative with respect to the temperature currently */
1577 ir->fepvals->separate_dvdl[i] = TRUE;
1583 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1585 nfep[i], efpt_names[i], max_n_lambda);
1588 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1589 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1591 /* the number of lambdas is the number we've read in, which is either zero
1592 or the same for all */
1593 fep->n_lambda = max_n_lambda;
1595 /* allocate space for the array of lambda values */
1596 snew(fep->all_lambda, efptNR);
1597 /* if init_lambda is defined, we need to set lambda */
1598 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1600 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1602 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1603 for (i = 0; i < efptNR; i++)
1605 snew(fep->all_lambda[i], fep->n_lambda);
1606 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1609 for (j = 0; j < fep->n_lambda; j++)
1611 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1613 sfree(count_fep_lambdas[i]);
1616 sfree(count_fep_lambdas);
1618 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1619 internal bookkeeping -- for now, init_lambda */
1621 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1623 for (i = 0; i < fep->n_lambda; i++)
1625 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1629 /* check to see if only a single component lambda is defined, and soft core is defined.
1630 In this case, turn on coulomb soft core */
1632 if (max_n_lambda == 0)
1638 for (i = 0; i < efptNR; i++)
1640 if ((nfep[i] != 0) && (i != efptFEP))
1646 if ((bOneLambda) && (fep->sc_alpha > 0))
1648 fep->bScCoul = TRUE;
1651 /* Fill in the others with the efptFEP if they are not explicitly
1652 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1653 they are all zero. */
1655 for (i = 0; i < efptNR; i++)
1657 if ((nfep[i] == 0) && (i != efptFEP))
1659 for (j = 0; j < fep->n_lambda; j++)
1661 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1667 /* now read in the weights */
1668 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1671 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1673 else if (nweights != fep->n_lambda)
1675 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1676 nweights, fep->n_lambda);
1678 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1680 expand->nstexpanded = fep->nstdhdl;
1681 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1686 static void do_simtemp_params(t_inputrec* ir)
1689 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1690 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1693 template<typename T>
1694 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1697 for (const auto& input : inputs)
1701 outputs[i] = gmx::fromStdString<T>(input);
1703 catch (gmx::GromacsException&)
1705 auto message = gmx::formatString(
1706 "Invalid value for mdp option %s. %s should only consist of integers separated "
1709 warning_error(wi, message);
1715 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1718 for (const auto& input : inputs)
1722 outputs[i] = gmx::fromString<real>(input);
1724 catch (gmx::GromacsException&)
1726 auto message = gmx::formatString(
1727 "Invalid value for mdp option %s. %s should only consist of real numbers "
1728 "separated by spaces.",
1730 warning_error(wi, message);
1736 static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
1739 for (const auto& input : inputs)
1743 outputs[i][d] = gmx::fromString<real>(input);
1745 catch (gmx::GromacsException&)
1747 auto message = gmx::formatString(
1748 "Invalid value for mdp option %s. %s should only consist of real numbers "
1749 "separated by spaces.",
1751 warning_error(wi, message);
1762 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1764 opts->wall_atomtype[0] = nullptr;
1765 opts->wall_atomtype[1] = nullptr;
1767 ir->wall_atomtype[0] = -1;
1768 ir->wall_atomtype[1] = -1;
1769 ir->wall_density[0] = 0;
1770 ir->wall_density[1] = 0;
1774 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1775 if (wallAtomTypes.size() != size_t(ir->nwall))
1777 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu", ir->nwall,
1778 wallAtomTypes.size());
1780 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1781 for (int i = 0; i < ir->nwall; i++)
1783 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1786 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1788 auto wallDensity = gmx::splitString(wall_density);
1789 if (wallDensity.size() != size_t(ir->nwall))
1791 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall,
1792 wallDensity.size());
1794 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1795 for (int i = 0; i < ir->nwall; i++)
1797 if (ir->wall_density[i] <= 0)
1799 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1806 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1810 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1811 for (int i = 0; i < nwall; i++)
1813 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1814 grps->emplace_back(groups->groupNames.size() - 1);
1819 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1821 /* read expanded ensemble parameters */
1822 printStringNewline(inp, "expanded ensemble variables");
1823 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1824 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1825 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1826 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1827 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1828 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1829 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1830 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1831 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1832 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1833 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1834 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1835 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1836 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1837 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1838 expand->bSymmetrizedTMatrix =
1839 (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1840 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1841 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1842 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1843 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1844 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1845 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1846 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1849 /*! \brief Return whether an end state with the given coupling-lambda
1850 * value describes fully-interacting VDW.
1852 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1853 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1855 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1857 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1863 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1866 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1868 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1870 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1873 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1874 std::string message = gmx::formatExceptionMessageToString(*ex);
1875 warning_error(wi_, message.c_str());
1880 std::string getOptionName(const gmx::KeyValueTreePath& context)
1882 if (mapping_ != nullptr)
1884 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1885 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1888 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1893 const gmx::IKeyValueTreeBackMapping* mapping_;
1898 void get_ir(const char* mdparin,
1899 const char* mdparout,
1900 gmx::MDModules* mdModules,
1903 WriteMdpHeader writeMdpHeader,
1907 double dumdub[2][6];
1909 char warn_buf[STRLEN];
1910 t_lambda* fep = ir->fepvals;
1911 t_expanded* expand = ir->expandedvals;
1913 const char* no_names[] = { "no", nullptr };
1915 init_inputrec_strings();
1916 gmx::TextInputFile stream(mdparin);
1917 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1919 snew(dumstr[0], STRLEN);
1920 snew(dumstr[1], STRLEN);
1922 /* ignore the following deprecated commands */
1923 replace_inp_entry(inp, "title", nullptr);
1924 replace_inp_entry(inp, "cpp", nullptr);
1925 replace_inp_entry(inp, "domain-decomposition", nullptr);
1926 replace_inp_entry(inp, "andersen-seed", nullptr);
1927 replace_inp_entry(inp, "dihre", nullptr);
1928 replace_inp_entry(inp, "dihre-fc", nullptr);
1929 replace_inp_entry(inp, "dihre-tau", nullptr);
1930 replace_inp_entry(inp, "nstdihreout", nullptr);
1931 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1932 replace_inp_entry(inp, "optimize-fft", nullptr);
1933 replace_inp_entry(inp, "adress_type", nullptr);
1934 replace_inp_entry(inp, "adress_const_wf", nullptr);
1935 replace_inp_entry(inp, "adress_ex_width", nullptr);
1936 replace_inp_entry(inp, "adress_hy_width", nullptr);
1937 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1938 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1939 replace_inp_entry(inp, "adress_site", nullptr);
1940 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1941 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1942 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1943 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1944 replace_inp_entry(inp, "rlistlong", nullptr);
1945 replace_inp_entry(inp, "nstcalclr", nullptr);
1946 replace_inp_entry(inp, "pull-print-com2", nullptr);
1947 replace_inp_entry(inp, "gb-algorithm", nullptr);
1948 replace_inp_entry(inp, "nstgbradii", nullptr);
1949 replace_inp_entry(inp, "rgbradii", nullptr);
1950 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1951 replace_inp_entry(inp, "gb-saltconc", nullptr);
1952 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1953 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1954 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1955 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1956 replace_inp_entry(inp, "sa-algorithm", nullptr);
1957 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1958 replace_inp_entry(inp, "ns-type", nullptr);
1960 /* replace the following commands with the clearer new versions*/
1961 replace_inp_entry(inp, "unconstrained-start", "continuation");
1962 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1963 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1964 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1965 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1966 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1967 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1969 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1970 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1971 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1972 setStringEntry(&inp, "include", opts->include, nullptr);
1973 printStringNoNewline(
1974 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1975 setStringEntry(&inp, "define", opts->define, nullptr);
1977 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1978 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1979 printStringNoNewline(&inp, "Start time and timestep in ps");
1980 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1981 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1982 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1983 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1984 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1985 printStringNoNewline(
1986 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1987 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1988 printStringNoNewline(&inp, "Multiple time-stepping");
1989 ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
1992 opts->numMtsLevels = get_eint(&inp, "mts-levels", 2, wi);
1993 ir->mtsLevels.resize(2);
1994 gmx::MtsLevel& mtsLevel = ir->mtsLevels[1];
1995 opts->mtsLevel2Forces = setStringEntry(&inp, "mts-level2-forces",
1996 "longrange-nonbonded nonbonded pair dihedral");
1997 mtsLevel.stepFactor = get_eint(&inp, "mts-level2-factor", 2, wi);
1999 // We clear after reading without dynamics to not force the user to remove MTS mdp options
2000 if (!EI_DYNAMICS(ir->eI))
2003 ir->mtsLevels.clear();
2006 printStringNoNewline(&inp, "mode for center of mass motion removal");
2007 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
2008 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
2009 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
2010 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
2011 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
2013 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2014 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2015 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2016 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2019 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2020 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2021 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2022 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2023 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2024 ir->niter = get_eint(&inp, "niter", 20, wi);
2025 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2026 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2027 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2028 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2029 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2031 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2032 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2034 /* Output options */
2035 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2036 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2037 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2038 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2039 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2040 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2041 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2042 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2043 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2044 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2045 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2046 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2047 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2048 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2049 printStringNoNewline(&inp, "default, all atoms will be written.");
2050 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2051 printStringNoNewline(&inp, "Selection of energy groups");
2052 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2054 /* Neighbor searching */
2055 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2056 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2057 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
2058 printStringNoNewline(&inp, "nblist update frequency");
2059 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2060 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2061 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2062 std::vector<const char*> pbcTypesNamesChar;
2063 for (const auto& pbcTypeName : c_pbcTypeNames)
2065 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2067 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2068 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
2069 printStringNoNewline(&inp,
2070 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2071 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2072 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2073 printStringNoNewline(&inp, "nblist cut-off");
2074 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2075 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2077 /* Electrostatics */
2078 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2079 printStringNoNewline(&inp, "Method for doing electrostatics");
2080 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
2081 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
2082 printStringNoNewline(&inp, "cut-off lengths");
2083 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2084 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2085 printStringNoNewline(&inp,
2086 "Relative dielectric constant for the medium and the reaction field");
2087 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2088 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2089 printStringNoNewline(&inp, "Method for doing Van der Waals");
2090 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
2091 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
2092 printStringNoNewline(&inp, "cut-off lengths");
2093 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2094 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2095 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2096 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
2097 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2098 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2099 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2100 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2101 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2102 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2103 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2104 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2105 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2106 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2107 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2108 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2109 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2110 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2111 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
2112 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
2113 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2115 /* Implicit solvation is no longer supported, but we need grompp
2116 to be able to refuse old .mdp files that would have built a tpr
2117 to run it. Thus, only "no" is accepted. */
2118 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2120 /* Coupling stuff */
2121 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2122 printStringNoNewline(&inp, "Temperature coupling");
2123 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
2124 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2125 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2126 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
2127 printStringNoNewline(&inp, "Groups to couple separately");
2128 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2129 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2130 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2131 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2132 printStringNoNewline(&inp, "pressure coupling");
2133 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
2134 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
2135 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2136 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2137 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2138 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2139 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2140 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2141 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
2144 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2145 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
2146 printStringNoNewline(&inp, "Groups treated with MiMiC");
2147 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2149 /* Simulated annealing */
2150 printStringNewline(&inp, "SIMULATED ANNEALING");
2151 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2152 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2153 printStringNoNewline(&inp,
2154 "Number of time points to use for specifying annealing in each group");
2155 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2156 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2157 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2158 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2159 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2162 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2163 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
2164 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2165 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2168 printStringNewline(&inp, "OPTIONS FOR BONDS");
2169 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2170 printStringNoNewline(&inp, "Type of constraint algorithm");
2171 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
2172 printStringNoNewline(&inp, "Do not constrain the start configuration");
2173 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
2174 printStringNoNewline(&inp,
2175 "Use successive overrelaxation to reduce the number of shake iterations");
2176 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
2177 printStringNoNewline(&inp, "Relative tolerance of shake");
2178 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2179 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2180 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2181 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2182 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2183 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2184 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2185 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2186 printStringNoNewline(&inp, "rotates over more degrees than");
2187 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2188 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2189 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
2191 /* Energy group exclusions */
2192 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2193 printStringNoNewline(
2194 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2195 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2198 printStringNewline(&inp, "WALLS");
2199 printStringNoNewline(
2200 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2201 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2202 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
2203 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2204 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2205 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2206 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2209 printStringNewline(&inp, "COM PULLING");
2210 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2213 ir->pull = std::make_unique<pull_params_t>();
2214 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
2218 for (int c = 0; c < ir->pull->ncoord; c++)
2220 if (ir->pull->coord[c].eType == epullCONSTRAINT)
2223 "Constraint COM pulling is not supported in combination with "
2224 "multiple time stepping");
2232 NOTE: needs COM pulling or free energy input */
2233 printStringNewline(&inp, "AWH biasing");
2234 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2237 ir->awhParams = gmx::readAwhParams(&inp, wi);
2240 /* Enforced rotation */
2241 printStringNewline(&inp, "ENFORCED ROTATION");
2242 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2243 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2247 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2250 /* Interactive MD */
2252 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2253 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2254 if (inputrecStrings->imd_grp[0] != '\0')
2261 printStringNewline(&inp, "NMR refinement stuff");
2262 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2263 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2264 printStringNoNewline(
2265 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2266 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2267 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2268 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2269 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2270 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2271 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2272 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2273 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2274 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2275 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2276 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2277 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2278 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2279 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2280 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2282 /* free energy variables */
2283 printStringNewline(&inp, "Free energy variables");
2284 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2285 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2286 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2287 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2288 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2290 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2292 it was not entered */
2293 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2294 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2295 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2296 setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
2297 setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
2298 setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
2299 setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
2300 setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
2301 setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
2302 setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
2303 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2304 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2305 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2306 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2307 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2308 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2309 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2310 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2311 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2312 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2313 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2314 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2315 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2316 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2318 /* Non-equilibrium MD stuff */
2319 printStringNewline(&inp, "Non-equilibrium MD stuff");
2320 setStringEntry(&inp, "acc-grps", inputrecStrings->accgrps, nullptr);
2321 setStringEntry(&inp, "accelerate", inputrecStrings->acc, nullptr);
2322 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2323 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2324 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2325 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2327 /* simulated tempering variables */
2328 printStringNewline(&inp, "simulated tempering variables");
2329 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2330 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2331 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2332 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2334 /* expanded ensemble variables */
2335 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2337 read_expandedparams(&inp, expand, wi);
2340 /* Electric fields */
2342 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2343 gmx::KeyValueTreeTransformer transform;
2344 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2345 mdModules->initMdpTransform(transform.rules());
2346 for (const auto& path : transform.mappedPaths())
2348 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2349 mark_einp_set(inp, path[0].c_str());
2351 MdpErrorHandler errorHandler(wi);
2352 auto result = transform.transform(convertedValues, &errorHandler);
2353 ir->params = new gmx::KeyValueTreeObject(result.object());
2354 mdModules->adjustInputrecBasedOnModules(ir);
2355 errorHandler.setBackMapping(result.backMapping());
2356 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2359 /* Ion/water position swapping ("computational electrophysiology") */
2360 printStringNewline(&inp,
2361 "Ion/water position swapping for computational electrophysiology setups");
2362 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2363 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2364 if (ir->eSwapCoords != eswapNO)
2371 printStringNoNewline(&inp, "Swap attempt frequency");
2372 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2373 printStringNoNewline(&inp, "Number of ion types to be controlled");
2374 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2377 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2379 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2380 snew(ir->swap->grp, ir->swap->ngrp);
2381 for (i = 0; i < ir->swap->ngrp; i++)
2383 snew(ir->swap->grp[i].molname, STRLEN);
2385 printStringNoNewline(&inp,
2386 "Two index groups that contain the compartment-partitioning atoms");
2387 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2388 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2389 printStringNoNewline(&inp,
2390 "Use center of mass of split groups (yes/no), otherwise center of "
2391 "geometry is used");
2392 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2393 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2395 printStringNoNewline(&inp, "Name of solvent molecules");
2396 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2398 printStringNoNewline(&inp,
2399 "Split cylinder: radius, upper and lower extension (nm) (this will "
2400 "define the channels)");
2401 printStringNoNewline(&inp,
2402 "Note that the split cylinder settings do not have an influence on "
2403 "the swapping protocol,");
2404 printStringNoNewline(
2406 "however, if correctly defined, the permeation events are recorded per channel");
2407 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2408 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2409 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2410 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2411 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2412 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2414 printStringNoNewline(
2416 "Average the number of ions per compartment over these many swap attempt steps");
2417 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2419 printStringNoNewline(
2420 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2421 printStringNoNewline(
2422 &inp, "and the requested number of ions of this type in compartments A and B");
2423 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2424 for (i = 0; i < nIonTypes; i++)
2426 int ig = eSwapFixedGrpNR + i;
2428 sprintf(buf, "iontype%d-name", i);
2429 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2430 sprintf(buf, "iontype%d-in-A", i);
2431 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2432 sprintf(buf, "iontype%d-in-B", i);
2433 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2436 printStringNoNewline(
2438 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2439 printStringNoNewline(
2441 "at maximum distance (= bulk concentration) to the split group layers. However,");
2442 printStringNoNewline(&inp,
2443 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2444 "layer from the middle at 0.0");
2445 printStringNoNewline(&inp,
2446 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2447 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2448 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2449 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2450 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2452 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2455 printStringNoNewline(
2456 &inp, "Start to swap ions if threshold difference to requested count is reached");
2457 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2460 /* AdResS is no longer supported, but we need grompp to be able to
2461 refuse to process old .mdp files that used it. */
2462 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2464 /* User defined thingies */
2465 printStringNewline(&inp, "User defined thingies");
2466 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2467 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2468 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2469 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2470 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2471 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2472 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2473 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2474 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2475 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2479 gmx::TextOutputFile stream(mdparout);
2480 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2482 // Transform module data into a flat key-value tree for output.
2483 gmx::KeyValueTreeBuilder builder;
2484 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2485 mdModules->buildMdpOutput(&builderObject);
2487 gmx::TextWriter writer(&stream);
2488 writeKeyValueTreeAsMdp(&writer, builder.build());
2493 /* Process options if necessary */
2494 for (m = 0; m < 2; m++)
2496 for (i = 0; i < 2 * DIM; i++)
2505 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2509 "Pressure coupling incorrect number of values (I need exactly 1)");
2511 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2513 case epctSEMIISOTROPIC:
2514 case epctSURFACETENSION:
2515 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2519 "Pressure coupling incorrect number of values (I need exactly 2)");
2521 dumdub[m][YY] = dumdub[m][XX];
2523 case epctANISOTROPIC:
2524 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf", &(dumdub[m][XX]), &(dumdub[m][YY]),
2525 &(dumdub[m][ZZ]), &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5]))
2530 "Pressure coupling incorrect number of values (I need exactly 6)");
2534 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2535 epcoupltype_names[ir->epct]);
2539 clear_mat(ir->ref_p);
2540 clear_mat(ir->compress);
2541 for (i = 0; i < DIM; i++)
2543 ir->ref_p[i][i] = dumdub[1][i];
2544 ir->compress[i][i] = dumdub[0][i];
2546 if (ir->epct == epctANISOTROPIC)
2548 ir->ref_p[XX][YY] = dumdub[1][3];
2549 ir->ref_p[XX][ZZ] = dumdub[1][4];
2550 ir->ref_p[YY][ZZ] = dumdub[1][5];
2551 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2554 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2555 "apply a threefold shear stress?\n");
2557 ir->compress[XX][YY] = dumdub[0][3];
2558 ir->compress[XX][ZZ] = dumdub[0][4];
2559 ir->compress[YY][ZZ] = dumdub[0][5];
2560 for (i = 0; i < DIM; i++)
2562 for (m = 0; m < i; m++)
2564 ir->ref_p[i][m] = ir->ref_p[m][i];
2565 ir->compress[i][m] = ir->compress[m][i];
2570 if (ir->comm_mode == ecmNO)
2575 opts->couple_moltype = nullptr;
2576 if (strlen(inputrecStrings->couple_moltype) > 0)
2578 if (ir->efep != efepNO)
2580 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2581 if (opts->couple_lam0 == opts->couple_lam1)
2583 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2585 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2589 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2596 "Free energy is turned off, so we will not decouple the molecule listed "
2600 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2601 if (ir->efep != efepNO)
2603 if (fep->delta_lambda != 0)
2605 ir->efep = efepSLOWGROWTH;
2609 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2611 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2613 "Old option for dhdl-print-energy given: "
2614 "changing \"yes\" to \"total\"\n");
2617 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2619 /* always print out the energy to dhdl if we are doing
2620 expanded ensemble, since we need the total energy for
2621 analysis if the temperature is changing. In some
2622 conditions one may only want the potential energy, so
2623 we will allow that if the appropriate mdp setting has
2624 been enabled. Otherwise, total it is:
2626 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2629 if ((ir->efep != efepNO) || ir->bSimTemp)
2631 ir->bExpanded = FALSE;
2632 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2634 ir->bExpanded = TRUE;
2636 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2637 if (ir->bSimTemp) /* done after fep params */
2639 do_simtemp_params(ir);
2642 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2643 * setup and not on the old way of specifying the free-energy setup,
2644 * we should check for using soft-core when not needed, since that
2645 * can complicate the sampling significantly.
2646 * Note that we only check for the automated coupling setup.
2647 * If the (advanced) user does FEP through manual topology changes,
2648 * this check will not be triggered.
2650 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 && ir->fepvals->sc_alpha != 0
2651 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2654 "You are using soft-core interactions while the Van der Waals interactions are "
2655 "not decoupled (note that the sc-coul option is only active when using lambda "
2656 "states). Although this will not lead to errors, you will need much more "
2657 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2662 ir->fepvals->n_lambda = 0;
2665 /* WALL PARAMETERS */
2667 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2669 /* ORIENTATION RESTRAINT PARAMETERS */
2671 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2673 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2676 /* DEFORMATION PARAMETERS */
2678 clear_mat(ir->deform);
2679 for (i = 0; i < 6; i++)
2684 double gmx_unused canary;
2685 int ndeform = sscanf(inputrecStrings->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]),
2686 &(dumdub[0][1]), &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]),
2687 &(dumdub[0][5]), &canary);
2689 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2693 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2694 inputrecStrings->deform)
2697 for (i = 0; i < 3; i++)
2699 ir->deform[i][i] = dumdub[0][i];
2701 ir->deform[YY][XX] = dumdub[0][3];
2702 ir->deform[ZZ][XX] = dumdub[0][4];
2703 ir->deform[ZZ][YY] = dumdub[0][5];
2704 if (ir->epc != epcNO)
2706 for (i = 0; i < 3; i++)
2708 for (j = 0; j <= i; j++)
2710 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2712 warning_error(wi, "A box element has deform set and compressibility > 0");
2716 for (i = 0; i < 3; i++)
2718 for (j = 0; j < i; j++)
2720 if (ir->deform[i][j] != 0)
2722 for (m = j; m < DIM; m++)
2724 if (ir->compress[m][j] != 0)
2727 "An off-diagonal box element has deform set while "
2728 "compressibility > 0 for the same component of another box "
2729 "vector, this might lead to spurious periodicity effects.");
2730 warning(wi, warn_buf);
2738 /* Ion/water position swapping checks */
2739 if (ir->eSwapCoords != eswapNO)
2741 if (ir->swap->nstswap < 1)
2743 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2745 if (ir->swap->nAverage < 1)
2747 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2749 if (ir->swap->threshold < 1.0)
2751 warning_error(wi, "Ion count threshold must be at least 1.\n");
2755 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2758 setupMtsLevels(ir->mtsLevels, *ir, *opts, wi);
2763 gmx::checkAwhParams(ir->awhParams, ir, wi);
2770 /* We would like gn to be const as well, but C doesn't allow this */
2771 /* TODO this is utility functionality (search for the index of a
2772 string in a collection), so should be refactored and located more
2774 int search_string(const char* s, int ng, char* gn[])
2778 for (i = 0; (i < ng); i++)
2780 if (gmx_strcasecmp(s, gn[i]) == 0)
2787 "Group %s referenced in the .mdp file was not found in the index file.\n"
2788 "Group names must match either [moleculetype] names or custom index group\n"
2789 "names, in which case you must supply an index file to the '-n' option\n"
2794 static void do_numbering(int natoms,
2795 SimulationGroups* groups,
2796 gmx::ArrayRef<std::string> groupsFromMdpFile,
2799 SimulationAtomGroupType gtype,
2805 unsigned short* cbuf;
2806 AtomGroupIndices* grps = &(groups->groups[gtype]);
2807 int j, gid, aj, ognr, ntot = 0;
2809 char warn_buf[STRLEN];
2811 title = shortName(gtype);
2814 /* Mark all id's as not set */
2815 for (int i = 0; (i < natoms); i++)
2820 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2822 /* Lookup the group name in the block structure */
2823 gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2824 if ((grptp != egrptpONE) || (i == 0))
2826 grps->emplace_back(gid);
2829 /* Now go over the atoms in the group */
2830 for (j = block->index[gid]; (j < block->index[gid + 1]); j++)
2835 /* Range checking */
2836 if ((aj < 0) || (aj >= natoms))
2838 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2840 /* Lookup up the old group number */
2844 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
2849 /* Store the group number in buffer */
2850 if (grptp == egrptpONE)
2863 /* Now check whether we have done all atoms */
2866 if (grptp == egrptpALL)
2868 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2870 else if (grptp == egrptpPART)
2872 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2873 warning_note(wi, warn_buf);
2875 /* Assign all atoms currently unassigned to a rest group */
2876 for (j = 0; (j < natoms); j++)
2878 if (cbuf[j] == NOGID)
2880 cbuf[j] = grps->size();
2883 if (grptp != egrptpPART)
2887 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title,
2890 /* Add group name "rest" */
2891 grps->emplace_back(restnm);
2893 /* Assign the rest name to all atoms not currently assigned to a group */
2894 for (j = 0; (j < natoms); j++)
2896 if (cbuf[j] == NOGID)
2898 // group size was not updated before this here, so need to use -1.
2899 cbuf[j] = grps->size() - 1;
2905 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2907 /* All atoms are part of one (or no) group, no index required */
2908 groups->groupNumbers[gtype].clear();
2912 for (int j = 0; (j < natoms); j++)
2914 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2921 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2924 pull_params_t* pull;
2925 int natoms, imin, jmin;
2926 int * nrdf2, *na_vcm, na_tot;
2927 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2932 * First calc 3xnr-atoms for each group
2933 * then subtract half a degree of freedom for each constraint
2935 * Only atoms and nuclei contribute to the degrees of freedom...
2940 const SimulationGroups& groups = mtop->groups;
2941 natoms = mtop->natoms;
2943 /* Allocate one more for a possible rest group */
2944 /* We need to sum degrees of freedom into doubles,
2945 * since floats give too low nrdf's above 3 million atoms.
2947 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
2948 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2949 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2950 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2951 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2953 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2957 for (gmx::index i = 0;
2958 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; i++)
2961 clear_ivec(dof_vcm[i]);
2963 nrdf_vcm_sub[i] = 0;
2965 snew(nrdf2, natoms);
2966 for (const AtomProxy atomP : AtomRange(*mtop))
2968 const t_atom& local = atomP.atom();
2969 int i = atomP.globalAtomNumber();
2971 if (local.ptype == eptAtom || local.ptype == eptNucleus)
2973 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
2974 for (int d = 0; d < DIM; d++)
2976 if (opts->nFreeze[g][d] == 0)
2978 /* Add one DOF for particle i (counted as 2*1) */
2980 /* VCM group i has dim d as a DOF */
2981 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
2985 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
2987 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
2993 for (const gmx_molblock_t& molb : mtop->molblock)
2995 const gmx_moltype_t& molt = mtop->moltype[molb.type];
2996 const t_atom* atom = molt.atoms.atom;
2997 for (int mol = 0; mol < molb.nmol; mol++)
2999 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
3001 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
3002 for (int i = 0; i < molt.ilist[ftype].size();)
3004 /* Subtract degrees of freedom for the constraints,
3005 * if the particles still have degrees of freedom left.
3006 * If one of the particles is a vsite or a shell, then all
3007 * constraint motion will go there, but since they do not
3008 * contribute to the constraints the degrees of freedom do not
3011 int ai = as + ia[i + 1];
3012 int aj = as + ia[i + 2];
3013 if (((atom[ia[i + 1]].ptype == eptNucleus) || (atom[ia[i + 1]].ptype == eptAtom))
3014 && ((atom[ia[i + 2]].ptype == eptNucleus) || (atom[ia[i + 2]].ptype == eptAtom)))
3032 imin = std::min(imin, nrdf2[ai]);
3033 jmin = std::min(jmin, nrdf2[aj]);
3036 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3038 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3040 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3042 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3045 i += interaction_function[ftype].nratoms + 1;
3048 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3049 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3051 /* Subtract 1 dof from every atom in the SETTLE */
3052 for (int j = 0; j < 3; j++)
3054 int ai = as + ia[i + 1 + j];
3055 imin = std::min(2, nrdf2[ai]);
3057 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3059 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3064 as += molt.atoms.nr;
3070 /* Correct nrdf for the COM constraints.
3071 * We correct using the TC and VCM group of the first atom
3072 * in the reference and pull group. If atoms in one pull group
3073 * belong to different TC or VCM groups it is anyhow difficult
3074 * to determine the optimal nrdf assignment.
3076 pull = ir->pull.get();
3078 for (int i = 0; i < pull->ncoord; i++)
3080 if (pull->coord[i].eType != epullCONSTRAINT)
3087 for (int j = 0; j < 2; j++)
3089 const t_pull_group* pgrp;
3091 pgrp = &pull->group[pull->coord[i].group[j]];
3093 if (!pgrp->ind.empty())
3095 /* Subtract 1/2 dof from each group */
3096 int ai = pgrp->ind[0];
3097 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3099 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3101 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3104 "Center of mass pulling constraints caused the number of degrees "
3105 "of freedom for temperature coupling group %s to be negative",
3106 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3107 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3112 /* We need to subtract the whole DOF from group j=1 */
3119 if (ir->nstcomm != 0)
3121 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3122 "Expect at least one group when removing COM motion");
3124 /* We remove COM motion up to dim ndof_com() */
3125 const int ndim_rm_vcm = ndof_com(ir);
3127 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3128 * the number of degrees of freedom in each vcm group when COM
3129 * translation is removed and 6 when rotation is removed as well.
3130 * Note that we do not and should not include the rest group here.
3132 for (gmx::index j = 0;
3133 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]); j++)
3135 switch (ir->comm_mode)
3138 case ecmLINEAR_ACCELERATION_CORRECTION:
3139 nrdf_vcm_sub[j] = 0;
3140 for (int d = 0; d < ndim_rm_vcm; d++)
3148 case ecmANGULAR: nrdf_vcm_sub[j] = 6; break;
3149 default: gmx_incons("Checking comm_mode");
3153 for (gmx::index i = 0;
3154 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
3156 /* Count the number of atoms of TC group i for every VCM group */
3157 for (gmx::index j = 0;
3158 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3163 for (int ai = 0; ai < natoms; ai++)
3165 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3167 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3171 /* Correct for VCM removal according to the fraction of each VCM
3172 * group present in this TC group.
3174 nrdf_uc = nrdf_tc[i];
3176 for (gmx::index j = 0;
3177 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3179 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3181 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3182 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3187 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3189 opts->nrdf[i] = nrdf_tc[i];
3190 if (opts->nrdf[i] < 0)
3194 fprintf(stderr, "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3195 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
3203 sfree(nrdf_vcm_sub);
3206 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3208 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3209 * But since this is much larger than STRLEN, such a line can not be parsed.
3210 * The real maximum is the number of names that fit in a string: STRLEN/2.
3212 #define EGP_MAX (STRLEN / 2)
3216 auto names = gmx::splitString(val);
3217 if (names.size() % 2 != 0)
3219 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3221 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3223 for (size_t i = 0; i < names.size() / 2; i++)
3225 // TODO this needs to be replaced by a solution using std::find_if
3229 names[2 * i].c_str(),
3230 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3236 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3241 names[2 * i + 1].c_str(),
3242 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3248 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3250 if ((j < nr) && (k < nr))
3252 ir->opts.egp_flags[nr * j + k] |= flag;
3253 ir->opts.egp_flags[nr * k + j] |= flag;
3262 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3264 int ig = -1, i = 0, gind;
3268 /* Just a quick check here, more thorough checks are in mdrun */
3269 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3271 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3274 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3275 for (ig = 0; ig < swap->ngrp; ig++)
3277 swapg = &swap->grp[ig];
3278 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3279 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3283 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3284 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap", swap->grp[ig].molname, swapg->nat);
3285 snew(swapg->ind, swapg->nat);
3286 for (i = 0; i < swapg->nat; i++)
3288 swapg->ind[i] = grps->a[grps->index[gind] + i];
3293 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3299 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3304 ig = search_string(IMDgname, grps->nr, gnames);
3305 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3307 if (IMDgroup->nat > 0)
3310 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3312 IMDgname, IMDgroup->nat);
3313 snew(IMDgroup->ind, IMDgroup->nat);
3314 for (i = 0; i < IMDgroup->nat; i++)
3316 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3321 /* Checks whether atoms are both part of a COM removal group and frozen.
3322 * If a fully frozen atom is part of a COM removal group, it is removed
3323 * from the COM removal group. A note is issued if such atoms are present.
3324 * A warning is issued for atom with one or two dimensions frozen that
3325 * are part of a COM removal group (mdrun would need to compute COM mass
3326 * per dimension to handle this correctly).
3327 * Also issues a warning when non-frozen atoms are not part of a COM
3328 * removal group while COM removal is active.
3330 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3332 const t_grpopts& opts,
3335 const int vcmRestGroup =
3336 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3338 int numFullyFrozenVcmAtoms = 0;
3339 int numPartiallyFrozenVcmAtoms = 0;
3340 int numNonVcmAtoms = 0;
3341 for (int a = 0; a < numAtoms; a++)
3343 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3344 int numFrozenDims = 0;
3345 for (int d = 0; d < DIM; d++)
3347 numFrozenDims += opts.nFreeze[freezeGroup][d];
3350 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3351 if (vcmGroup < vcmRestGroup)
3353 if (numFrozenDims == DIM)
3355 /* Do not remove COM motion for this fully frozen atom */
3356 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3358 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3361 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3362 numFullyFrozenVcmAtoms++;
3364 else if (numFrozenDims > 0)
3366 numPartiallyFrozenVcmAtoms++;
3369 else if (numFrozenDims < DIM)
3375 if (numFullyFrozenVcmAtoms > 0)
3377 std::string warningText = gmx::formatString(
3378 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3379 "removing these atoms from the COMM removal group(s)",
3380 numFullyFrozenVcmAtoms);
3381 warning_note(wi, warningText.c_str());
3383 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3385 std::string warningText = gmx::formatString(
3386 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3387 "removal group(s), due to limitations in the code these still contribute to the "
3388 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3390 numPartiallyFrozenVcmAtoms, DIM);
3391 warning(wi, warningText.c_str());
3393 if (numNonVcmAtoms > 0)
3395 std::string warningText = gmx::formatString(
3396 "%d atoms are not part of any center of mass motion removal group.\n"
3397 "This may lead to artifacts.\n"
3398 "In most cases one should use one group for the whole system.",
3400 warning(wi, warningText.c_str());
3404 void do_index(const char* mdparin,
3408 const gmx::MdModulesNotifier& notifier,
3412 t_blocka* defaultIndexGroups;
3420 int i, j, k, restnm;
3421 bool bExcl, bTable, bAnneal;
3422 char warn_buf[STRLEN];
3426 fprintf(stderr, "processing index file...\n");
3430 snew(defaultIndexGroups, 1);
3431 snew(defaultIndexGroups->index, 1);
3433 atoms_all = gmx_mtop_global_atoms(mtop);
3434 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3435 done_atom(&atoms_all);
3439 defaultIndexGroups = init_index(ndx, &gnames);
3442 SimulationGroups* groups = &mtop->groups;
3443 natoms = mtop->natoms;
3444 symtab = &mtop->symtab;
3446 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3448 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3450 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3451 restnm = groups->groupNames.size() - 1;
3452 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3453 srenew(gnames, defaultIndexGroups->nr + 1);
3454 gnames[restnm] = *(groups->groupNames.back());
3456 set_warning_line(wi, mdparin, -1);
3458 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3459 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3460 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3461 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3462 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3465 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3467 temperatureCouplingGroupNames.size(), temperatureCouplingReferenceValues.size(),
3468 temperatureCouplingTauValues.size());
3471 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3472 do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
3473 SimulationAtomGroupType::TemperatureCoupling, restnm,
3474 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3475 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3477 snew(ir->opts.nrdf, nr);
3478 snew(ir->opts.tau_t, nr);
3479 snew(ir->opts.ref_t, nr);
3480 if (ir->eI == eiBD && ir->bd_fric == 0)
3482 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3485 if (useReferenceTemperature)
3487 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3489 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3493 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3494 for (i = 0; (i < nr); i++)
3496 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3498 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3499 warning_error(wi, warn_buf);
3502 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3506 "tau-t = -1 is the value to signal that a group should not have "
3507 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3510 if (ir->opts.tau_t[i] >= 0)
3512 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3515 if (ir->etc != etcNO && ir->nsttcouple == -1)
3517 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3522 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3525 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3526 "md-vv; use either vrescale temperature with berendsen pressure or "
3527 "Nose-Hoover temperature with MTTK pressure");
3529 if (ir->epc == epcMTTK)
3531 if (ir->etc != etcNOSEHOOVER)
3534 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3539 if (ir->nstpcouple != ir->nsttcouple)
3541 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3542 ir->nstpcouple = ir->nsttcouple = mincouple;
3544 "for current Trotter decomposition methods with vv, nsttcouple and "
3545 "nstpcouple must be equal. Both have been reset to "
3546 "min(nsttcouple,nstpcouple) = %d",
3548 warning_note(wi, warn_buf);
3553 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3554 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3556 if (ETC_ANDERSEN(ir->etc))
3558 if (ir->nsttcouple != 1)
3562 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3563 "need for larger nsttcouple > 1, since no global parameters are computed. "
3564 "nsttcouple has been reset to 1");
3565 warning_note(wi, warn_buf);
3568 nstcmin = tcouple_min_integration_steps(ir->etc);
3571 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3574 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3575 "least %d times larger than nsttcouple*dt (%g)",
3576 ETCOUPLTYPE(ir->etc), tau_min, nstcmin, ir->nsttcouple * ir->delta_t);
3577 warning(wi, warn_buf);
3580 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3581 for (i = 0; (i < nr); i++)
3583 if (ir->opts.ref_t[i] < 0)
3585 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3588 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3589 if we are in this conditional) if mc_temp is negative */
3590 if (ir->expandedvals->mc_temp < 0)
3592 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3596 /* Simulated annealing for each group. There are nr groups */
3597 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3598 if (simulatedAnnealingGroupNames.size() == 1
3599 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3601 simulatedAnnealingGroupNames.resize(0);
3603 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3605 gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
3606 simulatedAnnealingGroupNames.size(), nr);
3610 snew(ir->opts.annealing, nr);
3611 snew(ir->opts.anneal_npoints, nr);
3612 snew(ir->opts.anneal_time, nr);
3613 snew(ir->opts.anneal_temp, nr);
3614 for (i = 0; i < nr; i++)
3616 ir->opts.annealing[i] = eannNO;
3617 ir->opts.anneal_npoints[i] = 0;
3618 ir->opts.anneal_time[i] = nullptr;
3619 ir->opts.anneal_temp[i] = nullptr;
3621 if (!simulatedAnnealingGroupNames.empty())
3624 for (i = 0; i < nr; i++)
3626 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3628 ir->opts.annealing[i] = eannNO;
3630 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3632 ir->opts.annealing[i] = eannSINGLE;
3635 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3637 ir->opts.annealing[i] = eannPERIODIC;
3643 /* Read the other fields too */
3644 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3645 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3647 gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
3648 simulatedAnnealingPoints.size(), simulatedAnnealingGroupNames.size());
3650 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3651 size_t numSimulatedAnnealingFields = 0;
3652 for (i = 0; i < nr; i++)
3654 if (ir->opts.anneal_npoints[i] == 1)
3658 "Please specify at least a start and an end point for annealing\n");
3660 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3661 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3662 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3665 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3667 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3669 gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
3670 simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
3672 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3673 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3675 gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
3676 simulatedAnnealingTemperatures.size(), numSimulatedAnnealingFields);
3679 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3680 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3681 convertReals(wi, simulatedAnnealingTimes, "anneal-time",
3682 allSimulatedAnnealingTimes.data());
3683 convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp",
3684 allSimulatedAnnealingTemperatures.data());
3685 for (i = 0, k = 0; i < nr; i++)
3687 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3689 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3690 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3693 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3695 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3701 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3704 "Annealing timepoints out of order: t=%f comes after "
3706 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j - 1]);
3709 if (ir->opts.anneal_temp[i][j] < 0)
3711 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n",
3712 ir->opts.anneal_temp[i][j]);
3717 /* Print out some summary information, to make sure we got it right */
3718 for (i = 0; i < nr; i++)
3720 if (ir->opts.annealing[i] != eannNO)
3722 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3723 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3724 *(groups->groupNames[j]), eann_names[ir->opts.annealing[i]],
3725 ir->opts.anneal_npoints[i]);
3726 fprintf(stderr, "Time (ps) Temperature (K)\n");
3727 /* All terms except the last one */
3728 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3730 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3731 ir->opts.anneal_temp[i][j]);
3734 /* Finally the last one */
3735 j = ir->opts.anneal_npoints[i] - 1;
3736 if (ir->opts.annealing[i] == eannSINGLE)
3738 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j],
3739 ir->opts.anneal_temp[i][j]);
3743 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3744 ir->opts.anneal_temp[i][j]);
3745 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3748 "There is a temperature jump when your annealing "
3760 process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3762 checkPullCoords(ir->pull->group, ir->pull->coord);
3767 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3770 if (ir->eSwapCoords != eswapNO)
3772 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3775 /* Make indices for IMD session */
3778 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3781 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3782 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3783 notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
3785 auto accelerations = gmx::splitString(inputrecStrings->acc);
3786 auto accelerationGroupNames = gmx::splitString(inputrecStrings->accgrps);
3787 if (accelerationGroupNames.size() * DIM != accelerations.size())
3789 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3790 accelerationGroupNames.size(), accelerations.size());
3792 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3793 SimulationAtomGroupType::Acceleration, restnm, egrptpALL_GENREST, bVerbose, wi);
3794 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3795 snew(ir->opts.acc, nr);
3796 ir->opts.ngacc = nr;
3798 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3800 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3801 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3802 if (freezeDims.size() != DIM * freezeGroupNames.size())
3804 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3805 freezeGroupNames.size(), freezeDims.size());
3807 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3808 SimulationAtomGroupType::Freeze, restnm, egrptpALL_GENREST, bVerbose, wi);
3809 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3810 ir->opts.ngfrz = nr;
3811 snew(ir->opts.nFreeze, nr);
3812 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3814 for (j = 0; (j < DIM); j++, k++)
3816 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3817 if (!ir->opts.nFreeze[i][j])
3819 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3822 "Please use Y(ES) or N(O) for freezedim only "
3824 freezeDims[k].c_str());
3825 warning(wi, warn_buf);
3830 for (; (i < nr); i++)
3832 for (j = 0; (j < DIM); j++)
3834 ir->opts.nFreeze[i][j] = 0;
3838 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3839 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3840 SimulationAtomGroupType::EnergyOutput, restnm, egrptpALL_GENREST, bVerbose, wi);
3841 add_wall_energrps(groups, ir->nwall, symtab);
3842 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3843 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3844 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3845 SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
3846 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3848 if (ir->comm_mode != ecmNO)
3850 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3853 /* Now we have filled the freeze struct, so we can calculate NRDF */
3854 calc_nrdf(mtop, ir, gnames);
3856 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3857 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3858 SimulationAtomGroupType::User1, restnm, egrptpALL_GENREST, bVerbose, wi);
3859 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3860 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3861 SimulationAtomGroupType::User2, restnm, egrptpALL_GENREST, bVerbose, wi);
3862 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3863 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3864 SimulationAtomGroupType::CompressedPositionOutput, restnm, egrptpONE, bVerbose, wi);
3865 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3866 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3867 SimulationAtomGroupType::OrientationRestraintsFit, restnm, egrptpALL_GENREST,
3870 /* MiMiC QMMM input processing */
3871 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
3872 if (qmGroupNames.size() > 1)
3874 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3876 /* group rest, if any, is always MM! */
3877 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3878 SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
3879 ir->opts.ngQM = qmGroupNames.size();
3881 /* end of MiMiC QMMM input */
3885 for (auto group : gmx::keysOf(groups->groups))
3887 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3888 for (const auto& entry : groups->groups[group])
3890 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3892 fprintf(stderr, "\n");
3896 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3897 snew(ir->opts.egp_flags, nr * nr);
3899 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
3900 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3902 warning_error(wi, "Energy group exclusions are currently not supported");
3904 if (bExcl && EEL_FULL(ir->coulombtype))
3906 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3909 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
3910 if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
3911 && !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
3914 "Can only have energy group pair tables in combination with user tables for VdW "
3918 /* final check before going out of scope if simulated tempering variables
3919 * need to be set to default values.
3921 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3923 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
3924 warning(wi, gmx::formatString(
3925 "the value for nstexpanded was not specified for "
3926 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3927 "by default, but it is recommended to set it to an explicit value!",
3928 ir->expandedvals->nstexpanded));
3930 for (i = 0; (i < defaultIndexGroups->nr); i++)
3935 done_blocka(defaultIndexGroups);
3936 sfree(defaultIndexGroups);
3940 static void check_disre(const gmx_mtop_t* mtop)
3942 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3944 const gmx_ffparams_t& ffparams = mtop->ffparams;
3947 for (int i = 0; i < ffparams.numTypes(); i++)
3949 int ftype = ffparams.functype[i];
3950 if (ftype == F_DISRES)
3952 int label = ffparams.iparams[i].disres.label;
3953 if (label == old_label)
3955 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3964 "Found %d double distance restraint indices,\n"
3965 "probably the parameters for multiple pairs in one restraint "
3966 "are not identical\n",
3972 static bool absolute_reference(const t_inputrec* ir, const gmx_mtop_t* sys, const bool posres_only, ivec AbsRef)
3975 gmx_mtop_ilistloop_t iloop;
3977 const t_iparams* pr;
3984 for (d = 0; d < DIM; d++)
3986 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3988 /* Check for freeze groups */
3989 for (g = 0; g < ir->opts.ngfrz; g++)
3991 for (d = 0; d < DIM; d++)
3993 if (ir->opts.nFreeze[g][d] != 0)
4001 /* Check for position restraints */
4002 iloop = gmx_mtop_ilistloop_init(sys);
4003 while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
4005 if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
4007 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
4009 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
4010 for (d = 0; d < DIM; d++)
4012 if (pr->posres.fcA[d] != 0)
4018 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
4020 /* Check for flat-bottom posres */
4021 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
4022 if (pr->fbposres.k != 0)
4024 switch (pr->fbposres.geom)
4026 case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
4027 case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
4028 case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
4029 case efbposresCYLINDER:
4030 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4031 case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
4032 case efbposresX: /* d=XX */
4033 case efbposresY: /* d=YY */
4034 case efbposresZ: /* d=ZZ */
4035 d = pr->fbposres.geom - efbposresX;
4040 " Invalid geometry for flat-bottom position restraint.\n"
4041 "Expected nr between 1 and %d. Found %d\n",
4042 efbposresNR - 1, pr->fbposres.geom);
4049 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
4052 static void check_combination_rule_differences(const gmx_mtop_t* mtop,
4054 bool* bC6ParametersWorkWithGeometricRules,
4055 bool* bC6ParametersWorkWithLBRules,
4056 bool* bLBRulesPossible)
4058 int ntypes, tpi, tpj;
4061 double c6i, c6j, c12i, c12j;
4062 double c6, c6_geometric, c6_LB;
4063 double sigmai, sigmaj, epsi, epsj;
4064 bool bCanDoLBRules, bCanDoGeometricRules;
4067 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4068 * force-field floating point parameters.
4071 ptr = getenv("GMX_LJCOMB_TOL");
4075 double gmx_unused canary;
4077 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4080 "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4085 *bC6ParametersWorkWithLBRules = TRUE;
4086 *bC6ParametersWorkWithGeometricRules = TRUE;
4087 bCanDoLBRules = TRUE;
4088 ntypes = mtop->ffparams.atnr;
4089 snew(typecount, ntypes);
4090 gmx_mtop_count_atomtypes(mtop, state, typecount);
4091 *bLBRulesPossible = TRUE;
4092 for (tpi = 0; tpi < ntypes; ++tpi)
4094 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4095 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4096 for (tpj = tpi; tpj < ntypes; ++tpj)
4098 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4099 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4100 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4101 c6_geometric = std::sqrt(c6i * c6j);
4102 if (!gmx_numzero(c6_geometric))
4104 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4106 sigmai = gmx::sixthroot(c12i / c6i);
4107 sigmaj = gmx::sixthroot(c12j / c6j);
4108 epsi = c6i * c6i / (4.0 * c12i);
4109 epsj = c6j * c6j / (4.0 * c12j);
4110 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4114 *bLBRulesPossible = FALSE;
4115 c6_LB = c6_geometric;
4117 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4122 *bC6ParametersWorkWithLBRules = FALSE;
4125 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4127 if (!bCanDoGeometricRules)
4129 *bC6ParametersWorkWithGeometricRules = FALSE;
4136 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop, warninp_t wi)
4138 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4140 check_combination_rule_differences(mtop, 0, &bC6ParametersWorkWithGeometricRules,
4141 &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4142 if (ir->ljpme_combination_rule == eljpmeLB)
4144 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4147 "You are using arithmetic-geometric combination rules "
4148 "in LJ-PME, but your non-bonded C6 parameters do not "
4149 "follow these rules.");
4154 if (!bC6ParametersWorkWithGeometricRules)
4156 if (ir->eDispCorr != edispcNO)
4159 "You are using geometric combination rules in "
4160 "LJ-PME, but your non-bonded C6 parameters do "
4161 "not follow these rules. "
4162 "This will introduce very small errors in the forces and energies in "
4163 "your simulations. Dispersion correction will correct total energy "
4164 "and/or pressure for isotropic systems, but not forces or surface "
4170 "You are using geometric combination rules in "
4171 "LJ-PME, but your non-bonded C6 parameters do "
4172 "not follow these rules. "
4173 "This will introduce very small errors in the forces and energies in "
4174 "your simulations. If your system is homogeneous, consider using "
4175 "dispersion correction "
4176 "for the total energy and pressure.");
4182 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4184 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4185 gmx::assertMtsRequirements(*ir);
4187 char err_buf[STRLEN];
4192 gmx_mtop_atomloop_block_t aloopb;
4194 char warn_buf[STRLEN];
4196 set_warning_line(wi, mdparin, -1);
4198 if (absolute_reference(ir, sys, false, AbsRef))
4201 "Removing center of mass motion in the presence of position restraints might "
4202 "cause artifacts. When you are using position restraints to equilibrate a "
4203 "macro-molecule, the artifacts are usually negligible.");
4206 if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 && ir->nstlist > 1
4207 && ((EI_MD(ir->eI) || EI_SD(ir->eI)) && (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4209 /* Check if a too small Verlet buffer might potentially
4210 * cause more drift than the thermostat can couple off.
4212 /* Temperature error fraction for warning and suggestion */
4213 const real T_error_warn = 0.002;
4214 const real T_error_suggest = 0.001;
4215 /* For safety: 2 DOF per atom (typical with constraints) */
4216 const real nrdf_at = 2;
4217 real T, tau, max_T_error;
4222 for (i = 0; i < ir->opts.ngtc; i++)
4224 T = std::max(T, ir->opts.ref_t[i]);
4225 tau = std::max(tau, ir->opts.tau_t[i]);
4229 /* This is a worst case estimate of the temperature error,
4230 * assuming perfect buffer estimation and no cancelation
4231 * of errors. The factor 0.5 is because energy distributes
4232 * equally over Ekin and Epot.
4234 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * BOLTZ * T);
4235 if (max_T_error > T_error_warn)
4238 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4239 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4240 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4241 "%.0e or decrease tau_t.",
4242 ir->verletbuf_tol, T, tau, 100 * max_T_error, 100 * T_error_suggest,
4243 ir->verletbuf_tol * T_error_suggest / max_T_error);
4244 warning(wi, warn_buf);
4249 if (ETC_ANDERSEN(ir->etc))
4253 for (i = 0; i < ir->opts.ngtc; i++)
4256 "all tau_t must currently be equal using Andersen temperature control, "
4257 "violated for group %d",
4259 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4261 "all tau_t must be positive using Andersen temperature control, "
4263 i, ir->opts.tau_t[i]);
4264 CHECK(ir->opts.tau_t[i] < 0);
4267 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
4269 for (i = 0; i < ir->opts.ngtc; i++)
4271 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4273 "tau_t/delta_t for group %d for temperature control method %s must be a "
4274 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4275 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4277 i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
4278 CHECK(nsteps % ir->nstcomm != 0);
4283 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD && ir->comm_mode == ecmNO
4284 && !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
4287 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4288 "rounding errors can lead to build up of kinetic energy of the center of mass");
4291 if (ir->epc == epcPARRINELLORAHMAN && ir->etc == etcNOSEHOOVER)
4294 for (int g = 0; g < ir->opts.ngtc; g++)
4296 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4298 if (ir->tau_p < 1.9 * tau_t_max)
4300 std::string message = gmx::formatString(
4301 "With %s T-coupling and %s p-coupling, "
4302 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4303 etcoupl_names[ir->etc], epcoupl_names[ir->epc], "tau-p", ir->tau_p, "tau-t",
4305 warning(wi, message.c_str());
4309 /* Check for pressure coupling with absolute position restraints */
4310 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4312 absolute_reference(ir, sys, TRUE, AbsRef);
4314 for (m = 0; m < DIM; m++)
4316 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4319 "You are using pressure coupling with absolute position restraints, "
4320 "this will give artifacts. Use the refcoord_scaling option.");
4328 aloopb = gmx_mtop_atomloop_block_init(sys);
4330 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4332 if (atom->q != 0 || atom->qB != 0)
4340 if (EEL_FULL(ir->coulombtype))
4343 "You are using full electrostatics treatment %s for a system without charges.\n"
4344 "This costs a lot of performance for just processing zeros, consider using %s "
4346 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4347 warning(wi, err_buf);
4352 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4355 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4356 "You might want to consider using %s electrostatics.\n",
4358 warning_note(wi, err_buf);
4362 /* Check if combination rules used in LJ-PME are the same as in the force field */
4363 if (EVDW_PME(ir->vdwtype))
4365 check_combination_rules(ir, sys, wi);
4368 /* Generalized reaction field */
4369 if (ir->coulombtype == eelGRF_NOTUSED)
4372 "Generalized reaction-field electrostatics is no longer supported. "
4373 "You can use normal reaction-field instead and compute the reaction-field "
4374 "constant by hand.");
4378 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4380 for (m = 0; (m < DIM); m++)
4382 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4391 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4392 for (const AtomProxy atomP : AtomRange(*sys))
4394 const t_atom& local = atomP.atom();
4395 int i = atomP.globalAtomNumber();
4396 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4399 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4401 for (m = 0; (m < DIM); m++)
4403 acc[m] += ir->opts.acc[i][m] * mgrp[i];
4407 for (m = 0; (m < DIM); m++)
4409 if (fabs(acc[m]) > 1e-6)
4411 const char* dim[DIM] = { "X", "Y", "Z" };
4412 fprintf(stderr, "Net Acceleration in %s direction, will %s be corrected\n", dim[m],
4413 ir->nstcomm != 0 ? "" : "not");
4414 if (ir->nstcomm != 0 && m < ndof_com(ir))
4418 (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4420 ir->opts.acc[i][m] -= acc[m];
4428 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0
4429 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4431 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4439 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4441 if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
4443 absolute_reference(ir, sys, FALSE, AbsRef);
4444 for (m = 0; m < DIM; m++)
4446 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4449 "You are using an absolute reference for pulling, but the rest of "
4450 "the system does not have an absolute reference. This will lead to "
4459 for (i = 0; i < 3; i++)
4461 for (m = 0; m <= i; m++)
4463 if ((ir->epc != epcNO && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4465 for (c = 0; c < ir->pull->ncoord; c++)
4467 if (ir->pull->coord[c].eGeom == epullgDIRPBC && ir->pull->coord[c].vec[m] != 0)
4470 "Can not have dynamic box while using pull geometry '%s' "
4472 EPULLGEOM(ir->pull->coord[c].eGeom), 'x' + m);
4483 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4485 char warn_buf[STRLEN];
4488 ptr = check_box(ir->pbcType, box);
4491 warning_error(wi, ptr);
4494 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4496 if (ir->shake_tol <= 0.0)
4498 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4499 warning_error(wi, warn_buf);
4503 if ((ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4505 /* If we have Lincs constraints: */
4506 if (ir->eI == eiMD && ir->etc == etcNO && ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4509 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4510 warning_note(wi, warn_buf);
4513 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4516 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4518 warning_note(wi, warn_buf);
4520 if (ir->epc == epcMTTK)
4522 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4526 if (bHasAnyConstraints && ir->epc == epcMTTK)
4528 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4531 if (ir->LincsWarnAngle > 90.0)
4533 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4534 warning(wi, warn_buf);
4535 ir->LincsWarnAngle = 90.0;
4538 if (ir->pbcType != PbcType::No)
4540 if (ir->nstlist == 0)
4543 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4544 "atoms might cause the simulation to crash.");
4546 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4549 "ERROR: The cut-off length is longer than half the shortest box vector or "
4550 "longer than the smallest box diagonal element. Increase the box size or "
4551 "decrease rlist.\n");
4552 warning_error(wi, warn_buf);