2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, The GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
50 #include "gromacs/applied_forces/awh/read_params.h"
51 #include "gromacs/fileio/readinp.h"
52 #include "gromacs/fileio/warninp.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/gmxpreprocess/toputil.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/calc_verletbuf.h"
59 #include "gromacs/mdrun/mdmodules.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/mdtypes/multipletimestepping.h"
63 #include "gromacs/mdtypes/pull_params.h"
64 #include "gromacs/options/options.h"
65 #include "gromacs/options/treesupport.h"
66 #include "gromacs/pbcutil/pbc.h"
67 #include "gromacs/selection/indexutil.h"
68 #include "gromacs/topology/block.h"
69 #include "gromacs/topology/ifunc.h"
70 #include "gromacs/topology/index.h"
71 #include "gromacs/topology/mtop_util.h"
72 #include "gromacs/topology/symtab.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/exceptions.h"
76 #include "gromacs/utility/fatalerror.h"
77 #include "gromacs/utility/filestream.h"
78 #include "gromacs/utility/gmxassert.h"
79 #include "gromacs/utility/ikeyvaluetreeerror.h"
80 #include "gromacs/utility/keyvaluetree.h"
81 #include "gromacs/utility/keyvaluetreebuilder.h"
82 #include "gromacs/utility/keyvaluetreemdpwriter.h"
83 #include "gromacs/utility/keyvaluetreetransform.h"
84 #include "gromacs/utility/mdmodulenotification.h"
85 #include "gromacs/utility/smalloc.h"
86 #include "gromacs/utility/strconvert.h"
87 #include "gromacs/utility/stringcompare.h"
88 #include "gromacs/utility/stringutil.h"
89 #include "gromacs/utility/textwriter.h"
94 /* Resource parameters
95 * Do not change any of these until you read the instruction
96 * in readinp.h. Some cpp's do not take spaces after the backslash
97 * (like the c-shell), which will give you a very weird compiler
101 struct gmx_inputrec_strings
103 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], acc[STRLEN], accgrps[STRLEN], freeze[STRLEN],
104 frdim[STRLEN], energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN],
105 x_compressed_groups[STRLEN], couple_moltype[STRLEN], orirefitgrp[STRLEN],
106 egptable[STRLEN], egpexcl[STRLEN], wall_atomtype[STRLEN], wall_density[STRLEN],
107 deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
108 char fep_lambda[efptNR][STRLEN];
109 char lambda_weights[STRLEN];
110 std::vector<std::string> pullGroupNames;
111 std::vector<std::string> rotateGroupNames;
112 char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
115 static gmx_inputrec_strings* inputrecStrings = nullptr;
117 void init_inputrec_strings()
122 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
123 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
125 inputrecStrings = new gmx_inputrec_strings();
128 void done_inputrec_strings()
130 delete inputrecStrings;
131 inputrecStrings = nullptr;
137 egrptpALL, /* All particles have to be a member of a group. */
138 egrptpALL_GENREST, /* A rest group with name is generated for particles *
139 * that are not part of any group. */
140 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
141 * for the rest group. */
142 egrptpONE /* Merge all selected groups into one group, *
143 * make a rest group for the remaining particles. */
146 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
147 "h-angles", "all-angles", nullptr };
149 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
151 static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lambdas)
156 for (i = 0; i < ntemps; i++)
158 /* simple linear scaling -- allows more control */
159 if (simtemp->eSimTempScale == esimtempLINEAR)
161 simtemp->temperatures[i] =
163 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
165 else if (simtemp->eSimTempScale
166 == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
168 simtemp->temperatures[i] = simtemp->simtemp_low
169 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
170 static_cast<real>((1.0 * i) / (ntemps - 1)));
172 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
174 simtemp->temperatures[i] = simtemp->simtemp_low
175 + (simtemp->simtemp_high - simtemp->simtemp_low)
176 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
181 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
182 gmx_fatal(FARGS, "%s", errorstr);
188 static void _low_check(bool b, const char* s, warninp_t wi)
192 warning_error(wi, s);
196 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
200 if (*p > 0 && *p % nst != 0)
202 /* Round up to the next multiple of nst */
203 *p = ((*p) / nst + 1) * nst;
204 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
209 static int lcd(int n1, int n2)
214 for (i = 2; (i <= n1 && i <= n2); i++)
216 if (n1 % i == 0 && n2 % i == 0)
225 //! Convert legacy mdp entries to modern ones.
226 static void process_interaction_modifier(int* eintmod)
228 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
230 *eintmod = eintmodPOTSHIFT;
234 static void checkMtsRequirement(const t_inputrec& ir, const char* param, const int nstValue, warninp_t wi)
236 GMX_RELEASE_ASSERT(ir.mtsLevels.size() >= 2, "Need at least two levels for MTS");
237 const int mtsFactor = ir.mtsLevels.back().stepFactor;
238 if (nstValue % mtsFactor != 0)
240 auto message = gmx::formatString(
241 "With MTS, %s = %d should be a multiple of mts-factor = %d", param, nstValue, mtsFactor);
242 warning_error(wi, message.c_str());
246 static void setupMtsLevels(gmx::ArrayRef<gmx::MtsLevel> mtsLevels,
247 const t_inputrec& ir,
248 const t_gromppopts& opts,
251 if (!(ir.eI == eiMD || ir.eI == eiSD1))
253 auto message = gmx::formatString(
254 "Multiple time stepping is only supported with integrators %s and %s",
255 ei_names[eiMD], ei_names[eiSD1]);
256 warning_error(wi, message.c_str());
258 if (opts.numMtsLevels != 2)
260 warning_error(wi, "Only mts-levels = 2 is supported");
264 const std::vector<std::string> inputForceGroups = gmx::splitString(opts.mtsLevel2Forces);
265 auto& forceGroups = mtsLevels[1].forceGroups;
266 for (const auto& inputForceGroup : inputForceGroups)
270 for (const auto& forceGroupName : gmx::mtsForceGroupNames)
272 if (gmx::equalCaseInsensitive(inputForceGroup, forceGroupName))
274 forceGroups.set(nameIndex);
282 gmx::formatString("Unknown MTS force group '%s'", inputForceGroup.c_str());
283 warning_error(wi, message.c_str());
287 if (mtsLevels[1].stepFactor <= 1)
289 gmx_fatal(FARGS, "mts-factor should be larger than 1");
292 // Make the level 0 use the complement of the force groups of group 1
293 mtsLevels[0].forceGroups = ~mtsLevels[1].forceGroups;
294 mtsLevels[0].stepFactor = 1;
296 if ((EEL_FULL(ir.coulombtype) || EVDW_PME(ir.vdwtype))
297 && !mtsLevels[1].forceGroups[static_cast<int>(gmx::MtsForceGroups::LongrangeNonbonded)])
300 "With long-range electrostatics and/or LJ treatment, the long-range part "
301 "has to be part of the mts-level2-forces");
304 if (ir.nstcalcenergy > 0)
306 checkMtsRequirement(ir, "nstcalcenergy", ir.nstcalcenergy, wi);
308 checkMtsRequirement(ir, "nstenergy", ir.nstenergy, wi);
309 checkMtsRequirement(ir, "nstlog", ir.nstlog, wi);
310 if (ir.efep != efepNO)
312 checkMtsRequirement(ir, "nstdhdl", ir.fepvals->nstdhdl, wi);
317 const int pullMtsLevel = gmx::forceGroupMtsLevel(ir.mtsLevels, gmx::MtsForceGroups::Pull);
318 if (ir.pull->nstxout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
320 warning_error(wi, "pull-nstxout should be a multiple of mts-factor");
322 if (ir.pull->nstfout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
324 warning_error(wi, "pull-nstfout should be a multiple of mts-factor");
330 void check_ir(const char* mdparin,
331 const gmx::MdModulesNotifier& mdModulesNotifier,
335 /* Check internal consistency.
336 * NOTE: index groups are not set here yet, don't check things
337 * like temperature coupling group options here, but in triple_check
340 /* Strange macro: first one fills the err_buf, and then one can check
341 * the condition, which will print the message and increase the error
344 #define CHECK(b) _low_check(b, err_buf, wi)
345 char err_buf[256], warn_buf[STRLEN];
348 t_lambda* fep = ir->fepvals;
349 t_expanded* expand = ir->expandedvals;
351 set_warning_line(wi, mdparin, -1);
353 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
355 sprintf(warn_buf, "%s electrostatics is no longer supported", eel_names[eelRF_NEC_UNSUPPORTED]);
356 warning_error(wi, warn_buf);
359 /* BASIC CUT-OFF STUFF */
360 if (ir->rcoulomb < 0)
362 warning_error(wi, "rcoulomb should be >= 0");
366 warning_error(wi, "rvdw should be >= 0");
368 if (ir->rlist < 0 && !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
370 warning_error(wi, "rlist should be >= 0");
373 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
374 "neighbour-list update scheme for efficient buffering for improved energy "
375 "conservation, please use the Verlet cut-off scheme instead.)");
376 CHECK(ir->nstlist < 0);
378 process_interaction_modifier(&ir->coulomb_modifier);
379 process_interaction_modifier(&ir->vdw_modifier);
381 if (ir->cutoff_scheme == ecutsGROUP)
384 "The group cutoff scheme has been removed since GROMACS 2020. "
385 "Please use the Verlet cutoff scheme.");
387 if (ir->cutoff_scheme == ecutsVERLET)
391 /* Normal Verlet type neighbor-list, currently only limited feature support */
392 if (inputrec2nboundeddim(ir) < 3)
394 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
397 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
398 if (ir->rcoulomb != ir->rvdw)
400 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
401 // for PME load balancing, we can support this exception.
402 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == evdwCUT
403 && ir->rcoulomb > ir->rvdw);
404 if (!bUsesPmeTwinRangeKernel)
407 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
408 "rcoulomb>rvdw with PME electrostatics)");
412 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
414 if (ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT)
416 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
419 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
421 evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
422 warning_note(wi, warn_buf);
424 ir->vdwtype = evdwCUT;
428 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
429 evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
430 warning_error(wi, warn_buf);
434 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
437 "With Verlet lists only cut-off and PME LJ interactions are supported");
439 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) || EEL_PME(ir->coulombtype)
440 || ir->coulombtype == eelEWALD))
443 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
444 "electrostatics are supported");
446 if (!(ir->coulomb_modifier == eintmodNONE || ir->coulomb_modifier == eintmodPOTSHIFT))
448 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
449 warning_error(wi, warn_buf);
452 if (EEL_USER(ir->coulombtype))
454 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme",
455 eel_names[ir->coulombtype]);
456 warning_error(wi, warn_buf);
459 if (ir->nstlist <= 0)
461 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
464 if (ir->nstlist < 10)
467 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
468 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
472 rc_max = std::max(ir->rvdw, ir->rcoulomb);
476 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
477 ir->verletbuf_tol = 0;
480 else if (ir->verletbuf_tol <= 0)
482 if (ir->verletbuf_tol == 0)
484 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
487 if (ir->rlist < rc_max)
490 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
493 if (ir->rlist == rc_max && ir->nstlist > 1)
497 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
498 "buffer. The cluster pair list does have a buffering effect, but choosing "
499 "a larger rlist might be necessary for good energy conservation.");
504 if (ir->rlist > rc_max)
507 "You have set rlist larger than the interaction cut-off, but you also "
508 "have verlet-buffer-tolerance > 0. Will set rlist using "
509 "verlet-buffer-tolerance.");
512 if (ir->nstlist == 1)
514 /* No buffer required */
519 if (EI_DYNAMICS(ir->eI))
521 if (inputrec2nboundeddim(ir) < 3)
524 "The box volume is required for calculating rlist from the "
525 "energy drift with verlet-buffer-tolerance > 0. You are "
526 "using at least one unbounded dimension, so no volume can be "
527 "computed. Either use a finite box, or set rlist yourself "
528 "together with verlet-buffer-tolerance = -1.");
530 /* Set rlist temporarily so we can continue processing */
535 /* Set the buffer to 5% of the cut-off */
536 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
542 /* GENERAL INTEGRATOR STUFF */
545 if (ir->etc != etcNO)
547 if (EI_RANDOM(ir->eI))
550 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
551 "implicitly. See the documentation for more information on which "
552 "parameters affect temperature for %s.",
553 etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
558 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
560 etcoupl_names[ir->etc], ei_names[ir->eI]);
562 warning_note(wi, warn_buf);
566 if (ir->eI == eiVVAK)
569 "Integrator method %s is implemented primarily for validation purposes; for "
570 "molecular dynamics, you should probably be using %s or %s",
571 ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
572 warning_note(wi, warn_buf);
574 if (!EI_DYNAMICS(ir->eI))
576 if (ir->epc != epcNO)
579 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
580 epcoupl_names[ir->epc], ei_names[ir->eI]);
581 warning_note(wi, warn_buf);
585 if (EI_DYNAMICS(ir->eI))
587 if (ir->nstcalcenergy < 0)
589 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
590 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
592 /* nstcalcenergy larger than nstener does not make sense.
593 * We ideally want nstcalcenergy=nstener.
597 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
601 ir->nstcalcenergy = ir->nstenergy;
605 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
606 || (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
607 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
610 const char* nsten = "nstenergy";
611 const char* nstdh = "nstdhdl";
612 const char* min_name = nsten;
613 int min_nst = ir->nstenergy;
615 /* find the smallest of ( nstenergy, nstdhdl ) */
616 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
617 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
619 min_nst = ir->fepvals->nstdhdl;
622 /* If the user sets nstenergy small, we should respect that */
623 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy,
625 warning_note(wi, warn_buf);
626 ir->nstcalcenergy = min_nst;
629 if (ir->epc != epcNO)
631 if (ir->nstpcouple < 0)
633 ir->nstpcouple = ir_optimal_nstpcouple(ir);
635 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
638 "With multiple time stepping, nstpcouple should be a mutiple of "
643 if (ir->nstcalcenergy > 0)
645 if (ir->efep != efepNO)
647 /* nstdhdl should be a multiple of nstcalcenergy */
648 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
652 /* nstexpanded should be a multiple of nstcalcenergy */
653 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded",
654 &ir->expandedvals->nstexpanded, wi);
656 /* for storing exact averages nstenergy should be
657 * a multiple of nstcalcenergy
659 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
662 // Inquire all MdModules, if their parameters match with the energy
663 // calculation frequency
664 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
665 mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
667 // Emit all errors from the energy calculation frequency checks
668 for (const std::string& energyFrequencyErrorMessage :
669 energyCalculationFrequencyErrors.errorMessages())
671 warning_error(wi, energyFrequencyErrorMessage);
675 if (ir->nsteps == 0 && !ir->bContinuation)
678 "For a correct single-point energy evaluation with nsteps = 0, use "
679 "continuation = yes to avoid constraining the input coordinates.");
683 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
686 "You are doing a continuation with SD or BD, make sure that ld_seed is "
687 "different from the previous run (using ld_seed=-1 will ensure this)");
693 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
694 CHECK(ir->pbcType != PbcType::Xyz);
695 sprintf(err_buf, "with TPI nstlist should be larger than zero");
696 CHECK(ir->nstlist <= 0);
697 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
698 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
702 if ((opts->nshake > 0) && (opts->bMorse))
704 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
705 warning(wi, warn_buf);
708 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
711 "You are doing a continuation with SD or BD, make sure that ld_seed is "
712 "different from the previous run (using ld_seed=-1 will ensure this)");
714 /* verify simulated tempering options */
718 bool bAllTempZero = TRUE;
719 for (i = 0; i < fep->n_lambda; i++)
721 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
722 efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
723 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
724 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
726 bAllTempZero = FALSE;
729 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
730 CHECK(bAllTempZero == TRUE);
732 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
733 CHECK(ir->eI != eiVV);
735 /* check compatability of the temperature coupling with simulated tempering */
737 if (ir->etc == etcNOSEHOOVER)
740 "Nose-Hoover based temperature control such as [%s] my not be "
741 "entirelyconsistent with simulated tempering",
742 etcoupl_names[ir->etc]);
743 warning_note(wi, warn_buf);
746 /* check that the temperatures make sense */
749 "Higher simulated tempering temperature (%g) must be >= than the simulated "
750 "tempering lower temperature (%g)",
751 ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
752 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
754 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero",
755 ir->simtempvals->simtemp_high);
756 CHECK(ir->simtempvals->simtemp_high <= 0);
758 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero",
759 ir->simtempvals->simtemp_low);
760 CHECK(ir->simtempvals->simtemp_low <= 0);
763 /* verify free energy options */
765 if (ir->efep != efepNO)
768 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
769 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
772 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
774 static_cast<int>(fep->sc_r_power));
775 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
778 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
781 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
783 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
785 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
787 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
788 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
790 sprintf(err_buf, "Free-energy not implemented for Ewald");
791 CHECK(ir->coulombtype == eelEWALD);
793 /* check validty of lambda inputs */
794 if (fep->n_lambda == 0)
796 /* Clear output in case of no states:*/
797 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.",
798 fep->init_fep_state);
799 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
803 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d",
804 fep->init_fep_state, fep->n_lambda - 1);
805 CHECK((fep->init_fep_state >= fep->n_lambda));
809 "Lambda state must be set, either with init-lambda-state or with init-lambda");
810 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
813 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
814 "init-lambda-state or with init-lambda, but not both",
815 fep->init_lambda, fep->init_fep_state);
816 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
819 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
823 for (i = 0; i < efptNR; i++)
825 if (fep->separate_dvdl[i])
830 if (n_lambda_terms > 1)
833 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
834 "use init-lambda to set lambda state (except for slow growth). Use "
835 "init-lambda-state instead.");
836 warning(wi, warn_buf);
839 if (n_lambda_terms < 2 && fep->n_lambda > 0)
842 "init-lambda is deprecated for setting lambda state (except for slow "
843 "growth). Use init-lambda-state instead.");
847 for (j = 0; j < efptNR; j++)
849 for (i = 0; i < fep->n_lambda; i++)
851 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
852 efpt_names[j], fep->all_lambda[j][i]);
853 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
857 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
859 for (i = 0; i < fep->n_lambda; i++)
862 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
863 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
864 "crashes, and is not supported.",
865 i, fep->all_lambda[efptVDW][i], fep->all_lambda[efptCOUL][i]);
866 CHECK((fep->sc_alpha > 0)
867 && (((fep->all_lambda[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL][i] < 1.0))
868 && ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0))));
872 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
874 real sigma, lambda, r_sc;
877 /* Maximum estimate for A and B charges equal with lambda power 1 */
879 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
880 1.0 / fep->sc_r_power);
882 "With PME there is a minor soft core effect present at the cut-off, "
883 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
884 "energy conservation, but usually other effects dominate. With a common sigma "
885 "value of %g nm the fraction of the particle-particle potential at the cut-off "
886 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
887 fep->sc_r_power, sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
888 warning_note(wi, warn_buf);
891 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
892 be treated differently, but that's the next step */
894 for (i = 0; i < efptNR; i++)
896 for (j = 0; j < fep->n_lambda; j++)
898 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
899 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
904 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
908 /* checking equilibration of weights inputs for validity */
911 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
913 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
914 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
917 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
919 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
920 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
923 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
924 expand->equil_steps, elmceq_names[elmceqSTEPS]);
925 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
928 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
929 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
930 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
933 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
934 expand->equil_ratio, elmceq_names[elmceqRATIO]);
935 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
938 "weight-equil-number-all-lambda (%d) must be a positive integer if "
939 "lmc-weights-equil=%s",
940 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
941 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
944 "weight-equil-number-samples (%d) must be a positive integer if "
945 "lmc-weights-equil=%s",
946 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
947 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
950 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
951 expand->equil_steps, elmceq_names[elmceqSTEPS]);
952 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
954 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
955 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
956 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
958 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
959 expand->equil_ratio, elmceq_names[elmceqRATIO]);
960 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
962 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
963 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
964 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
966 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
967 CHECK((expand->lmc_repeats <= 0));
968 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
969 CHECK((expand->minvarmin <= 0));
970 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
971 CHECK((expand->c_range < 0));
973 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
975 fep->init_fep_state, expand->lmc_forced_nstart);
976 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
977 && (expand->elmcmove != elmcmoveNO));
978 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
979 CHECK((expand->lmc_forced_nstart < 0));
980 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
981 fep->init_fep_state);
982 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
984 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
985 CHECK((expand->init_wl_delta < 0));
986 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
987 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
988 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
989 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
991 /* if there is no temperature control, we need to specify an MC temperature */
992 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO)
993 && (expand->mc_temp <= 0.0))
996 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
997 "to a positive number");
998 warning_error(wi, err_buf);
1000 if (expand->nstTij > 0)
1002 sprintf(err_buf, "nstlog must be non-zero");
1003 CHECK(ir->nstlog == 0);
1004 // Avoid modulus by zero in the case that already triggered an error exit.
1005 if (ir->nstlog != 0)
1008 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
1009 expand->nstTij, ir->nstlog);
1010 CHECK((expand->nstTij % ir->nstlog) != 0);
1016 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
1017 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
1020 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
1022 if (ir->pbcType == PbcType::No)
1024 if (ir->epc != epcNO)
1026 warning(wi, "Turning off pressure coupling for vacuum system");
1032 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
1033 c_pbcTypeNames[ir->pbcType].c_str());
1034 CHECK(ir->epc != epcNO);
1036 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
1037 CHECK(EEL_FULL(ir->coulombtype));
1039 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
1040 c_pbcTypeNames[ir->pbcType].c_str());
1041 CHECK(ir->eDispCorr != edispcNO);
1044 if (ir->rlist == 0.0)
1047 "can only have neighborlist cut-off zero (=infinite)\n"
1048 "with coulombtype = %s or coulombtype = %s\n"
1049 "without periodic boundary conditions (pbc = %s) and\n"
1050 "rcoulomb and rvdw set to zero",
1051 eel_names[eelCUT], eel_names[eelUSER], c_pbcTypeNames[PbcType::No].c_str());
1052 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
1053 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1055 if (ir->nstlist > 0)
1058 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1059 "nstype=simple and only one MPI rank");
1064 if (ir->nstcomm == 0)
1066 // TODO Change this behaviour. There should be exactly one way
1067 // to turn off an algorithm.
1068 ir->comm_mode = ecmNO;
1070 if (ir->comm_mode != ecmNO)
1072 if (ir->nstcomm < 0)
1074 // TODO Such input was once valid. Now that we've been
1075 // helpful for a few years, we should reject such input,
1076 // lest we have to support every historical decision
1079 "If you want to remove the rotation around the center of mass, you should set "
1080 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1081 "its absolute value");
1082 ir->nstcomm = abs(ir->nstcomm);
1085 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
1088 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
1089 "nstcomm to nstcalcenergy");
1090 ir->nstcomm = ir->nstcalcenergy;
1093 if (ir->comm_mode == ecmANGULAR)
1096 "Can not remove the rotation around the center of mass with periodic "
1098 CHECK(ir->bPeriodicMols);
1099 if (ir->pbcType != PbcType::No)
1102 "Removing the rotation around the center of mass in a periodic system, "
1103 "this can lead to artifacts. Only use this on a single (cluster of) "
1104 "molecules. This cluster should not cross periodic boundaries.");
1109 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
1112 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1113 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1116 warning_note(wi, warn_buf);
1119 /* TEMPERATURE COUPLING */
1120 if (ir->etc == etcYES)
1122 ir->etc = etcBERENDSEN;
1124 "Old option for temperature coupling given: "
1125 "changing \"yes\" to \"Berendsen\"\n");
1128 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
1130 if (ir->opts.nhchainlength < 1)
1133 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1135 ir->opts.nhchainlength);
1136 ir->opts.nhchainlength = 1;
1137 warning(wi, warn_buf);
1140 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1144 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1145 ir->opts.nhchainlength = 1;
1150 ir->opts.nhchainlength = 0;
1153 if (ir->eI == eiVVAK)
1156 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1159 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1162 if (ETC_ANDERSEN(ir->etc))
1164 sprintf(err_buf, "%s temperature control not supported for integrator %s.",
1165 etcoupl_names[ir->etc], ei_names[ir->eI]);
1166 CHECK(!(EI_VV(ir->eI)));
1168 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
1171 "Center of mass removal not necessary for %s. All velocities of coupled "
1172 "groups are rerandomized periodically, so flying ice cube errors will not "
1174 etcoupl_names[ir->etc]);
1175 warning_note(wi, warn_buf);
1179 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1180 "randomized every time step",
1181 ir->nstcomm, etcoupl_names[ir->etc]);
1182 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
1185 if (ir->etc == etcBERENDSEN)
1188 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1189 "might want to consider using the %s thermostat.",
1190 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
1191 warning_note(wi, warn_buf);
1194 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc)) && ir->epc == epcBERENDSEN)
1197 "Using Berendsen pressure coupling invalidates the "
1198 "true ensemble for the thermostat");
1199 warning(wi, warn_buf);
1202 /* PRESSURE COUPLING */
1203 if (ir->epc == epcISOTROPIC)
1205 ir->epc = epcBERENDSEN;
1207 "Old option for pressure coupling given: "
1208 "changing \"Isotropic\" to \"Berendsen\"\n");
1211 if (ir->epc != epcNO)
1213 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1215 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1216 CHECK(ir->tau_p <= 0);
1218 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1221 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1222 "times larger than nstpcouple*dt (%g)",
1223 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1224 warning(wi, warn_buf);
1228 "compressibility must be > 0 when using pressure"
1230 EPCOUPLTYPE(ir->epc));
1231 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1232 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1233 && ir->compress[ZZ][YY] <= 0));
1235 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1238 "You are generating velocities so I am assuming you "
1239 "are equilibrating a system. You are using "
1240 "%s pressure coupling, but this can be "
1241 "unstable for equilibration. If your system crashes, try "
1242 "equilibrating first with Berendsen pressure coupling. If "
1243 "you are not equilibrating the system, you can probably "
1244 "ignore this warning.",
1245 epcoupl_names[ir->epc]);
1246 warning(wi, warn_buf);
1252 if (ir->epc == epcMTTK)
1254 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1258 /* ELECTROSTATICS */
1259 /* More checks are in triple check (grompp.c) */
1261 if (ir->coulombtype == eelSWITCH)
1264 "coulombtype = %s is only for testing purposes and can lead to serious "
1265 "artifacts, advice: use coulombtype = %s",
1266 eel_names[ir->coulombtype], eel_names[eelRF_ZERO]);
1267 warning(wi, warn_buf);
1270 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1273 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1274 "format and exchanging epsilon-r and epsilon-rf",
1276 warning(wi, warn_buf);
1277 ir->epsilon_rf = ir->epsilon_r;
1278 ir->epsilon_r = 1.0;
1281 if (ir->epsilon_r == 0)
1284 "It is pointless to use long-range electrostatics with infinite relative "
1286 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1288 CHECK(EEL_FULL(ir->coulombtype));
1291 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1293 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1294 CHECK(ir->epsilon_r < 0);
1297 if (EEL_RF(ir->coulombtype))
1299 /* reaction field (at the cut-off) */
1301 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1304 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1305 eel_names[ir->coulombtype]);
1306 warning(wi, warn_buf);
1307 ir->epsilon_rf = 0.0;
1310 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1311 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1312 if (ir->epsilon_rf == ir->epsilon_r)
1314 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1315 eel_names[ir->coulombtype]);
1316 warning(wi, warn_buf);
1319 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1320 * means the interaction is zero outside rcoulomb, but it helps to
1321 * provide accurate energy conservation.
1323 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1325 if (ir_coulomb_switched(ir))
1328 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1329 "potential modifier options!",
1330 eel_names[ir->coulombtype]);
1331 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1335 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1338 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1339 "secondary coulomb-modifier.");
1340 CHECK(ir->coulomb_modifier != eintmodNONE);
1342 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1345 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1346 "secondary vdw-modifier.");
1347 CHECK(ir->vdw_modifier != eintmodNONE);
1350 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->vdwtype == evdwSWITCH
1351 || ir->vdwtype == evdwSHIFT)
1354 "The switch/shift interaction settings are just for compatibility; you will get "
1356 "performance from applying potential modifiers to your interactions!\n");
1357 warning_note(wi, warn_buf);
1360 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1362 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1364 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1366 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1367 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1368 "will be good regardless, since ewald_rtol = %g.",
1369 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1370 warning(wi, warn_buf);
1374 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1376 if (ir->rvdw_switch == 0)
1379 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1380 "potential. This suggests it was not set in the mdp, which can lead to large "
1381 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1382 "switching range.");
1383 warning(wi, warn_buf);
1387 if (EEL_FULL(ir->coulombtype))
1389 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER
1390 || ir->coulombtype == eelPMEUSERSWITCH)
1392 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1393 eel_names[ir->coulombtype]);
1394 CHECK(ir->rcoulomb > ir->rlist);
1398 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1400 // TODO: Move these checks into the ewald module with the options class
1402 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1404 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1406 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d",
1407 eel_names[ir->coulombtype], orderMin, orderMax);
1408 warning_error(wi, warn_buf);
1412 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1414 if (ir->ewald_geometry == eewg3D)
1416 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1417 c_pbcTypeNames[ir->pbcType].c_str(), eewg_names[eewg3DC]);
1418 warning(wi, warn_buf);
1420 /* This check avoids extra pbc coding for exclusion corrections */
1421 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1422 CHECK(ir->wall_ewald_zfac < 2);
1424 if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
1426 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1427 eel_names[ir->coulombtype], eewg_names[eewg3DC], c_pbcTypeNames[PbcType::XY].c_str());
1428 warning(wi, warn_buf);
1430 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1432 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1433 CHECK(ir->bPeriodicMols);
1434 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1435 warning_note(wi, warn_buf);
1437 "With epsilon_surface > 0 you can only use domain decomposition "
1438 "when there are only small molecules with all bonds constrained (mdrun will check "
1440 warning_note(wi, warn_buf);
1443 if (ir_vdw_switched(ir))
1445 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1446 CHECK(ir->rvdw_switch >= ir->rvdw);
1448 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1451 "You are applying a switch function to vdw forces or potentials from %g to %g "
1452 "nm, which is more than half the interaction range, whereas switch functions "
1453 "are intended to act only close to the cut-off.",
1454 ir->rvdw_switch, ir->rvdw);
1455 warning_note(wi, warn_buf);
1459 if (ir->vdwtype == evdwPME)
1461 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1463 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1464 evdw_names[ir->vdwtype], eintmod_names[eintmodPOTSHIFT], eintmod_names[eintmodNONE]);
1465 warning_error(wi, err_buf);
1469 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1472 "You have selected user tables with dispersion correction, the dispersion "
1473 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1474 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1475 "really want dispersion correction to -C6/r^6.");
1478 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT) && ir->rvdw != 0)
1480 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1483 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1485 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1488 /* IMPLICIT SOLVENT */
1489 if (ir->coulombtype == eelGB_NOTUSED)
1491 sprintf(warn_buf, "Invalid option %s for coulombtype", eel_names[ir->coulombtype]);
1492 warning_error(wi, warn_buf);
1497 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1502 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1506 /* interpret a number of doubles from a string and put them in an array,
1507 after allocating space for them.
1508 str = the input string
1509 n = the (pre-allocated) number of doubles read
1510 r = the output array of doubles. */
1511 static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
1513 auto values = gmx::splitString(str);
1517 for (int i = 0; i < *n; i++)
1521 (*r)[i] = gmx::fromString<real>(values[i]);
1523 catch (gmx::GromacsException&)
1525 warning_error(wi, "Invalid value " + values[i]
1526 + " in string in mdp file. Expected a real number.");
1532 static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1535 int i, j, max_n_lambda, nweights, nfep[efptNR];
1536 t_lambda* fep = ir->fepvals;
1537 t_expanded* expand = ir->expandedvals;
1538 real** count_fep_lambdas;
1539 bool bOneLambda = TRUE;
1541 snew(count_fep_lambdas, efptNR);
1543 /* FEP input processing */
1544 /* first, identify the number of lambda values for each type.
1545 All that are nonzero must have the same number */
1547 for (i = 0; i < efptNR; i++)
1549 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1552 /* now, determine the number of components. All must be either zero, or equal. */
1555 for (i = 0; i < efptNR; i++)
1557 if (nfep[i] > max_n_lambda)
1559 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1560 must have the same number if its not zero.*/
1565 for (i = 0; i < efptNR; i++)
1569 ir->fepvals->separate_dvdl[i] = FALSE;
1571 else if (nfep[i] == max_n_lambda)
1573 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute
1574 the derivative with respect to the temperature currently */
1576 ir->fepvals->separate_dvdl[i] = TRUE;
1582 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1584 nfep[i], efpt_names[i], max_n_lambda);
1587 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1588 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1590 /* the number of lambdas is the number we've read in, which is either zero
1591 or the same for all */
1592 fep->n_lambda = max_n_lambda;
1594 /* allocate space for the array of lambda values */
1595 snew(fep->all_lambda, efptNR);
1596 /* if init_lambda is defined, we need to set lambda */
1597 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1599 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1601 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1602 for (i = 0; i < efptNR; i++)
1604 snew(fep->all_lambda[i], fep->n_lambda);
1605 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1608 for (j = 0; j < fep->n_lambda; j++)
1610 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1612 sfree(count_fep_lambdas[i]);
1615 sfree(count_fep_lambdas);
1617 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1618 internal bookkeeping -- for now, init_lambda */
1620 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1622 for (i = 0; i < fep->n_lambda; i++)
1624 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1628 /* check to see if only a single component lambda is defined, and soft core is defined.
1629 In this case, turn on coulomb soft core */
1631 if (max_n_lambda == 0)
1637 for (i = 0; i < efptNR; i++)
1639 if ((nfep[i] != 0) && (i != efptFEP))
1645 if ((bOneLambda) && (fep->sc_alpha > 0))
1647 fep->bScCoul = TRUE;
1650 /* Fill in the others with the efptFEP if they are not explicitly
1651 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1652 they are all zero. */
1654 for (i = 0; i < efptNR; i++)
1656 if ((nfep[i] == 0) && (i != efptFEP))
1658 for (j = 0; j < fep->n_lambda; j++)
1660 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1666 /* now read in the weights */
1667 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1670 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1672 else if (nweights != fep->n_lambda)
1674 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1675 nweights, fep->n_lambda);
1677 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1679 expand->nstexpanded = fep->nstdhdl;
1680 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1685 static void do_simtemp_params(t_inputrec* ir)
1688 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1689 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1692 template<typename T>
1693 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1696 for (const auto& input : inputs)
1700 outputs[i] = gmx::fromStdString<T>(input);
1702 catch (gmx::GromacsException&)
1704 auto message = gmx::formatString(
1705 "Invalid value for mdp option %s. %s should only consist of integers separated "
1708 warning_error(wi, message);
1714 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1717 for (const auto& input : inputs)
1721 outputs[i] = gmx::fromString<real>(input);
1723 catch (gmx::GromacsException&)
1725 auto message = gmx::formatString(
1726 "Invalid value for mdp option %s. %s should only consist of real numbers "
1727 "separated by spaces.",
1729 warning_error(wi, message);
1735 static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
1738 for (const auto& input : inputs)
1742 outputs[i][d] = gmx::fromString<real>(input);
1744 catch (gmx::GromacsException&)
1746 auto message = gmx::formatString(
1747 "Invalid value for mdp option %s. %s should only consist of real numbers "
1748 "separated by spaces.",
1750 warning_error(wi, message);
1761 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1763 opts->wall_atomtype[0] = nullptr;
1764 opts->wall_atomtype[1] = nullptr;
1766 ir->wall_atomtype[0] = -1;
1767 ir->wall_atomtype[1] = -1;
1768 ir->wall_density[0] = 0;
1769 ir->wall_density[1] = 0;
1773 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1774 if (wallAtomTypes.size() != size_t(ir->nwall))
1776 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu", ir->nwall,
1777 wallAtomTypes.size());
1779 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1780 for (int i = 0; i < ir->nwall; i++)
1782 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1785 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1787 auto wallDensity = gmx::splitString(wall_density);
1788 if (wallDensity.size() != size_t(ir->nwall))
1790 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall,
1791 wallDensity.size());
1793 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1794 for (int i = 0; i < ir->nwall; i++)
1796 if (ir->wall_density[i] <= 0)
1798 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1805 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1809 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1810 for (int i = 0; i < nwall; i++)
1812 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1813 grps->emplace_back(groups->groupNames.size() - 1);
1818 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1820 /* read expanded ensemble parameters */
1821 printStringNewline(inp, "expanded ensemble variables");
1822 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1823 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1824 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1825 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1826 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1827 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1828 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1829 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1830 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1831 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1832 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1833 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1834 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1835 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1836 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1837 expand->bSymmetrizedTMatrix =
1838 (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1839 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1840 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1841 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1842 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1843 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1844 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1845 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1848 /*! \brief Return whether an end state with the given coupling-lambda
1849 * value describes fully-interacting VDW.
1851 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1852 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1854 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1856 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1862 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1865 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1867 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1869 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1872 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1873 std::string message = gmx::formatExceptionMessageToString(*ex);
1874 warning_error(wi_, message.c_str());
1879 std::string getOptionName(const gmx::KeyValueTreePath& context)
1881 if (mapping_ != nullptr)
1883 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1884 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1887 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1892 const gmx::IKeyValueTreeBackMapping* mapping_;
1897 void get_ir(const char* mdparin,
1898 const char* mdparout,
1899 gmx::MDModules* mdModules,
1902 WriteMdpHeader writeMdpHeader,
1906 double dumdub[2][6];
1908 char warn_buf[STRLEN];
1909 t_lambda* fep = ir->fepvals;
1910 t_expanded* expand = ir->expandedvals;
1912 const char* no_names[] = { "no", nullptr };
1914 init_inputrec_strings();
1915 gmx::TextInputFile stream(mdparin);
1916 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1918 snew(dumstr[0], STRLEN);
1919 snew(dumstr[1], STRLEN);
1921 /* ignore the following deprecated commands */
1922 replace_inp_entry(inp, "title", nullptr);
1923 replace_inp_entry(inp, "cpp", nullptr);
1924 replace_inp_entry(inp, "domain-decomposition", nullptr);
1925 replace_inp_entry(inp, "andersen-seed", nullptr);
1926 replace_inp_entry(inp, "dihre", nullptr);
1927 replace_inp_entry(inp, "dihre-fc", nullptr);
1928 replace_inp_entry(inp, "dihre-tau", nullptr);
1929 replace_inp_entry(inp, "nstdihreout", nullptr);
1930 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1931 replace_inp_entry(inp, "optimize-fft", nullptr);
1932 replace_inp_entry(inp, "adress_type", nullptr);
1933 replace_inp_entry(inp, "adress_const_wf", nullptr);
1934 replace_inp_entry(inp, "adress_ex_width", nullptr);
1935 replace_inp_entry(inp, "adress_hy_width", nullptr);
1936 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1937 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1938 replace_inp_entry(inp, "adress_site", nullptr);
1939 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1940 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1941 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1942 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1943 replace_inp_entry(inp, "rlistlong", nullptr);
1944 replace_inp_entry(inp, "nstcalclr", nullptr);
1945 replace_inp_entry(inp, "pull-print-com2", nullptr);
1946 replace_inp_entry(inp, "gb-algorithm", nullptr);
1947 replace_inp_entry(inp, "nstgbradii", nullptr);
1948 replace_inp_entry(inp, "rgbradii", nullptr);
1949 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1950 replace_inp_entry(inp, "gb-saltconc", nullptr);
1951 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1952 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1953 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1954 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1955 replace_inp_entry(inp, "sa-algorithm", nullptr);
1956 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1957 replace_inp_entry(inp, "ns-type", nullptr);
1959 /* replace the following commands with the clearer new versions*/
1960 replace_inp_entry(inp, "unconstrained-start", "continuation");
1961 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1962 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1963 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1964 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1965 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1966 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1968 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1969 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1970 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1971 setStringEntry(&inp, "include", opts->include, nullptr);
1972 printStringNoNewline(
1973 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1974 setStringEntry(&inp, "define", opts->define, nullptr);
1976 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1977 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1978 printStringNoNewline(&inp, "Start time and timestep in ps");
1979 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1980 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1981 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1982 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1983 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1984 printStringNoNewline(
1985 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1986 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1987 printStringNoNewline(&inp, "Multiple time-stepping");
1988 ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
1991 opts->numMtsLevels = get_eint(&inp, "mts-levels", 2, wi);
1992 ir->mtsLevels.resize(2);
1993 gmx::MtsLevel& mtsLevel = ir->mtsLevels[1];
1994 opts->mtsLevel2Forces = setStringEntry(&inp, "mts-level2-forces",
1995 "longrange-nonbonded nonbonded pair dihedral");
1996 mtsLevel.stepFactor = get_eint(&inp, "mts-level2-factor", 2, wi);
1998 // We clear after reading without dynamics to not force the user to remove MTS mdp options
1999 if (!EI_DYNAMICS(ir->eI))
2002 ir->mtsLevels.clear();
2005 printStringNoNewline(&inp, "mode for center of mass motion removal");
2006 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
2007 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
2008 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
2009 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
2010 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
2012 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2013 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2014 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2015 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2018 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2019 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2020 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2021 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2022 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2023 ir->niter = get_eint(&inp, "niter", 20, wi);
2024 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2025 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2026 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2027 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2028 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2030 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2031 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2033 /* Output options */
2034 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2035 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2036 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2037 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2038 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2039 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2040 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2041 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2042 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2043 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2044 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2045 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2046 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2047 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2048 printStringNoNewline(&inp, "default, all atoms will be written.");
2049 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2050 printStringNoNewline(&inp, "Selection of energy groups");
2051 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2053 /* Neighbor searching */
2054 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2055 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2056 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
2057 printStringNoNewline(&inp, "nblist update frequency");
2058 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2059 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2060 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2061 std::vector<const char*> pbcTypesNamesChar;
2062 for (const auto& pbcTypeName : c_pbcTypeNames)
2064 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2066 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2067 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
2068 printStringNoNewline(&inp,
2069 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2070 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2071 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2072 printStringNoNewline(&inp, "nblist cut-off");
2073 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2074 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2076 /* Electrostatics */
2077 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2078 printStringNoNewline(&inp, "Method for doing electrostatics");
2079 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
2080 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
2081 printStringNoNewline(&inp, "cut-off lengths");
2082 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2083 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2084 printStringNoNewline(&inp,
2085 "Relative dielectric constant for the medium and the reaction field");
2086 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2087 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2088 printStringNoNewline(&inp, "Method for doing Van der Waals");
2089 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
2090 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
2091 printStringNoNewline(&inp, "cut-off lengths");
2092 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2093 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2094 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2095 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
2096 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2097 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2098 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2099 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2100 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2101 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2102 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2103 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2104 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2105 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2106 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2107 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2108 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2109 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2110 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
2111 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
2112 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2114 /* Implicit solvation is no longer supported, but we need grompp
2115 to be able to refuse old .mdp files that would have built a tpr
2116 to run it. Thus, only "no" is accepted. */
2117 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2119 /* Coupling stuff */
2120 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2121 printStringNoNewline(&inp, "Temperature coupling");
2122 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
2123 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2124 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2125 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
2126 printStringNoNewline(&inp, "Groups to couple separately");
2127 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2128 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2129 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2130 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2131 printStringNoNewline(&inp, "pressure coupling");
2132 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
2133 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
2134 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2135 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2136 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2137 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2138 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2139 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2140 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
2143 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2144 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
2145 printStringNoNewline(&inp, "Groups treated with MiMiC");
2146 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2148 /* Simulated annealing */
2149 printStringNewline(&inp, "SIMULATED ANNEALING");
2150 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2151 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2152 printStringNoNewline(&inp,
2153 "Number of time points to use for specifying annealing in each group");
2154 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2155 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2156 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2157 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2158 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2161 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2162 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
2163 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2164 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2167 printStringNewline(&inp, "OPTIONS FOR BONDS");
2168 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2169 printStringNoNewline(&inp, "Type of constraint algorithm");
2170 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
2171 printStringNoNewline(&inp, "Do not constrain the start configuration");
2172 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
2173 printStringNoNewline(&inp,
2174 "Use successive overrelaxation to reduce the number of shake iterations");
2175 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
2176 printStringNoNewline(&inp, "Relative tolerance of shake");
2177 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2178 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2179 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2180 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2181 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2182 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2183 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2184 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2185 printStringNoNewline(&inp, "rotates over more degrees than");
2186 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2187 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2188 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
2190 /* Energy group exclusions */
2191 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2192 printStringNoNewline(
2193 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2194 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2197 printStringNewline(&inp, "WALLS");
2198 printStringNoNewline(
2199 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2200 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2201 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
2202 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2203 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2204 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2205 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2208 printStringNewline(&inp, "COM PULLING");
2209 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2213 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull, wi);
2217 for (int c = 0; c < ir->pull->ncoord; c++)
2219 if (ir->pull->coord[c].eType == epullCONSTRAINT)
2222 "Constraint COM pulling is not supported in combination with "
2223 "multiple time stepping");
2231 NOTE: needs COM pulling or free energy input */
2232 printStringNewline(&inp, "AWH biasing");
2233 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2236 ir->awhParams = gmx::readAwhParams(&inp, wi);
2239 /* Enforced rotation */
2240 printStringNewline(&inp, "ENFORCED ROTATION");
2241 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2242 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2246 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2249 /* Interactive MD */
2251 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2252 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2253 if (inputrecStrings->imd_grp[0] != '\0')
2260 printStringNewline(&inp, "NMR refinement stuff");
2261 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2262 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2263 printStringNoNewline(
2264 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2265 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2266 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2267 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2268 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2269 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2270 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2271 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2272 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2273 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2274 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2275 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2276 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2277 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2278 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2279 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2281 /* free energy variables */
2282 printStringNewline(&inp, "Free energy variables");
2283 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2284 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2285 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2286 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2287 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2289 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2291 it was not entered */
2292 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2293 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2294 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2295 setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
2296 setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
2297 setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
2298 setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
2299 setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
2300 setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
2301 setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
2302 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2303 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2304 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2305 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2306 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2307 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2308 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2309 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2310 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2311 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2312 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2313 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2314 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2315 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2317 /* Non-equilibrium MD stuff */
2318 printStringNewline(&inp, "Non-equilibrium MD stuff");
2319 setStringEntry(&inp, "acc-grps", inputrecStrings->accgrps, nullptr);
2320 setStringEntry(&inp, "accelerate", inputrecStrings->acc, nullptr);
2321 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2322 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2323 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2324 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2326 /* simulated tempering variables */
2327 printStringNewline(&inp, "simulated tempering variables");
2328 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2329 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2330 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2331 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2333 /* expanded ensemble variables */
2334 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2336 read_expandedparams(&inp, expand, wi);
2339 /* Electric fields */
2341 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2342 gmx::KeyValueTreeTransformer transform;
2343 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2344 mdModules->initMdpTransform(transform.rules());
2345 for (const auto& path : transform.mappedPaths())
2347 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2348 mark_einp_set(inp, path[0].c_str());
2350 MdpErrorHandler errorHandler(wi);
2351 auto result = transform.transform(convertedValues, &errorHandler);
2352 ir->params = new gmx::KeyValueTreeObject(result.object());
2353 mdModules->adjustInputrecBasedOnModules(ir);
2354 errorHandler.setBackMapping(result.backMapping());
2355 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2358 /* Ion/water position swapping ("computational electrophysiology") */
2359 printStringNewline(&inp,
2360 "Ion/water position swapping for computational electrophysiology setups");
2361 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2362 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2363 if (ir->eSwapCoords != eswapNO)
2370 printStringNoNewline(&inp, "Swap attempt frequency");
2371 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2372 printStringNoNewline(&inp, "Number of ion types to be controlled");
2373 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2376 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2378 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2379 snew(ir->swap->grp, ir->swap->ngrp);
2380 for (i = 0; i < ir->swap->ngrp; i++)
2382 snew(ir->swap->grp[i].molname, STRLEN);
2384 printStringNoNewline(&inp,
2385 "Two index groups that contain the compartment-partitioning atoms");
2386 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2387 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2388 printStringNoNewline(&inp,
2389 "Use center of mass of split groups (yes/no), otherwise center of "
2390 "geometry is used");
2391 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2392 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2394 printStringNoNewline(&inp, "Name of solvent molecules");
2395 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2397 printStringNoNewline(&inp,
2398 "Split cylinder: radius, upper and lower extension (nm) (this will "
2399 "define the channels)");
2400 printStringNoNewline(&inp,
2401 "Note that the split cylinder settings do not have an influence on "
2402 "the swapping protocol,");
2403 printStringNoNewline(
2405 "however, if correctly defined, the permeation events are recorded per channel");
2406 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2407 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2408 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2409 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2410 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2411 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2413 printStringNoNewline(
2415 "Average the number of ions per compartment over these many swap attempt steps");
2416 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2418 printStringNoNewline(
2419 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2420 printStringNoNewline(
2421 &inp, "and the requested number of ions of this type in compartments A and B");
2422 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2423 for (i = 0; i < nIonTypes; i++)
2425 int ig = eSwapFixedGrpNR + i;
2427 sprintf(buf, "iontype%d-name", i);
2428 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2429 sprintf(buf, "iontype%d-in-A", i);
2430 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2431 sprintf(buf, "iontype%d-in-B", i);
2432 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2435 printStringNoNewline(
2437 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2438 printStringNoNewline(
2440 "at maximum distance (= bulk concentration) to the split group layers. However,");
2441 printStringNoNewline(&inp,
2442 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2443 "layer from the middle at 0.0");
2444 printStringNoNewline(&inp,
2445 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2446 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2447 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2448 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2449 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2451 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2454 printStringNoNewline(
2455 &inp, "Start to swap ions if threshold difference to requested count is reached");
2456 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2459 /* AdResS is no longer supported, but we need grompp to be able to
2460 refuse to process old .mdp files that used it. */
2461 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2463 /* User defined thingies */
2464 printStringNewline(&inp, "User defined thingies");
2465 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2466 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2467 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2468 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2469 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2470 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2471 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2472 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2473 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2474 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2478 gmx::TextOutputFile stream(mdparout);
2479 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2481 // Transform module data into a flat key-value tree for output.
2482 gmx::KeyValueTreeBuilder builder;
2483 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2484 mdModules->buildMdpOutput(&builderObject);
2486 gmx::TextWriter writer(&stream);
2487 writeKeyValueTreeAsMdp(&writer, builder.build());
2492 /* Process options if necessary */
2493 for (m = 0; m < 2; m++)
2495 for (i = 0; i < 2 * DIM; i++)
2504 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2508 "Pressure coupling incorrect number of values (I need exactly 1)");
2510 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2512 case epctSEMIISOTROPIC:
2513 case epctSURFACETENSION:
2514 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2518 "Pressure coupling incorrect number of values (I need exactly 2)");
2520 dumdub[m][YY] = dumdub[m][XX];
2522 case epctANISOTROPIC:
2523 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf", &(dumdub[m][XX]), &(dumdub[m][YY]),
2524 &(dumdub[m][ZZ]), &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5]))
2529 "Pressure coupling incorrect number of values (I need exactly 6)");
2533 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2534 epcoupltype_names[ir->epct]);
2538 clear_mat(ir->ref_p);
2539 clear_mat(ir->compress);
2540 for (i = 0; i < DIM; i++)
2542 ir->ref_p[i][i] = dumdub[1][i];
2543 ir->compress[i][i] = dumdub[0][i];
2545 if (ir->epct == epctANISOTROPIC)
2547 ir->ref_p[XX][YY] = dumdub[1][3];
2548 ir->ref_p[XX][ZZ] = dumdub[1][4];
2549 ir->ref_p[YY][ZZ] = dumdub[1][5];
2550 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2553 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2554 "apply a threefold shear stress?\n");
2556 ir->compress[XX][YY] = dumdub[0][3];
2557 ir->compress[XX][ZZ] = dumdub[0][4];
2558 ir->compress[YY][ZZ] = dumdub[0][5];
2559 for (i = 0; i < DIM; i++)
2561 for (m = 0; m < i; m++)
2563 ir->ref_p[i][m] = ir->ref_p[m][i];
2564 ir->compress[i][m] = ir->compress[m][i];
2569 if (ir->comm_mode == ecmNO)
2574 opts->couple_moltype = nullptr;
2575 if (strlen(inputrecStrings->couple_moltype) > 0)
2577 if (ir->efep != efepNO)
2579 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2580 if (opts->couple_lam0 == opts->couple_lam1)
2582 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2584 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2588 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2595 "Free energy is turned off, so we will not decouple the molecule listed "
2599 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2600 if (ir->efep != efepNO)
2602 if (fep->delta_lambda != 0)
2604 ir->efep = efepSLOWGROWTH;
2608 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2610 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2612 "Old option for dhdl-print-energy given: "
2613 "changing \"yes\" to \"total\"\n");
2616 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2618 /* always print out the energy to dhdl if we are doing
2619 expanded ensemble, since we need the total energy for
2620 analysis if the temperature is changing. In some
2621 conditions one may only want the potential energy, so
2622 we will allow that if the appropriate mdp setting has
2623 been enabled. Otherwise, total it is:
2625 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2628 if ((ir->efep != efepNO) || ir->bSimTemp)
2630 ir->bExpanded = FALSE;
2631 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2633 ir->bExpanded = TRUE;
2635 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2636 if (ir->bSimTemp) /* done after fep params */
2638 do_simtemp_params(ir);
2641 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2642 * setup and not on the old way of specifying the free-energy setup,
2643 * we should check for using soft-core when not needed, since that
2644 * can complicate the sampling significantly.
2645 * Note that we only check for the automated coupling setup.
2646 * If the (advanced) user does FEP through manual topology changes,
2647 * this check will not be triggered.
2649 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 && ir->fepvals->sc_alpha != 0
2650 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2653 "You are using soft-core interactions while the Van der Waals interactions are "
2654 "not decoupled (note that the sc-coul option is only active when using lambda "
2655 "states). Although this will not lead to errors, you will need much more "
2656 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2661 ir->fepvals->n_lambda = 0;
2664 /* WALL PARAMETERS */
2666 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2668 /* ORIENTATION RESTRAINT PARAMETERS */
2670 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2672 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2675 /* DEFORMATION PARAMETERS */
2677 clear_mat(ir->deform);
2678 for (i = 0; i < 6; i++)
2683 double gmx_unused canary;
2684 int ndeform = sscanf(inputrecStrings->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]),
2685 &(dumdub[0][1]), &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]),
2686 &(dumdub[0][5]), &canary);
2688 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2692 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2693 inputrecStrings->deform)
2696 for (i = 0; i < 3; i++)
2698 ir->deform[i][i] = dumdub[0][i];
2700 ir->deform[YY][XX] = dumdub[0][3];
2701 ir->deform[ZZ][XX] = dumdub[0][4];
2702 ir->deform[ZZ][YY] = dumdub[0][5];
2703 if (ir->epc != epcNO)
2705 for (i = 0; i < 3; i++)
2707 for (j = 0; j <= i; j++)
2709 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2711 warning_error(wi, "A box element has deform set and compressibility > 0");
2715 for (i = 0; i < 3; i++)
2717 for (j = 0; j < i; j++)
2719 if (ir->deform[i][j] != 0)
2721 for (m = j; m < DIM; m++)
2723 if (ir->compress[m][j] != 0)
2726 "An off-diagonal box element has deform set while "
2727 "compressibility > 0 for the same component of another box "
2728 "vector, this might lead to spurious periodicity effects.");
2729 warning(wi, warn_buf);
2737 /* Ion/water position swapping checks */
2738 if (ir->eSwapCoords != eswapNO)
2740 if (ir->swap->nstswap < 1)
2742 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2744 if (ir->swap->nAverage < 1)
2746 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2748 if (ir->swap->threshold < 1.0)
2750 warning_error(wi, "Ion count threshold must be at least 1.\n");
2754 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2757 setupMtsLevels(ir->mtsLevels, *ir, *opts, wi);
2762 gmx::checkAwhParams(ir->awhParams, ir, wi);
2769 /* We would like gn to be const as well, but C doesn't allow this */
2770 /* TODO this is utility functionality (search for the index of a
2771 string in a collection), so should be refactored and located more
2773 int search_string(const char* s, int ng, char* gn[])
2777 for (i = 0; (i < ng); i++)
2779 if (gmx_strcasecmp(s, gn[i]) == 0)
2786 "Group %s referenced in the .mdp file was not found in the index file.\n"
2787 "Group names must match either [moleculetype] names or custom index group\n"
2788 "names, in which case you must supply an index file to the '-n' option\n"
2793 static void do_numbering(int natoms,
2794 SimulationGroups* groups,
2795 gmx::ArrayRef<std::string> groupsFromMdpFile,
2798 SimulationAtomGroupType gtype,
2804 unsigned short* cbuf;
2805 AtomGroupIndices* grps = &(groups->groups[gtype]);
2806 int j, gid, aj, ognr, ntot = 0;
2808 char warn_buf[STRLEN];
2810 title = shortName(gtype);
2813 /* Mark all id's as not set */
2814 for (int i = 0; (i < natoms); i++)
2819 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2821 /* Lookup the group name in the block structure */
2822 gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2823 if ((grptp != egrptpONE) || (i == 0))
2825 grps->emplace_back(gid);
2828 /* Now go over the atoms in the group */
2829 for (j = block->index[gid]; (j < block->index[gid + 1]); j++)
2834 /* Range checking */
2835 if ((aj < 0) || (aj >= natoms))
2837 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2839 /* Lookup up the old group number */
2843 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
2848 /* Store the group number in buffer */
2849 if (grptp == egrptpONE)
2862 /* Now check whether we have done all atoms */
2865 if (grptp == egrptpALL)
2867 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2869 else if (grptp == egrptpPART)
2871 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2872 warning_note(wi, warn_buf);
2874 /* Assign all atoms currently unassigned to a rest group */
2875 for (j = 0; (j < natoms); j++)
2877 if (cbuf[j] == NOGID)
2879 cbuf[j] = grps->size();
2882 if (grptp != egrptpPART)
2886 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title,
2889 /* Add group name "rest" */
2890 grps->emplace_back(restnm);
2892 /* Assign the rest name to all atoms not currently assigned to a group */
2893 for (j = 0; (j < natoms); j++)
2895 if (cbuf[j] == NOGID)
2897 // group size was not updated before this here, so need to use -1.
2898 cbuf[j] = grps->size() - 1;
2904 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2906 /* All atoms are part of one (or no) group, no index required */
2907 groups->groupNumbers[gtype].clear();
2911 for (int j = 0; (j < natoms); j++)
2913 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2920 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2923 pull_params_t* pull;
2924 int natoms, imin, jmin;
2925 int * nrdf2, *na_vcm, na_tot;
2926 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2931 * First calc 3xnr-atoms for each group
2932 * then subtract half a degree of freedom for each constraint
2934 * Only atoms and nuclei contribute to the degrees of freedom...
2939 const SimulationGroups& groups = mtop->groups;
2940 natoms = mtop->natoms;
2942 /* Allocate one more for a possible rest group */
2943 /* We need to sum degrees of freedom into doubles,
2944 * since floats give too low nrdf's above 3 million atoms.
2946 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
2947 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2948 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2949 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2950 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2952 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2956 for (gmx::index i = 0;
2957 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; i++)
2960 clear_ivec(dof_vcm[i]);
2962 nrdf_vcm_sub[i] = 0;
2964 snew(nrdf2, natoms);
2965 for (const AtomProxy atomP : AtomRange(*mtop))
2967 const t_atom& local = atomP.atom();
2968 int i = atomP.globalAtomNumber();
2970 if (local.ptype == eptAtom || local.ptype == eptNucleus)
2972 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
2973 for (int d = 0; d < DIM; d++)
2975 if (opts->nFreeze[g][d] == 0)
2977 /* Add one DOF for particle i (counted as 2*1) */
2979 /* VCM group i has dim d as a DOF */
2980 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
2984 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
2986 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
2992 for (const gmx_molblock_t& molb : mtop->molblock)
2994 const gmx_moltype_t& molt = mtop->moltype[molb.type];
2995 const t_atom* atom = molt.atoms.atom;
2996 for (int mol = 0; mol < molb.nmol; mol++)
2998 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
3000 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
3001 for (int i = 0; i < molt.ilist[ftype].size();)
3003 /* Subtract degrees of freedom for the constraints,
3004 * if the particles still have degrees of freedom left.
3005 * If one of the particles is a vsite or a shell, then all
3006 * constraint motion will go there, but since they do not
3007 * contribute to the constraints the degrees of freedom do not
3010 int ai = as + ia[i + 1];
3011 int aj = as + ia[i + 2];
3012 if (((atom[ia[i + 1]].ptype == eptNucleus) || (atom[ia[i + 1]].ptype == eptAtom))
3013 && ((atom[ia[i + 2]].ptype == eptNucleus) || (atom[ia[i + 2]].ptype == eptAtom)))
3031 imin = std::min(imin, nrdf2[ai]);
3032 jmin = std::min(jmin, nrdf2[aj]);
3035 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3037 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3039 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3041 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3044 i += interaction_function[ftype].nratoms + 1;
3047 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3048 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3050 /* Subtract 1 dof from every atom in the SETTLE */
3051 for (int j = 0; j < 3; j++)
3053 int ai = as + ia[i + 1 + j];
3054 imin = std::min(2, nrdf2[ai]);
3056 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3058 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3063 as += molt.atoms.nr;
3069 /* Correct nrdf for the COM constraints.
3070 * We correct using the TC and VCM group of the first atom
3071 * in the reference and pull group. If atoms in one pull group
3072 * belong to different TC or VCM groups it is anyhow difficult
3073 * to determine the optimal nrdf assignment.
3077 for (int i = 0; i < pull->ncoord; i++)
3079 if (pull->coord[i].eType != epullCONSTRAINT)
3086 for (int j = 0; j < 2; j++)
3088 const t_pull_group* pgrp;
3090 pgrp = &pull->group[pull->coord[i].group[j]];
3094 /* Subtract 1/2 dof from each group */
3095 int ai = pgrp->ind[0];
3096 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3098 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3100 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3103 "Center of mass pulling constraints caused the number of degrees "
3104 "of freedom for temperature coupling group %s to be negative",
3105 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3106 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3111 /* We need to subtract the whole DOF from group j=1 */
3118 if (ir->nstcomm != 0)
3120 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3121 "Expect at least one group when removing COM motion");
3123 /* We remove COM motion up to dim ndof_com() */
3124 const int ndim_rm_vcm = ndof_com(ir);
3126 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3127 * the number of degrees of freedom in each vcm group when COM
3128 * translation is removed and 6 when rotation is removed as well.
3129 * Note that we do not and should not include the rest group here.
3131 for (gmx::index j = 0;
3132 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]); j++)
3134 switch (ir->comm_mode)
3137 case ecmLINEAR_ACCELERATION_CORRECTION:
3138 nrdf_vcm_sub[j] = 0;
3139 for (int d = 0; d < ndim_rm_vcm; d++)
3147 case ecmANGULAR: nrdf_vcm_sub[j] = 6; break;
3148 default: gmx_incons("Checking comm_mode");
3152 for (gmx::index i = 0;
3153 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
3155 /* Count the number of atoms of TC group i for every VCM group */
3156 for (gmx::index j = 0;
3157 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3162 for (int ai = 0; ai < natoms; ai++)
3164 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3166 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3170 /* Correct for VCM removal according to the fraction of each VCM
3171 * group present in this TC group.
3173 nrdf_uc = nrdf_tc[i];
3175 for (gmx::index j = 0;
3176 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3178 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3180 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3181 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3186 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3188 opts->nrdf[i] = nrdf_tc[i];
3189 if (opts->nrdf[i] < 0)
3193 fprintf(stderr, "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3194 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
3202 sfree(nrdf_vcm_sub);
3205 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3207 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3208 * But since this is much larger than STRLEN, such a line can not be parsed.
3209 * The real maximum is the number of names that fit in a string: STRLEN/2.
3211 #define EGP_MAX (STRLEN / 2)
3215 auto names = gmx::splitString(val);
3216 if (names.size() % 2 != 0)
3218 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3220 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3222 for (size_t i = 0; i < names.size() / 2; i++)
3224 // TODO this needs to be replaced by a solution using std::find_if
3228 names[2 * i].c_str(),
3229 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3235 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3240 names[2 * i + 1].c_str(),
3241 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3247 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3249 if ((j < nr) && (k < nr))
3251 ir->opts.egp_flags[nr * j + k] |= flag;
3252 ir->opts.egp_flags[nr * k + j] |= flag;
3261 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3263 int ig = -1, i = 0, gind;
3267 /* Just a quick check here, more thorough checks are in mdrun */
3268 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3270 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3273 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3274 for (ig = 0; ig < swap->ngrp; ig++)
3276 swapg = &swap->grp[ig];
3277 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3278 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3282 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3283 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap", swap->grp[ig].molname, swapg->nat);
3284 snew(swapg->ind, swapg->nat);
3285 for (i = 0; i < swapg->nat; i++)
3287 swapg->ind[i] = grps->a[grps->index[gind] + i];
3292 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3298 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3303 ig = search_string(IMDgname, grps->nr, gnames);
3304 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3306 if (IMDgroup->nat > 0)
3309 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3311 IMDgname, IMDgroup->nat);
3312 snew(IMDgroup->ind, IMDgroup->nat);
3313 for (i = 0; i < IMDgroup->nat; i++)
3315 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3320 /* Checks whether atoms are both part of a COM removal group and frozen.
3321 * If a fully frozen atom is part of a COM removal group, it is removed
3322 * from the COM removal group. A note is issued if such atoms are present.
3323 * A warning is issued for atom with one or two dimensions frozen that
3324 * are part of a COM removal group (mdrun would need to compute COM mass
3325 * per dimension to handle this correctly).
3326 * Also issues a warning when non-frozen atoms are not part of a COM
3327 * removal group while COM removal is active.
3329 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3331 const t_grpopts& opts,
3334 const int vcmRestGroup =
3335 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3337 int numFullyFrozenVcmAtoms = 0;
3338 int numPartiallyFrozenVcmAtoms = 0;
3339 int numNonVcmAtoms = 0;
3340 for (int a = 0; a < numAtoms; a++)
3342 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3343 int numFrozenDims = 0;
3344 for (int d = 0; d < DIM; d++)
3346 numFrozenDims += opts.nFreeze[freezeGroup][d];
3349 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3350 if (vcmGroup < vcmRestGroup)
3352 if (numFrozenDims == DIM)
3354 /* Do not remove COM motion for this fully frozen atom */
3355 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3357 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3360 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3361 numFullyFrozenVcmAtoms++;
3363 else if (numFrozenDims > 0)
3365 numPartiallyFrozenVcmAtoms++;
3368 else if (numFrozenDims < DIM)
3374 if (numFullyFrozenVcmAtoms > 0)
3376 std::string warningText = gmx::formatString(
3377 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3378 "removing these atoms from the COMM removal group(s)",
3379 numFullyFrozenVcmAtoms);
3380 warning_note(wi, warningText.c_str());
3382 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3384 std::string warningText = gmx::formatString(
3385 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3386 "removal group(s), due to limitations in the code these still contribute to the "
3387 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3389 numPartiallyFrozenVcmAtoms, DIM);
3390 warning(wi, warningText.c_str());
3392 if (numNonVcmAtoms > 0)
3394 std::string warningText = gmx::formatString(
3395 "%d atoms are not part of any center of mass motion removal group.\n"
3396 "This may lead to artifacts.\n"
3397 "In most cases one should use one group for the whole system.",
3399 warning(wi, warningText.c_str());
3403 void do_index(const char* mdparin,
3407 const gmx::MdModulesNotifier& notifier,
3411 t_blocka* defaultIndexGroups;
3419 int i, j, k, restnm;
3420 bool bExcl, bTable, bAnneal;
3421 char warn_buf[STRLEN];
3425 fprintf(stderr, "processing index file...\n");
3429 snew(defaultIndexGroups, 1);
3430 snew(defaultIndexGroups->index, 1);
3432 atoms_all = gmx_mtop_global_atoms(mtop);
3433 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3434 done_atom(&atoms_all);
3438 defaultIndexGroups = init_index(ndx, &gnames);
3441 SimulationGroups* groups = &mtop->groups;
3442 natoms = mtop->natoms;
3443 symtab = &mtop->symtab;
3445 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3447 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3449 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3450 restnm = groups->groupNames.size() - 1;
3451 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3452 srenew(gnames, defaultIndexGroups->nr + 1);
3453 gnames[restnm] = *(groups->groupNames.back());
3455 set_warning_line(wi, mdparin, -1);
3457 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3458 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3459 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3460 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3461 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3464 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3466 temperatureCouplingGroupNames.size(), temperatureCouplingReferenceValues.size(),
3467 temperatureCouplingTauValues.size());
3470 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3471 do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
3472 SimulationAtomGroupType::TemperatureCoupling, restnm,
3473 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3474 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3476 snew(ir->opts.nrdf, nr);
3477 snew(ir->opts.tau_t, nr);
3478 snew(ir->opts.ref_t, nr);
3479 if (ir->eI == eiBD && ir->bd_fric == 0)
3481 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3484 if (useReferenceTemperature)
3486 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3488 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3492 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3493 for (i = 0; (i < nr); i++)
3495 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3497 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3498 warning_error(wi, warn_buf);
3501 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3505 "tau-t = -1 is the value to signal that a group should not have "
3506 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3509 if (ir->opts.tau_t[i] >= 0)
3511 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3514 if (ir->etc != etcNO && ir->nsttcouple == -1)
3516 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3521 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3524 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3525 "md-vv; use either vrescale temperature with berendsen pressure or "
3526 "Nose-Hoover temperature with MTTK pressure");
3528 if (ir->epc == epcMTTK)
3530 if (ir->etc != etcNOSEHOOVER)
3533 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3538 if (ir->nstpcouple != ir->nsttcouple)
3540 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3541 ir->nstpcouple = ir->nsttcouple = mincouple;
3543 "for current Trotter decomposition methods with vv, nsttcouple and "
3544 "nstpcouple must be equal. Both have been reset to "
3545 "min(nsttcouple,nstpcouple) = %d",
3547 warning_note(wi, warn_buf);
3552 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3553 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3555 if (ETC_ANDERSEN(ir->etc))
3557 if (ir->nsttcouple != 1)
3561 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3562 "need for larger nsttcouple > 1, since no global parameters are computed. "
3563 "nsttcouple has been reset to 1");
3564 warning_note(wi, warn_buf);
3567 nstcmin = tcouple_min_integration_steps(ir->etc);
3570 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3573 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3574 "least %d times larger than nsttcouple*dt (%g)",
3575 ETCOUPLTYPE(ir->etc), tau_min, nstcmin, ir->nsttcouple * ir->delta_t);
3576 warning(wi, warn_buf);
3579 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3580 for (i = 0; (i < nr); i++)
3582 if (ir->opts.ref_t[i] < 0)
3584 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3587 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3588 if we are in this conditional) if mc_temp is negative */
3589 if (ir->expandedvals->mc_temp < 0)
3591 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3595 /* Simulated annealing for each group. There are nr groups */
3596 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3597 if (simulatedAnnealingGroupNames.size() == 1
3598 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3600 simulatedAnnealingGroupNames.resize(0);
3602 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3604 gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
3605 simulatedAnnealingGroupNames.size(), nr);
3609 snew(ir->opts.annealing, nr);
3610 snew(ir->opts.anneal_npoints, nr);
3611 snew(ir->opts.anneal_time, nr);
3612 snew(ir->opts.anneal_temp, nr);
3613 for (i = 0; i < nr; i++)
3615 ir->opts.annealing[i] = eannNO;
3616 ir->opts.anneal_npoints[i] = 0;
3617 ir->opts.anneal_time[i] = nullptr;
3618 ir->opts.anneal_temp[i] = nullptr;
3620 if (!simulatedAnnealingGroupNames.empty())
3623 for (i = 0; i < nr; i++)
3625 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3627 ir->opts.annealing[i] = eannNO;
3629 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3631 ir->opts.annealing[i] = eannSINGLE;
3634 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3636 ir->opts.annealing[i] = eannPERIODIC;
3642 /* Read the other fields too */
3643 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3644 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3646 gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
3647 simulatedAnnealingPoints.size(), simulatedAnnealingGroupNames.size());
3649 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3650 size_t numSimulatedAnnealingFields = 0;
3651 for (i = 0; i < nr; i++)
3653 if (ir->opts.anneal_npoints[i] == 1)
3657 "Please specify at least a start and an end point for annealing\n");
3659 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3660 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3661 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3664 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3666 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3668 gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
3669 simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
3671 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3672 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3674 gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
3675 simulatedAnnealingTemperatures.size(), numSimulatedAnnealingFields);
3678 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3679 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3680 convertReals(wi, simulatedAnnealingTimes, "anneal-time",
3681 allSimulatedAnnealingTimes.data());
3682 convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp",
3683 allSimulatedAnnealingTemperatures.data());
3684 for (i = 0, k = 0; i < nr; i++)
3686 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3688 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3689 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3692 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3694 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3700 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3703 "Annealing timepoints out of order: t=%f comes after "
3705 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j - 1]);
3708 if (ir->opts.anneal_temp[i][j] < 0)
3710 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n",
3711 ir->opts.anneal_temp[i][j]);
3716 /* Print out some summary information, to make sure we got it right */
3717 for (i = 0; i < nr; i++)
3719 if (ir->opts.annealing[i] != eannNO)
3721 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3722 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3723 *(groups->groupNames[j]), eann_names[ir->opts.annealing[i]],
3724 ir->opts.anneal_npoints[i]);
3725 fprintf(stderr, "Time (ps) Temperature (K)\n");
3726 /* All terms except the last one */
3727 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3729 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3730 ir->opts.anneal_temp[i][j]);
3733 /* Finally the last one */
3734 j = ir->opts.anneal_npoints[i] - 1;
3735 if (ir->opts.annealing[i] == eannSINGLE)
3737 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j],
3738 ir->opts.anneal_temp[i][j]);
3742 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3743 ir->opts.anneal_temp[i][j]);
3744 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3747 "There is a temperature jump when your annealing "
3759 make_pull_groups(ir->pull, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3761 make_pull_coords(ir->pull);
3766 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3769 if (ir->eSwapCoords != eswapNO)
3771 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3774 /* Make indices for IMD session */
3777 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3780 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3781 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3782 notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
3784 auto accelerations = gmx::splitString(inputrecStrings->acc);
3785 auto accelerationGroupNames = gmx::splitString(inputrecStrings->accgrps);
3786 if (accelerationGroupNames.size() * DIM != accelerations.size())
3788 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3789 accelerationGroupNames.size(), accelerations.size());
3791 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3792 SimulationAtomGroupType::Acceleration, restnm, egrptpALL_GENREST, bVerbose, wi);
3793 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3794 snew(ir->opts.acc, nr);
3795 ir->opts.ngacc = nr;
3797 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3799 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3800 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3801 if (freezeDims.size() != DIM * freezeGroupNames.size())
3803 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3804 freezeGroupNames.size(), freezeDims.size());
3806 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3807 SimulationAtomGroupType::Freeze, restnm, egrptpALL_GENREST, bVerbose, wi);
3808 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3809 ir->opts.ngfrz = nr;
3810 snew(ir->opts.nFreeze, nr);
3811 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3813 for (j = 0; (j < DIM); j++, k++)
3815 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3816 if (!ir->opts.nFreeze[i][j])
3818 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3821 "Please use Y(ES) or N(O) for freezedim only "
3823 freezeDims[k].c_str());
3824 warning(wi, warn_buf);
3829 for (; (i < nr); i++)
3831 for (j = 0; (j < DIM); j++)
3833 ir->opts.nFreeze[i][j] = 0;
3837 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3838 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3839 SimulationAtomGroupType::EnergyOutput, restnm, egrptpALL_GENREST, bVerbose, wi);
3840 add_wall_energrps(groups, ir->nwall, symtab);
3841 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3842 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3843 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3844 SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
3845 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3847 if (ir->comm_mode != ecmNO)
3849 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3852 /* Now we have filled the freeze struct, so we can calculate NRDF */
3853 calc_nrdf(mtop, ir, gnames);
3855 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3856 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3857 SimulationAtomGroupType::User1, restnm, egrptpALL_GENREST, bVerbose, wi);
3858 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3859 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3860 SimulationAtomGroupType::User2, restnm, egrptpALL_GENREST, bVerbose, wi);
3861 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3862 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3863 SimulationAtomGroupType::CompressedPositionOutput, restnm, egrptpONE, bVerbose, wi);
3864 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3865 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3866 SimulationAtomGroupType::OrientationRestraintsFit, restnm, egrptpALL_GENREST,
3869 /* MiMiC QMMM input processing */
3870 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
3871 if (qmGroupNames.size() > 1)
3873 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3875 /* group rest, if any, is always MM! */
3876 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3877 SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
3878 ir->opts.ngQM = qmGroupNames.size();
3880 /* end of MiMiC QMMM input */
3884 for (auto group : gmx::keysOf(groups->groups))
3886 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3887 for (const auto& entry : groups->groups[group])
3889 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3891 fprintf(stderr, "\n");
3895 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3896 snew(ir->opts.egp_flags, nr * nr);
3898 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
3899 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3901 warning_error(wi, "Energy group exclusions are currently not supported");
3903 if (bExcl && EEL_FULL(ir->coulombtype))
3905 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3908 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
3909 if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
3910 && !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
3913 "Can only have energy group pair tables in combination with user tables for VdW "
3917 /* final check before going out of scope if simulated tempering variables
3918 * need to be set to default values.
3920 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3922 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
3923 warning(wi, gmx::formatString(
3924 "the value for nstexpanded was not specified for "
3925 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3926 "by default, but it is recommended to set it to an explicit value!",
3927 ir->expandedvals->nstexpanded));
3929 for (i = 0; (i < defaultIndexGroups->nr); i++)
3934 done_blocka(defaultIndexGroups);
3935 sfree(defaultIndexGroups);
3939 static void check_disre(const gmx_mtop_t* mtop)
3941 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3943 const gmx_ffparams_t& ffparams = mtop->ffparams;
3946 for (int i = 0; i < ffparams.numTypes(); i++)
3948 int ftype = ffparams.functype[i];
3949 if (ftype == F_DISRES)
3951 int label = ffparams.iparams[i].disres.label;
3952 if (label == old_label)
3954 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3963 "Found %d double distance restraint indices,\n"
3964 "probably the parameters for multiple pairs in one restraint "
3965 "are not identical\n",
3971 static bool absolute_reference(const t_inputrec* ir, const gmx_mtop_t* sys, const bool posres_only, ivec AbsRef)
3974 gmx_mtop_ilistloop_t iloop;
3976 const t_iparams* pr;
3983 for (d = 0; d < DIM; d++)
3985 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3987 /* Check for freeze groups */
3988 for (g = 0; g < ir->opts.ngfrz; g++)
3990 for (d = 0; d < DIM; d++)
3992 if (ir->opts.nFreeze[g][d] != 0)
4000 /* Check for position restraints */
4001 iloop = gmx_mtop_ilistloop_init(sys);
4002 while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
4004 if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
4006 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
4008 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
4009 for (d = 0; d < DIM; d++)
4011 if (pr->posres.fcA[d] != 0)
4017 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
4019 /* Check for flat-bottom posres */
4020 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
4021 if (pr->fbposres.k != 0)
4023 switch (pr->fbposres.geom)
4025 case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
4026 case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
4027 case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
4028 case efbposresCYLINDER:
4029 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4030 case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
4031 case efbposresX: /* d=XX */
4032 case efbposresY: /* d=YY */
4033 case efbposresZ: /* d=ZZ */
4034 d = pr->fbposres.geom - efbposresX;
4039 " Invalid geometry for flat-bottom position restraint.\n"
4040 "Expected nr between 1 and %d. Found %d\n",
4041 efbposresNR - 1, pr->fbposres.geom);
4048 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
4051 static void check_combination_rule_differences(const gmx_mtop_t* mtop,
4053 bool* bC6ParametersWorkWithGeometricRules,
4054 bool* bC6ParametersWorkWithLBRules,
4055 bool* bLBRulesPossible)
4057 int ntypes, tpi, tpj;
4060 double c6i, c6j, c12i, c12j;
4061 double c6, c6_geometric, c6_LB;
4062 double sigmai, sigmaj, epsi, epsj;
4063 bool bCanDoLBRules, bCanDoGeometricRules;
4066 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4067 * force-field floating point parameters.
4070 ptr = getenv("GMX_LJCOMB_TOL");
4074 double gmx_unused canary;
4076 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4079 "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4084 *bC6ParametersWorkWithLBRules = TRUE;
4085 *bC6ParametersWorkWithGeometricRules = TRUE;
4086 bCanDoLBRules = TRUE;
4087 ntypes = mtop->ffparams.atnr;
4088 snew(typecount, ntypes);
4089 gmx_mtop_count_atomtypes(mtop, state, typecount);
4090 *bLBRulesPossible = TRUE;
4091 for (tpi = 0; tpi < ntypes; ++tpi)
4093 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4094 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4095 for (tpj = tpi; tpj < ntypes; ++tpj)
4097 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4098 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4099 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4100 c6_geometric = std::sqrt(c6i * c6j);
4101 if (!gmx_numzero(c6_geometric))
4103 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4105 sigmai = gmx::sixthroot(c12i / c6i);
4106 sigmaj = gmx::sixthroot(c12j / c6j);
4107 epsi = c6i * c6i / (4.0 * c12i);
4108 epsj = c6j * c6j / (4.0 * c12j);
4109 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4113 *bLBRulesPossible = FALSE;
4114 c6_LB = c6_geometric;
4116 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4121 *bC6ParametersWorkWithLBRules = FALSE;
4124 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4126 if (!bCanDoGeometricRules)
4128 *bC6ParametersWorkWithGeometricRules = FALSE;
4135 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop, warninp_t wi)
4137 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4139 check_combination_rule_differences(mtop, 0, &bC6ParametersWorkWithGeometricRules,
4140 &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4141 if (ir->ljpme_combination_rule == eljpmeLB)
4143 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4146 "You are using arithmetic-geometric combination rules "
4147 "in LJ-PME, but your non-bonded C6 parameters do not "
4148 "follow these rules.");
4153 if (!bC6ParametersWorkWithGeometricRules)
4155 if (ir->eDispCorr != edispcNO)
4158 "You are using geometric combination rules in "
4159 "LJ-PME, but your non-bonded C6 parameters do "
4160 "not follow these rules. "
4161 "This will introduce very small errors in the forces and energies in "
4162 "your simulations. Dispersion correction will correct total energy "
4163 "and/or pressure for isotropic systems, but not forces or surface "
4169 "You are using geometric combination rules in "
4170 "LJ-PME, but your non-bonded C6 parameters do "
4171 "not follow these rules. "
4172 "This will introduce very small errors in the forces and energies in "
4173 "your simulations. If your system is homogeneous, consider using "
4174 "dispersion correction "
4175 "for the total energy and pressure.");
4181 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4183 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4184 gmx::assertMtsRequirements(*ir);
4186 char err_buf[STRLEN];
4191 gmx_mtop_atomloop_block_t aloopb;
4193 char warn_buf[STRLEN];
4195 set_warning_line(wi, mdparin, -1);
4197 if (absolute_reference(ir, sys, false, AbsRef))
4200 "Removing center of mass motion in the presence of position restraints might "
4201 "cause artifacts. When you are using position restraints to equilibrate a "
4202 "macro-molecule, the artifacts are usually negligible.");
4205 if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 && ir->nstlist > 1
4206 && ((EI_MD(ir->eI) || EI_SD(ir->eI)) && (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4208 /* Check if a too small Verlet buffer might potentially
4209 * cause more drift than the thermostat can couple off.
4211 /* Temperature error fraction for warning and suggestion */
4212 const real T_error_warn = 0.002;
4213 const real T_error_suggest = 0.001;
4214 /* For safety: 2 DOF per atom (typical with constraints) */
4215 const real nrdf_at = 2;
4216 real T, tau, max_T_error;
4221 for (i = 0; i < ir->opts.ngtc; i++)
4223 T = std::max(T, ir->opts.ref_t[i]);
4224 tau = std::max(tau, ir->opts.tau_t[i]);
4228 /* This is a worst case estimate of the temperature error,
4229 * assuming perfect buffer estimation and no cancelation
4230 * of errors. The factor 0.5 is because energy distributes
4231 * equally over Ekin and Epot.
4233 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * BOLTZ * T);
4234 if (max_T_error > T_error_warn)
4237 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4238 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4239 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4240 "%.0e or decrease tau_t.",
4241 ir->verletbuf_tol, T, tau, 100 * max_T_error, 100 * T_error_suggest,
4242 ir->verletbuf_tol * T_error_suggest / max_T_error);
4243 warning(wi, warn_buf);
4248 if (ETC_ANDERSEN(ir->etc))
4252 for (i = 0; i < ir->opts.ngtc; i++)
4255 "all tau_t must currently be equal using Andersen temperature control, "
4256 "violated for group %d",
4258 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4260 "all tau_t must be positive using Andersen temperature control, "
4262 i, ir->opts.tau_t[i]);
4263 CHECK(ir->opts.tau_t[i] < 0);
4266 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
4268 for (i = 0; i < ir->opts.ngtc; i++)
4270 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4272 "tau_t/delta_t for group %d for temperature control method %s must be a "
4273 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4274 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4276 i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
4277 CHECK(nsteps % ir->nstcomm != 0);
4282 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD && ir->comm_mode == ecmNO
4283 && !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
4286 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4287 "rounding errors can lead to build up of kinetic energy of the center of mass");
4290 if (ir->epc == epcPARRINELLORAHMAN && ir->etc == etcNOSEHOOVER)
4293 for (int g = 0; g < ir->opts.ngtc; g++)
4295 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4297 if (ir->tau_p < 1.9 * tau_t_max)
4299 std::string message = gmx::formatString(
4300 "With %s T-coupling and %s p-coupling, "
4301 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4302 etcoupl_names[ir->etc], epcoupl_names[ir->epc], "tau-p", ir->tau_p, "tau-t",
4304 warning(wi, message.c_str());
4308 /* Check for pressure coupling with absolute position restraints */
4309 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4311 absolute_reference(ir, sys, TRUE, AbsRef);
4313 for (m = 0; m < DIM; m++)
4315 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4318 "You are using pressure coupling with absolute position restraints, "
4319 "this will give artifacts. Use the refcoord_scaling option.");
4327 aloopb = gmx_mtop_atomloop_block_init(sys);
4329 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4331 if (atom->q != 0 || atom->qB != 0)
4339 if (EEL_FULL(ir->coulombtype))
4342 "You are using full electrostatics treatment %s for a system without charges.\n"
4343 "This costs a lot of performance for just processing zeros, consider using %s "
4345 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4346 warning(wi, err_buf);
4351 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4354 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4355 "You might want to consider using %s electrostatics.\n",
4357 warning_note(wi, err_buf);
4361 /* Check if combination rules used in LJ-PME are the same as in the force field */
4362 if (EVDW_PME(ir->vdwtype))
4364 check_combination_rules(ir, sys, wi);
4367 /* Generalized reaction field */
4368 if (ir->coulombtype == eelGRF_NOTUSED)
4371 "Generalized reaction-field electrostatics is no longer supported. "
4372 "You can use normal reaction-field instead and compute the reaction-field "
4373 "constant by hand.");
4377 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4379 for (m = 0; (m < DIM); m++)
4381 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4390 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4391 for (const AtomProxy atomP : AtomRange(*sys))
4393 const t_atom& local = atomP.atom();
4394 int i = atomP.globalAtomNumber();
4395 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4398 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4400 for (m = 0; (m < DIM); m++)
4402 acc[m] += ir->opts.acc[i][m] * mgrp[i];
4406 for (m = 0; (m < DIM); m++)
4408 if (fabs(acc[m]) > 1e-6)
4410 const char* dim[DIM] = { "X", "Y", "Z" };
4411 fprintf(stderr, "Net Acceleration in %s direction, will %s be corrected\n", dim[m],
4412 ir->nstcomm != 0 ? "" : "not");
4413 if (ir->nstcomm != 0 && m < ndof_com(ir))
4417 (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4419 ir->opts.acc[i][m] -= acc[m];
4427 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0
4428 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4430 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4438 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4440 if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
4442 absolute_reference(ir, sys, FALSE, AbsRef);
4443 for (m = 0; m < DIM; m++)
4445 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4448 "You are using an absolute reference for pulling, but the rest of "
4449 "the system does not have an absolute reference. This will lead to "
4458 for (i = 0; i < 3; i++)
4460 for (m = 0; m <= i; m++)
4462 if ((ir->epc != epcNO && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4464 for (c = 0; c < ir->pull->ncoord; c++)
4466 if (ir->pull->coord[c].eGeom == epullgDIRPBC && ir->pull->coord[c].vec[m] != 0)
4469 "Can not have dynamic box while using pull geometry '%s' "
4471 EPULLGEOM(ir->pull->coord[c].eGeom), 'x' + m);
4482 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4484 char warn_buf[STRLEN];
4487 ptr = check_box(ir->pbcType, box);
4490 warning_error(wi, ptr);
4493 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4495 if (ir->shake_tol <= 0.0)
4497 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4498 warning_error(wi, warn_buf);
4502 if ((ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4504 /* If we have Lincs constraints: */
4505 if (ir->eI == eiMD && ir->etc == etcNO && ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4508 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4509 warning_note(wi, warn_buf);
4512 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4515 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4517 warning_note(wi, warn_buf);
4519 if (ir->epc == epcMTTK)
4521 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4525 if (bHasAnyConstraints && ir->epc == epcMTTK)
4527 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4530 if (ir->LincsWarnAngle > 90.0)
4532 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4533 warning(wi, warn_buf);
4534 ir->LincsWarnAngle = 90.0;
4537 if (ir->pbcType != PbcType::No)
4539 if (ir->nstlist == 0)
4542 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4543 "atoms might cause the simulation to crash.");
4545 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4548 "ERROR: The cut-off length is longer than half the shortest box vector or "
4549 "longer than the smallest box diagonal element. Increase the box size or "
4550 "decrease rlist.\n");
4551 warning_error(wi, warn_buf);