2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, The GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
51 #include "gromacs/applied_forces/awh/read_params.h"
52 #include "gromacs/fileio/readinp.h"
53 #include "gromacs/fileio/warninp.h"
54 #include "gromacs/gmxlib/network.h"
55 #include "gromacs/gmxpreprocess/toputil.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/calc_verletbuf.h"
60 #include "gromacs/mdrun/mdmodules.h"
61 #include "gromacs/mdtypes/inputrec.h"
62 #include "gromacs/mdtypes/md_enums.h"
63 #include "gromacs/mdtypes/multipletimestepping.h"
64 #include "gromacs/mdtypes/pull_params.h"
65 #include "gromacs/options/options.h"
66 #include "gromacs/options/treesupport.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/selection/indexutil.h"
69 #include "gromacs/topology/block.h"
70 #include "gromacs/topology/ifunc.h"
71 #include "gromacs/topology/index.h"
72 #include "gromacs/topology/mtop_util.h"
73 #include "gromacs/topology/symtab.h"
74 #include "gromacs/topology/topology.h"
75 #include "gromacs/utility/cstringutil.h"
76 #include "gromacs/utility/exceptions.h"
77 #include "gromacs/utility/fatalerror.h"
78 #include "gromacs/utility/filestream.h"
79 #include "gromacs/utility/gmxassert.h"
80 #include "gromacs/utility/ikeyvaluetreeerror.h"
81 #include "gromacs/utility/keyvaluetree.h"
82 #include "gromacs/utility/keyvaluetreebuilder.h"
83 #include "gromacs/utility/keyvaluetreemdpwriter.h"
84 #include "gromacs/utility/keyvaluetreetransform.h"
85 #include "gromacs/utility/mdmodulenotification.h"
86 #include "gromacs/utility/smalloc.h"
87 #include "gromacs/utility/strconvert.h"
88 #include "gromacs/utility/stringcompare.h"
89 #include "gromacs/utility/stringutil.h"
90 #include "gromacs/utility/textwriter.h"
95 /* Resource parameters
96 * Do not change any of these until you read the instruction
97 * in readinp.h. Some cpp's do not take spaces after the backslash
98 * (like the c-shell), which will give you a very weird compiler
102 struct gmx_inputrec_strings
104 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], acc[STRLEN], accgrps[STRLEN], freeze[STRLEN],
105 frdim[STRLEN], energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN],
106 x_compressed_groups[STRLEN], couple_moltype[STRLEN], orirefitgrp[STRLEN],
107 egptable[STRLEN], egpexcl[STRLEN], wall_atomtype[STRLEN], wall_density[STRLEN],
108 deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
109 char fep_lambda[efptNR][STRLEN];
110 char lambda_weights[STRLEN];
111 std::vector<std::string> pullGroupNames;
112 std::vector<std::string> rotateGroupNames;
113 char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
116 static gmx_inputrec_strings* inputrecStrings = nullptr;
118 void init_inputrec_strings()
123 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
124 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
126 inputrecStrings = new gmx_inputrec_strings();
129 void done_inputrec_strings()
131 delete inputrecStrings;
132 inputrecStrings = nullptr;
138 egrptpALL, /* All particles have to be a member of a group. */
139 egrptpALL_GENREST, /* A rest group with name is generated for particles *
140 * that are not part of any group. */
141 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
142 * for the rest group. */
143 egrptpONE /* Merge all selected groups into one group, *
144 * make a rest group for the remaining particles. */
147 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
148 "h-angles", "all-angles", nullptr };
150 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
152 static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lambdas)
157 for (i = 0; i < ntemps; i++)
159 /* simple linear scaling -- allows more control */
160 if (simtemp->eSimTempScale == esimtempLINEAR)
162 simtemp->temperatures[i] =
164 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
166 else if (simtemp->eSimTempScale
167 == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
169 simtemp->temperatures[i] = simtemp->simtemp_low
170 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
171 static_cast<real>((1.0 * i) / (ntemps - 1)));
173 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
175 simtemp->temperatures[i] = simtemp->simtemp_low
176 + (simtemp->simtemp_high - simtemp->simtemp_low)
177 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
182 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
183 gmx_fatal(FARGS, "%s", errorstr);
189 static void _low_check(bool b, const char* s, warninp_t wi)
193 warning_error(wi, s);
197 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
201 if (*p > 0 && *p % nst != 0)
203 /* Round up to the next multiple of nst */
204 *p = ((*p) / nst + 1) * nst;
205 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
210 static int lcd(int n1, int n2)
215 for (i = 2; (i <= n1 && i <= n2); i++)
217 if (n1 % i == 0 && n2 % i == 0)
226 //! Convert legacy mdp entries to modern ones.
227 static void process_interaction_modifier(int* eintmod)
229 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
231 *eintmod = eintmodPOTSHIFT;
235 static void checkMtsRequirement(const t_inputrec& ir, const char* param, const int nstValue, warninp_t wi)
237 GMX_RELEASE_ASSERT(ir.mtsLevels.size() >= 2, "Need at least two levels for MTS");
238 const int mtsFactor = ir.mtsLevels.back().stepFactor;
239 if (nstValue % mtsFactor != 0)
241 auto message = gmx::formatString(
242 "With MTS, %s = %d should be a multiple of mts-factor = %d", param, nstValue, mtsFactor);
243 warning_error(wi, message.c_str());
247 static void setupMtsLevels(gmx::ArrayRef<gmx::MtsLevel> mtsLevels,
248 const t_inputrec& ir,
249 const t_gromppopts& opts,
252 /* MD-VV has no MTS support yet.
253 * SD1 needs different scaling coefficients for the different MTS forces
254 * and the different forces are currently not available in ForceBuffers.
258 auto message = gmx::formatString(
259 "Multiple time stepping is only supported with integrator %s", ei_names[eiMD]);
260 warning_error(wi, message.c_str());
262 if (opts.numMtsLevels != 2)
264 warning_error(wi, "Only mts-levels = 2 is supported");
268 const std::vector<std::string> inputForceGroups = gmx::splitString(opts.mtsLevel2Forces);
269 auto& forceGroups = mtsLevels[1].forceGroups;
270 for (const auto& inputForceGroup : inputForceGroups)
274 for (const auto& forceGroupName : gmx::mtsForceGroupNames)
276 if (gmx::equalCaseInsensitive(inputForceGroup, forceGroupName))
278 forceGroups.set(nameIndex);
286 gmx::formatString("Unknown MTS force group '%s'", inputForceGroup.c_str());
287 warning_error(wi, message.c_str());
291 if (mtsLevels[1].stepFactor <= 1)
293 gmx_fatal(FARGS, "mts-factor should be larger than 1");
296 // Make the level 0 use the complement of the force groups of group 1
297 mtsLevels[0].forceGroups = ~mtsLevels[1].forceGroups;
298 mtsLevels[0].stepFactor = 1;
300 if ((EEL_FULL(ir.coulombtype) || EVDW_PME(ir.vdwtype))
301 && !mtsLevels[1].forceGroups[static_cast<int>(gmx::MtsForceGroups::LongrangeNonbonded)])
304 "With long-range electrostatics and/or LJ treatment, the long-range part "
305 "has to be part of the mts-level2-forces");
308 if (ir.nstcalcenergy > 0)
310 checkMtsRequirement(ir, "nstcalcenergy", ir.nstcalcenergy, wi);
312 checkMtsRequirement(ir, "nstenergy", ir.nstenergy, wi);
313 checkMtsRequirement(ir, "nstlog", ir.nstlog, wi);
314 if (ir.efep != efepNO)
316 checkMtsRequirement(ir, "nstdhdl", ir.fepvals->nstdhdl, wi);
321 const int pullMtsLevel = gmx::forceGroupMtsLevel(ir.mtsLevels, gmx::MtsForceGroups::Pull);
322 if (ir.pull->nstxout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
324 warning_error(wi, "pull-nstxout should be a multiple of mts-factor");
326 if (ir.pull->nstfout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
328 warning_error(wi, "pull-nstfout should be a multiple of mts-factor");
334 void check_ir(const char* mdparin,
335 const gmx::MdModulesNotifier& mdModulesNotifier,
339 /* Check internal consistency.
340 * NOTE: index groups are not set here yet, don't check things
341 * like temperature coupling group options here, but in triple_check
344 /* Strange macro: first one fills the err_buf, and then one can check
345 * the condition, which will print the message and increase the error
348 #define CHECK(b) _low_check(b, err_buf, wi)
349 char err_buf[256], warn_buf[STRLEN];
352 t_lambda* fep = ir->fepvals;
353 t_expanded* expand = ir->expandedvals;
355 set_warning_line(wi, mdparin, -1);
357 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
359 sprintf(warn_buf, "%s electrostatics is no longer supported", eel_names[eelRF_NEC_UNSUPPORTED]);
360 warning_error(wi, warn_buf);
363 /* BASIC CUT-OFF STUFF */
364 if (ir->rcoulomb < 0)
366 warning_error(wi, "rcoulomb should be >= 0");
370 warning_error(wi, "rvdw should be >= 0");
372 if (ir->rlist < 0 && !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
374 warning_error(wi, "rlist should be >= 0");
377 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
378 "neighbour-list update scheme for efficient buffering for improved energy "
379 "conservation, please use the Verlet cut-off scheme instead.)");
380 CHECK(ir->nstlist < 0);
382 process_interaction_modifier(&ir->coulomb_modifier);
383 process_interaction_modifier(&ir->vdw_modifier);
385 if (ir->cutoff_scheme == ecutsGROUP)
388 "The group cutoff scheme has been removed since GROMACS 2020. "
389 "Please use the Verlet cutoff scheme.");
391 if (ir->cutoff_scheme == ecutsVERLET)
395 /* Normal Verlet type neighbor-list, currently only limited feature support */
396 if (inputrec2nboundeddim(ir) < 3)
398 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
401 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
402 if (ir->rcoulomb != ir->rvdw)
404 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
405 // for PME load balancing, we can support this exception.
406 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == evdwCUT
407 && ir->rcoulomb > ir->rvdw);
408 if (!bUsesPmeTwinRangeKernel)
411 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
412 "rcoulomb>rvdw with PME electrostatics)");
416 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
418 if (ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT)
420 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
423 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
425 evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
426 warning_note(wi, warn_buf);
428 ir->vdwtype = evdwCUT;
432 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
433 evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
434 warning_error(wi, warn_buf);
438 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
441 "With Verlet lists only cut-off and PME LJ interactions are supported");
443 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) || EEL_PME(ir->coulombtype)
444 || ir->coulombtype == eelEWALD))
447 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
448 "electrostatics are supported");
450 if (!(ir->coulomb_modifier == eintmodNONE || ir->coulomb_modifier == eintmodPOTSHIFT))
452 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
453 warning_error(wi, warn_buf);
456 if (EEL_USER(ir->coulombtype))
458 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme",
459 eel_names[ir->coulombtype]);
460 warning_error(wi, warn_buf);
463 if (ir->nstlist <= 0)
465 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
468 if (ir->nstlist < 10)
471 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
472 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
476 rc_max = std::max(ir->rvdw, ir->rcoulomb);
480 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
481 ir->verletbuf_tol = 0;
484 else if (ir->verletbuf_tol <= 0)
486 if (ir->verletbuf_tol == 0)
488 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
491 if (ir->rlist < rc_max)
494 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
497 if (ir->rlist == rc_max && ir->nstlist > 1)
501 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
502 "buffer. The cluster pair list does have a buffering effect, but choosing "
503 "a larger rlist might be necessary for good energy conservation.");
508 if (ir->rlist > rc_max)
511 "You have set rlist larger than the interaction cut-off, but you also "
512 "have verlet-buffer-tolerance > 0. Will set rlist using "
513 "verlet-buffer-tolerance.");
516 if (ir->nstlist == 1)
518 /* No buffer required */
523 if (EI_DYNAMICS(ir->eI))
525 if (inputrec2nboundeddim(ir) < 3)
528 "The box volume is required for calculating rlist from the "
529 "energy drift with verlet-buffer-tolerance > 0. You are "
530 "using at least one unbounded dimension, so no volume can be "
531 "computed. Either use a finite box, or set rlist yourself "
532 "together with verlet-buffer-tolerance = -1.");
534 /* Set rlist temporarily so we can continue processing */
539 /* Set the buffer to 5% of the cut-off */
540 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
546 /* GENERAL INTEGRATOR STUFF */
549 if (ir->etc != etcNO)
551 if (EI_RANDOM(ir->eI))
554 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
555 "implicitly. See the documentation for more information on which "
556 "parameters affect temperature for %s.",
557 etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
562 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
564 etcoupl_names[ir->etc], ei_names[ir->eI]);
566 warning_note(wi, warn_buf);
570 if (ir->eI == eiVVAK)
573 "Integrator method %s is implemented primarily for validation purposes; for "
574 "molecular dynamics, you should probably be using %s or %s",
575 ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
576 warning_note(wi, warn_buf);
578 if (!EI_DYNAMICS(ir->eI))
580 if (ir->epc != epcNO)
583 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
584 epcoupl_names[ir->epc], ei_names[ir->eI]);
585 warning_note(wi, warn_buf);
589 if (EI_DYNAMICS(ir->eI))
591 if (ir->nstcalcenergy < 0)
593 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
594 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
596 /* nstcalcenergy larger than nstener does not make sense.
597 * We ideally want nstcalcenergy=nstener.
601 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
605 ir->nstcalcenergy = ir->nstenergy;
609 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
610 || (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
611 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
614 const char* nsten = "nstenergy";
615 const char* nstdh = "nstdhdl";
616 const char* min_name = nsten;
617 int min_nst = ir->nstenergy;
619 /* find the smallest of ( nstenergy, nstdhdl ) */
620 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
621 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
623 min_nst = ir->fepvals->nstdhdl;
626 /* If the user sets nstenergy small, we should respect that */
627 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy,
629 warning_note(wi, warn_buf);
630 ir->nstcalcenergy = min_nst;
633 if (ir->epc != epcNO)
635 if (ir->nstpcouple < 0)
637 ir->nstpcouple = ir_optimal_nstpcouple(ir);
639 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
642 "With multiple time stepping, nstpcouple should be a mutiple of "
647 if (ir->nstcalcenergy > 0)
649 if (ir->efep != efepNO)
651 /* nstdhdl should be a multiple of nstcalcenergy */
652 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
656 /* nstexpanded should be a multiple of nstcalcenergy */
657 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded",
658 &ir->expandedvals->nstexpanded, wi);
660 /* for storing exact averages nstenergy should be
661 * a multiple of nstcalcenergy
663 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
666 // Inquire all MdModules, if their parameters match with the energy
667 // calculation frequency
668 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
669 mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
671 // Emit all errors from the energy calculation frequency checks
672 for (const std::string& energyFrequencyErrorMessage :
673 energyCalculationFrequencyErrors.errorMessages())
675 warning_error(wi, energyFrequencyErrorMessage);
679 if (ir->nsteps == 0 && !ir->bContinuation)
682 "For a correct single-point energy evaluation with nsteps = 0, use "
683 "continuation = yes to avoid constraining the input coordinates.");
687 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
690 "You are doing a continuation with SD or BD, make sure that ld_seed is "
691 "different from the previous run (using ld_seed=-1 will ensure this)");
697 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
698 CHECK(ir->pbcType != PbcType::Xyz);
699 sprintf(err_buf, "with TPI nstlist should be larger than zero");
700 CHECK(ir->nstlist <= 0);
701 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
702 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
706 if ((opts->nshake > 0) && (opts->bMorse))
708 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
709 warning(wi, warn_buf);
712 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
715 "You are doing a continuation with SD or BD, make sure that ld_seed is "
716 "different from the previous run (using ld_seed=-1 will ensure this)");
718 /* verify simulated tempering options */
722 bool bAllTempZero = TRUE;
723 for (i = 0; i < fep->n_lambda; i++)
725 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
726 efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
727 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
728 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
730 bAllTempZero = FALSE;
733 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
734 CHECK(bAllTempZero == TRUE);
736 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
737 CHECK(ir->eI != eiVV);
739 /* check compatability of the temperature coupling with simulated tempering */
741 if (ir->etc == etcNOSEHOOVER)
744 "Nose-Hoover based temperature control such as [%s] my not be "
745 "entirelyconsistent with simulated tempering",
746 etcoupl_names[ir->etc]);
747 warning_note(wi, warn_buf);
750 /* check that the temperatures make sense */
753 "Higher simulated tempering temperature (%g) must be >= than the simulated "
754 "tempering lower temperature (%g)",
755 ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
756 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
758 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero",
759 ir->simtempvals->simtemp_high);
760 CHECK(ir->simtempvals->simtemp_high <= 0);
762 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero",
763 ir->simtempvals->simtemp_low);
764 CHECK(ir->simtempvals->simtemp_low <= 0);
767 /* verify free energy options */
769 if (ir->efep != efepNO)
772 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
773 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
776 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
778 static_cast<int>(fep->sc_r_power));
779 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
782 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
785 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
787 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
789 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
791 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
792 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
794 sprintf(err_buf, "Free-energy not implemented for Ewald");
795 CHECK(ir->coulombtype == eelEWALD);
797 /* check validty of lambda inputs */
798 if (fep->n_lambda == 0)
800 /* Clear output in case of no states:*/
801 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.",
802 fep->init_fep_state);
803 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
807 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d",
808 fep->init_fep_state, fep->n_lambda - 1);
809 CHECK((fep->init_fep_state >= fep->n_lambda));
813 "Lambda state must be set, either with init-lambda-state or with init-lambda");
814 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
817 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
818 "init-lambda-state or with init-lambda, but not both",
819 fep->init_lambda, fep->init_fep_state);
820 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
823 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
827 for (i = 0; i < efptNR; i++)
829 if (fep->separate_dvdl[i])
834 if (n_lambda_terms > 1)
837 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
838 "use init-lambda to set lambda state (except for slow growth). Use "
839 "init-lambda-state instead.");
840 warning(wi, warn_buf);
843 if (n_lambda_terms < 2 && fep->n_lambda > 0)
846 "init-lambda is deprecated for setting lambda state (except for slow "
847 "growth). Use init-lambda-state instead.");
851 for (j = 0; j < efptNR; j++)
853 for (i = 0; i < fep->n_lambda; i++)
855 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
856 efpt_names[j], fep->all_lambda[j][i]);
857 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
861 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
863 for (i = 0; i < fep->n_lambda; i++)
866 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
867 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
868 "crashes, and is not supported.",
869 i, fep->all_lambda[efptVDW][i], fep->all_lambda[efptCOUL][i]);
870 CHECK((fep->sc_alpha > 0)
871 && (((fep->all_lambda[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL][i] < 1.0))
872 && ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0))));
876 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
878 real sigma, lambda, r_sc;
881 /* Maximum estimate for A and B charges equal with lambda power 1 */
883 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
884 1.0 / fep->sc_r_power);
886 "With PME there is a minor soft core effect present at the cut-off, "
887 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
888 "energy conservation, but usually other effects dominate. With a common sigma "
889 "value of %g nm the fraction of the particle-particle potential at the cut-off "
890 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
891 fep->sc_r_power, sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
892 warning_note(wi, warn_buf);
895 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
896 be treated differently, but that's the next step */
898 for (i = 0; i < efptNR; i++)
900 for (j = 0; j < fep->n_lambda; j++)
902 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
903 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
908 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
912 /* checking equilibration of weights inputs for validity */
915 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
917 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
918 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
921 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
923 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
924 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
927 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
928 expand->equil_steps, elmceq_names[elmceqSTEPS]);
929 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
932 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
933 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
934 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
937 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
938 expand->equil_ratio, elmceq_names[elmceqRATIO]);
939 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
942 "weight-equil-number-all-lambda (%d) must be a positive integer if "
943 "lmc-weights-equil=%s",
944 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
945 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
948 "weight-equil-number-samples (%d) must be a positive integer if "
949 "lmc-weights-equil=%s",
950 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
951 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
954 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
955 expand->equil_steps, elmceq_names[elmceqSTEPS]);
956 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
958 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
959 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
960 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
962 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
963 expand->equil_ratio, elmceq_names[elmceqRATIO]);
964 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
966 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
967 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
968 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
970 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
971 CHECK((expand->lmc_repeats <= 0));
972 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
973 CHECK((expand->minvarmin <= 0));
974 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
975 CHECK((expand->c_range < 0));
977 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
979 fep->init_fep_state, expand->lmc_forced_nstart);
980 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
981 && (expand->elmcmove != elmcmoveNO));
982 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
983 CHECK((expand->lmc_forced_nstart < 0));
984 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
985 fep->init_fep_state);
986 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
988 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
989 CHECK((expand->init_wl_delta < 0));
990 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
991 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
992 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
993 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
995 /* if there is no temperature control, we need to specify an MC temperature */
996 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO)
997 && (expand->mc_temp <= 0.0))
1000 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
1001 "to a positive number");
1002 warning_error(wi, err_buf);
1004 if (expand->nstTij > 0)
1006 sprintf(err_buf, "nstlog must be non-zero");
1007 CHECK(ir->nstlog == 0);
1008 // Avoid modulus by zero in the case that already triggered an error exit.
1009 if (ir->nstlog != 0)
1012 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
1013 expand->nstTij, ir->nstlog);
1014 CHECK((expand->nstTij % ir->nstlog) != 0);
1020 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
1021 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
1024 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
1026 if (ir->pbcType == PbcType::No)
1028 if (ir->epc != epcNO)
1030 warning(wi, "Turning off pressure coupling for vacuum system");
1036 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
1037 c_pbcTypeNames[ir->pbcType].c_str());
1038 CHECK(ir->epc != epcNO);
1040 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
1041 CHECK(EEL_FULL(ir->coulombtype));
1043 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
1044 c_pbcTypeNames[ir->pbcType].c_str());
1045 CHECK(ir->eDispCorr != edispcNO);
1048 if (ir->rlist == 0.0)
1051 "can only have neighborlist cut-off zero (=infinite)\n"
1052 "with coulombtype = %s or coulombtype = %s\n"
1053 "without periodic boundary conditions (pbc = %s) and\n"
1054 "rcoulomb and rvdw set to zero",
1055 eel_names[eelCUT], eel_names[eelUSER], c_pbcTypeNames[PbcType::No].c_str());
1056 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
1057 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1059 if (ir->nstlist > 0)
1062 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1063 "nstype=simple and only one MPI rank");
1068 if (ir->nstcomm == 0)
1070 // TODO Change this behaviour. There should be exactly one way
1071 // to turn off an algorithm.
1072 ir->comm_mode = ecmNO;
1074 if (ir->comm_mode != ecmNO)
1076 if (ir->nstcomm < 0)
1078 // TODO Such input was once valid. Now that we've been
1079 // helpful for a few years, we should reject such input,
1080 // lest we have to support every historical decision
1083 "If you want to remove the rotation around the center of mass, you should set "
1084 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1085 "its absolute value");
1086 ir->nstcomm = abs(ir->nstcomm);
1089 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
1092 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
1093 "nstcomm to nstcalcenergy");
1094 ir->nstcomm = ir->nstcalcenergy;
1097 if (ir->comm_mode == ecmANGULAR)
1100 "Can not remove the rotation around the center of mass with periodic "
1102 CHECK(ir->bPeriodicMols);
1103 if (ir->pbcType != PbcType::No)
1106 "Removing the rotation around the center of mass in a periodic system, "
1107 "this can lead to artifacts. Only use this on a single (cluster of) "
1108 "molecules. This cluster should not cross periodic boundaries.");
1113 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
1116 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1117 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1120 warning_note(wi, warn_buf);
1123 /* TEMPERATURE COUPLING */
1124 if (ir->etc == etcYES)
1126 ir->etc = etcBERENDSEN;
1128 "Old option for temperature coupling given: "
1129 "changing \"yes\" to \"Berendsen\"\n");
1132 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
1134 if (ir->opts.nhchainlength < 1)
1137 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1139 ir->opts.nhchainlength);
1140 ir->opts.nhchainlength = 1;
1141 warning(wi, warn_buf);
1144 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1148 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1149 ir->opts.nhchainlength = 1;
1154 ir->opts.nhchainlength = 0;
1157 if (ir->eI == eiVVAK)
1160 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1163 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1166 if (ETC_ANDERSEN(ir->etc))
1168 sprintf(err_buf, "%s temperature control not supported for integrator %s.",
1169 etcoupl_names[ir->etc], ei_names[ir->eI]);
1170 CHECK(!(EI_VV(ir->eI)));
1172 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
1175 "Center of mass removal not necessary for %s. All velocities of coupled "
1176 "groups are rerandomized periodically, so flying ice cube errors will not "
1178 etcoupl_names[ir->etc]);
1179 warning_note(wi, warn_buf);
1183 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1184 "randomized every time step",
1185 ir->nstcomm, etcoupl_names[ir->etc]);
1186 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
1189 if (ir->etc == etcBERENDSEN)
1192 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1193 "might want to consider using the %s thermostat.",
1194 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
1195 warning_note(wi, warn_buf);
1198 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc)) && ir->epc == epcBERENDSEN)
1201 "Using Berendsen pressure coupling invalidates the "
1202 "true ensemble for the thermostat");
1203 warning(wi, warn_buf);
1206 /* PRESSURE COUPLING */
1207 if (ir->epc == epcISOTROPIC)
1209 ir->epc = epcBERENDSEN;
1211 "Old option for pressure coupling given: "
1212 "changing \"Isotropic\" to \"Berendsen\"\n");
1215 if (ir->epc != epcNO)
1217 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1219 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1220 CHECK(ir->tau_p <= 0);
1222 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1225 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1226 "times larger than nstpcouple*dt (%g)",
1227 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1228 warning(wi, warn_buf);
1232 "compressibility must be > 0 when using pressure"
1234 EPCOUPLTYPE(ir->epc));
1235 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1236 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1237 && ir->compress[ZZ][YY] <= 0));
1239 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1242 "You are generating velocities so I am assuming you "
1243 "are equilibrating a system. You are using "
1244 "%s pressure coupling, but this can be "
1245 "unstable for equilibration. If your system crashes, try "
1246 "equilibrating first with Berendsen pressure coupling. If "
1247 "you are not equilibrating the system, you can probably "
1248 "ignore this warning.",
1249 epcoupl_names[ir->epc]);
1250 warning(wi, warn_buf);
1256 if (ir->epc == epcMTTK)
1258 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1262 /* ELECTROSTATICS */
1263 /* More checks are in triple check (grompp.c) */
1265 if (ir->coulombtype == eelSWITCH)
1268 "coulombtype = %s is only for testing purposes and can lead to serious "
1269 "artifacts, advice: use coulombtype = %s",
1270 eel_names[ir->coulombtype], eel_names[eelRF_ZERO]);
1271 warning(wi, warn_buf);
1274 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1277 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1278 "format and exchanging epsilon-r and epsilon-rf",
1280 warning(wi, warn_buf);
1281 ir->epsilon_rf = ir->epsilon_r;
1282 ir->epsilon_r = 1.0;
1285 if (ir->epsilon_r == 0)
1288 "It is pointless to use long-range electrostatics with infinite relative "
1290 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1292 CHECK(EEL_FULL(ir->coulombtype));
1295 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1297 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1298 CHECK(ir->epsilon_r < 0);
1301 if (EEL_RF(ir->coulombtype))
1303 /* reaction field (at the cut-off) */
1305 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1308 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1309 eel_names[ir->coulombtype]);
1310 warning(wi, warn_buf);
1311 ir->epsilon_rf = 0.0;
1314 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1315 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1316 if (ir->epsilon_rf == ir->epsilon_r)
1318 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1319 eel_names[ir->coulombtype]);
1320 warning(wi, warn_buf);
1323 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1324 * means the interaction is zero outside rcoulomb, but it helps to
1325 * provide accurate energy conservation.
1327 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1329 if (ir_coulomb_switched(ir))
1332 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1333 "potential modifier options!",
1334 eel_names[ir->coulombtype]);
1335 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1339 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1342 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1343 "secondary coulomb-modifier.");
1344 CHECK(ir->coulomb_modifier != eintmodNONE);
1346 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1349 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1350 "secondary vdw-modifier.");
1351 CHECK(ir->vdw_modifier != eintmodNONE);
1354 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->vdwtype == evdwSWITCH
1355 || ir->vdwtype == evdwSHIFT)
1358 "The switch/shift interaction settings are just for compatibility; you will get "
1360 "performance from applying potential modifiers to your interactions!\n");
1361 warning_note(wi, warn_buf);
1364 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1366 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1368 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1370 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1371 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1372 "will be good regardless, since ewald_rtol = %g.",
1373 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1374 warning(wi, warn_buf);
1378 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1380 if (ir->rvdw_switch == 0)
1383 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1384 "potential. This suggests it was not set in the mdp, which can lead to large "
1385 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1386 "switching range.");
1387 warning(wi, warn_buf);
1391 if (EEL_FULL(ir->coulombtype))
1393 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER
1394 || ir->coulombtype == eelPMEUSERSWITCH)
1396 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1397 eel_names[ir->coulombtype]);
1398 CHECK(ir->rcoulomb > ir->rlist);
1402 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1404 // TODO: Move these checks into the ewald module with the options class
1406 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1408 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1410 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d",
1411 eel_names[ir->coulombtype], orderMin, orderMax);
1412 warning_error(wi, warn_buf);
1416 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1418 if (ir->ewald_geometry == eewg3D)
1420 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1421 c_pbcTypeNames[ir->pbcType].c_str(), eewg_names[eewg3DC]);
1422 warning(wi, warn_buf);
1424 /* This check avoids extra pbc coding for exclusion corrections */
1425 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1426 CHECK(ir->wall_ewald_zfac < 2);
1428 if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
1430 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1431 eel_names[ir->coulombtype], eewg_names[eewg3DC], c_pbcTypeNames[PbcType::XY].c_str());
1432 warning(wi, warn_buf);
1434 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1436 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1437 CHECK(ir->bPeriodicMols);
1438 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1439 warning_note(wi, warn_buf);
1441 "With epsilon_surface > 0 you can only use domain decomposition "
1442 "when there are only small molecules with all bonds constrained (mdrun will check "
1444 warning_note(wi, warn_buf);
1447 if (ir_vdw_switched(ir))
1449 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1450 CHECK(ir->rvdw_switch >= ir->rvdw);
1452 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1455 "You are applying a switch function to vdw forces or potentials from %g to %g "
1456 "nm, which is more than half the interaction range, whereas switch functions "
1457 "are intended to act only close to the cut-off.",
1458 ir->rvdw_switch, ir->rvdw);
1459 warning_note(wi, warn_buf);
1463 if (ir->vdwtype == evdwPME)
1465 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1467 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1468 evdw_names[ir->vdwtype], eintmod_names[eintmodPOTSHIFT], eintmod_names[eintmodNONE]);
1469 warning_error(wi, err_buf);
1473 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1476 "You have selected user tables with dispersion correction, the dispersion "
1477 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1478 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1479 "really want dispersion correction to -C6/r^6.");
1482 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT) && ir->rvdw != 0)
1484 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1487 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1489 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1492 /* IMPLICIT SOLVENT */
1493 if (ir->coulombtype == eelGB_NOTUSED)
1495 sprintf(warn_buf, "Invalid option %s for coulombtype", eel_names[ir->coulombtype]);
1496 warning_error(wi, warn_buf);
1501 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1506 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1510 /* interpret a number of doubles from a string and put them in an array,
1511 after allocating space for them.
1512 str = the input string
1513 n = the (pre-allocated) number of doubles read
1514 r = the output array of doubles. */
1515 static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
1517 auto values = gmx::splitString(str);
1521 for (int i = 0; i < *n; i++)
1525 (*r)[i] = gmx::fromString<real>(values[i]);
1527 catch (gmx::GromacsException&)
1529 warning_error(wi, "Invalid value " + values[i]
1530 + " in string in mdp file. Expected a real number.");
1536 static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1539 int i, j, max_n_lambda, nweights, nfep[efptNR];
1540 t_lambda* fep = ir->fepvals;
1541 t_expanded* expand = ir->expandedvals;
1542 real** count_fep_lambdas;
1543 bool bOneLambda = TRUE;
1545 snew(count_fep_lambdas, efptNR);
1547 /* FEP input processing */
1548 /* first, identify the number of lambda values for each type.
1549 All that are nonzero must have the same number */
1551 for (i = 0; i < efptNR; i++)
1553 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1556 /* now, determine the number of components. All must be either zero, or equal. */
1559 for (i = 0; i < efptNR; i++)
1561 if (nfep[i] > max_n_lambda)
1563 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1564 must have the same number if its not zero.*/
1569 for (i = 0; i < efptNR; i++)
1573 ir->fepvals->separate_dvdl[i] = FALSE;
1575 else if (nfep[i] == max_n_lambda)
1577 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute
1578 the derivative with respect to the temperature currently */
1580 ir->fepvals->separate_dvdl[i] = TRUE;
1586 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1588 nfep[i], efpt_names[i], max_n_lambda);
1591 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1592 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1594 /* the number of lambdas is the number we've read in, which is either zero
1595 or the same for all */
1596 fep->n_lambda = max_n_lambda;
1598 /* allocate space for the array of lambda values */
1599 snew(fep->all_lambda, efptNR);
1600 /* if init_lambda is defined, we need to set lambda */
1601 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1603 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1605 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1606 for (i = 0; i < efptNR; i++)
1608 snew(fep->all_lambda[i], fep->n_lambda);
1609 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1612 for (j = 0; j < fep->n_lambda; j++)
1614 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1616 sfree(count_fep_lambdas[i]);
1619 sfree(count_fep_lambdas);
1621 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1622 internal bookkeeping -- for now, init_lambda */
1624 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1626 for (i = 0; i < fep->n_lambda; i++)
1628 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1632 /* check to see if only a single component lambda is defined, and soft core is defined.
1633 In this case, turn on coulomb soft core */
1635 if (max_n_lambda == 0)
1641 for (i = 0; i < efptNR; i++)
1643 if ((nfep[i] != 0) && (i != efptFEP))
1649 if ((bOneLambda) && (fep->sc_alpha > 0))
1651 fep->bScCoul = TRUE;
1654 /* Fill in the others with the efptFEP if they are not explicitly
1655 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1656 they are all zero. */
1658 for (i = 0; i < efptNR; i++)
1660 if ((nfep[i] == 0) && (i != efptFEP))
1662 for (j = 0; j < fep->n_lambda; j++)
1664 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1670 /* now read in the weights */
1671 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1674 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1676 else if (nweights != fep->n_lambda)
1678 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1679 nweights, fep->n_lambda);
1681 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1683 expand->nstexpanded = fep->nstdhdl;
1684 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1689 static void do_simtemp_params(t_inputrec* ir)
1692 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1693 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1696 template<typename T>
1697 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1700 for (const auto& input : inputs)
1704 outputs[i] = gmx::fromStdString<T>(input);
1706 catch (gmx::GromacsException&)
1708 auto message = gmx::formatString(
1709 "Invalid value for mdp option %s. %s should only consist of integers separated "
1712 warning_error(wi, message);
1718 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1721 for (const auto& input : inputs)
1725 outputs[i] = gmx::fromString<real>(input);
1727 catch (gmx::GromacsException&)
1729 auto message = gmx::formatString(
1730 "Invalid value for mdp option %s. %s should only consist of real numbers "
1731 "separated by spaces.",
1733 warning_error(wi, message);
1739 static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
1742 for (const auto& input : inputs)
1746 outputs[i][d] = gmx::fromString<real>(input);
1748 catch (gmx::GromacsException&)
1750 auto message = gmx::formatString(
1751 "Invalid value for mdp option %s. %s should only consist of real numbers "
1752 "separated by spaces.",
1754 warning_error(wi, message);
1765 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1767 opts->wall_atomtype[0] = nullptr;
1768 opts->wall_atomtype[1] = nullptr;
1770 ir->wall_atomtype[0] = -1;
1771 ir->wall_atomtype[1] = -1;
1772 ir->wall_density[0] = 0;
1773 ir->wall_density[1] = 0;
1777 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1778 if (wallAtomTypes.size() != size_t(ir->nwall))
1780 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu", ir->nwall,
1781 wallAtomTypes.size());
1783 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1784 for (int i = 0; i < ir->nwall; i++)
1786 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1789 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1791 auto wallDensity = gmx::splitString(wall_density);
1792 if (wallDensity.size() != size_t(ir->nwall))
1794 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall,
1795 wallDensity.size());
1797 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1798 for (int i = 0; i < ir->nwall; i++)
1800 if (ir->wall_density[i] <= 0)
1802 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1809 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1813 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1814 for (int i = 0; i < nwall; i++)
1816 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1817 grps->emplace_back(groups->groupNames.size() - 1);
1822 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1824 /* read expanded ensemble parameters */
1825 printStringNewline(inp, "expanded ensemble variables");
1826 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1827 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1828 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1829 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1830 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1831 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1832 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1833 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1834 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1835 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1836 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1837 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1838 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1839 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1840 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1841 expand->bSymmetrizedTMatrix =
1842 (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1843 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1844 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1845 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1846 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1847 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1848 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1849 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1852 /*! \brief Return whether an end state with the given coupling-lambda
1853 * value describes fully-interacting VDW.
1855 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1856 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1858 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1860 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1866 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1869 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1871 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1873 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1876 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1877 std::string message = gmx::formatExceptionMessageToString(*ex);
1878 warning_error(wi_, message.c_str());
1883 std::string getOptionName(const gmx::KeyValueTreePath& context)
1885 if (mapping_ != nullptr)
1887 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1888 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1891 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1896 const gmx::IKeyValueTreeBackMapping* mapping_;
1901 void get_ir(const char* mdparin,
1902 const char* mdparout,
1903 gmx::MDModules* mdModules,
1906 WriteMdpHeader writeMdpHeader,
1910 double dumdub[2][6];
1912 char warn_buf[STRLEN];
1913 t_lambda* fep = ir->fepvals;
1914 t_expanded* expand = ir->expandedvals;
1916 const char* no_names[] = { "no", nullptr };
1918 init_inputrec_strings();
1919 gmx::TextInputFile stream(mdparin);
1920 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1922 snew(dumstr[0], STRLEN);
1923 snew(dumstr[1], STRLEN);
1925 /* ignore the following deprecated commands */
1926 replace_inp_entry(inp, "title", nullptr);
1927 replace_inp_entry(inp, "cpp", nullptr);
1928 replace_inp_entry(inp, "domain-decomposition", nullptr);
1929 replace_inp_entry(inp, "andersen-seed", nullptr);
1930 replace_inp_entry(inp, "dihre", nullptr);
1931 replace_inp_entry(inp, "dihre-fc", nullptr);
1932 replace_inp_entry(inp, "dihre-tau", nullptr);
1933 replace_inp_entry(inp, "nstdihreout", nullptr);
1934 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1935 replace_inp_entry(inp, "optimize-fft", nullptr);
1936 replace_inp_entry(inp, "adress_type", nullptr);
1937 replace_inp_entry(inp, "adress_const_wf", nullptr);
1938 replace_inp_entry(inp, "adress_ex_width", nullptr);
1939 replace_inp_entry(inp, "adress_hy_width", nullptr);
1940 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1941 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1942 replace_inp_entry(inp, "adress_site", nullptr);
1943 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1944 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1945 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1946 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1947 replace_inp_entry(inp, "rlistlong", nullptr);
1948 replace_inp_entry(inp, "nstcalclr", nullptr);
1949 replace_inp_entry(inp, "pull-print-com2", nullptr);
1950 replace_inp_entry(inp, "gb-algorithm", nullptr);
1951 replace_inp_entry(inp, "nstgbradii", nullptr);
1952 replace_inp_entry(inp, "rgbradii", nullptr);
1953 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1954 replace_inp_entry(inp, "gb-saltconc", nullptr);
1955 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1956 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1957 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1958 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1959 replace_inp_entry(inp, "sa-algorithm", nullptr);
1960 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1961 replace_inp_entry(inp, "ns-type", nullptr);
1963 /* replace the following commands with the clearer new versions*/
1964 replace_inp_entry(inp, "unconstrained-start", "continuation");
1965 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1966 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1967 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1968 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1969 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1970 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1972 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1973 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1974 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1975 setStringEntry(&inp, "include", opts->include, nullptr);
1976 printStringNoNewline(
1977 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1978 setStringEntry(&inp, "define", opts->define, nullptr);
1980 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1981 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1982 printStringNoNewline(&inp, "Start time and timestep in ps");
1983 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1984 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1985 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1986 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1987 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1988 printStringNoNewline(
1989 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1990 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1991 printStringNoNewline(&inp, "Multiple time-stepping");
1992 ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
1995 opts->numMtsLevels = get_eint(&inp, "mts-levels", 2, wi);
1996 ir->mtsLevels.resize(2);
1997 gmx::MtsLevel& mtsLevel = ir->mtsLevels[1];
1998 opts->mtsLevel2Forces = setStringEntry(&inp, "mts-level2-forces",
1999 "longrange-nonbonded nonbonded pair dihedral");
2000 mtsLevel.stepFactor = get_eint(&inp, "mts-level2-factor", 2, wi);
2002 // We clear after reading without dynamics to not force the user to remove MTS mdp options
2003 if (!EI_DYNAMICS(ir->eI))
2006 ir->mtsLevels.clear();
2009 printStringNoNewline(&inp, "mode for center of mass motion removal");
2010 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
2011 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
2012 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
2013 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
2014 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
2016 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2017 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2018 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2019 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2022 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2023 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2024 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2025 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2026 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2027 ir->niter = get_eint(&inp, "niter", 20, wi);
2028 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2029 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2030 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2031 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2032 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2034 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2035 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2037 /* Output options */
2038 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2039 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2040 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2041 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2042 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2043 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2044 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2045 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2046 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2047 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2048 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2049 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2050 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2051 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2052 printStringNoNewline(&inp, "default, all atoms will be written.");
2053 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2054 printStringNoNewline(&inp, "Selection of energy groups");
2055 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2057 /* Neighbor searching */
2058 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2059 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2060 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
2061 printStringNoNewline(&inp, "nblist update frequency");
2062 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2063 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2064 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2065 std::vector<const char*> pbcTypesNamesChar;
2066 for (const auto& pbcTypeName : c_pbcTypeNames)
2068 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2070 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2071 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
2072 printStringNoNewline(&inp,
2073 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2074 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2075 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2076 printStringNoNewline(&inp, "nblist cut-off");
2077 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2078 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2080 /* Electrostatics */
2081 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2082 printStringNoNewline(&inp, "Method for doing electrostatics");
2083 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
2084 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
2085 printStringNoNewline(&inp, "cut-off lengths");
2086 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2087 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2088 printStringNoNewline(&inp,
2089 "Relative dielectric constant for the medium and the reaction field");
2090 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2091 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2092 printStringNoNewline(&inp, "Method for doing Van der Waals");
2093 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
2094 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
2095 printStringNoNewline(&inp, "cut-off lengths");
2096 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2097 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2098 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2099 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
2100 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2101 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2102 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2103 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2104 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2105 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2106 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2107 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2108 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2109 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2110 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2111 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2112 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2113 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2114 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
2115 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
2116 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2118 /* Implicit solvation is no longer supported, but we need grompp
2119 to be able to refuse old .mdp files that would have built a tpr
2120 to run it. Thus, only "no" is accepted. */
2121 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2123 /* Coupling stuff */
2124 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2125 printStringNoNewline(&inp, "Temperature coupling");
2126 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
2127 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2128 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2129 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
2130 printStringNoNewline(&inp, "Groups to couple separately");
2131 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2132 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2133 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2134 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2135 printStringNoNewline(&inp, "pressure coupling");
2136 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
2137 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
2138 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2139 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2140 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2141 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2142 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2143 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2144 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
2147 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2148 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
2149 printStringNoNewline(&inp, "Groups treated with MiMiC");
2150 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2152 /* Simulated annealing */
2153 printStringNewline(&inp, "SIMULATED ANNEALING");
2154 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2155 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2156 printStringNoNewline(&inp,
2157 "Number of time points to use for specifying annealing in each group");
2158 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2159 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2160 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2161 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2162 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2165 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2166 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
2167 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2168 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2171 printStringNewline(&inp, "OPTIONS FOR BONDS");
2172 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2173 printStringNoNewline(&inp, "Type of constraint algorithm");
2174 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
2175 printStringNoNewline(&inp, "Do not constrain the start configuration");
2176 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
2177 printStringNoNewline(&inp,
2178 "Use successive overrelaxation to reduce the number of shake iterations");
2179 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
2180 printStringNoNewline(&inp, "Relative tolerance of shake");
2181 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2182 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2183 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2184 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2185 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2186 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2187 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2188 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2189 printStringNoNewline(&inp, "rotates over more degrees than");
2190 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2191 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2192 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
2194 /* Energy group exclusions */
2195 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2196 printStringNoNewline(
2197 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2198 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2201 printStringNewline(&inp, "WALLS");
2202 printStringNoNewline(
2203 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2204 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2205 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
2206 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2207 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2208 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2209 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2212 printStringNewline(&inp, "COM PULLING");
2213 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2216 ir->pull = std::make_unique<pull_params_t>();
2217 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
2221 for (int c = 0; c < ir->pull->ncoord; c++)
2223 if (ir->pull->coord[c].eType == epullCONSTRAINT)
2226 "Constraint COM pulling is not supported in combination with "
2227 "multiple time stepping");
2235 NOTE: needs COM pulling or free energy input */
2236 printStringNewline(&inp, "AWH biasing");
2237 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2240 ir->awhParams = gmx::readAwhParams(&inp, wi);
2243 /* Enforced rotation */
2244 printStringNewline(&inp, "ENFORCED ROTATION");
2245 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2246 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2250 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2253 /* Interactive MD */
2255 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2256 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2257 if (inputrecStrings->imd_grp[0] != '\0')
2264 printStringNewline(&inp, "NMR refinement stuff");
2265 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2266 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2267 printStringNoNewline(
2268 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2269 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2270 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2271 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2272 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2273 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2274 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2275 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2276 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2277 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2278 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2279 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2280 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2281 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2282 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2283 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2285 /* free energy variables */
2286 printStringNewline(&inp, "Free energy variables");
2287 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2288 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2289 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2290 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2291 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2293 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2295 it was not entered */
2296 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2297 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2298 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2299 setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
2300 setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
2301 setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
2302 setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
2303 setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
2304 setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
2305 setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
2306 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2307 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2308 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2309 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2310 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2311 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2312 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2313 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
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);
2316 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2317 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2318 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2319 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2321 /* Non-equilibrium MD stuff */
2322 printStringNewline(&inp, "Non-equilibrium MD stuff");
2323 setStringEntry(&inp, "acc-grps", inputrecStrings->accgrps, nullptr);
2324 setStringEntry(&inp, "accelerate", inputrecStrings->acc, nullptr);
2325 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2326 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2327 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2328 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2330 /* simulated tempering variables */
2331 printStringNewline(&inp, "simulated tempering variables");
2332 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2333 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2334 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2335 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2337 /* expanded ensemble variables */
2338 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2340 read_expandedparams(&inp, expand, wi);
2343 /* Electric fields */
2345 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2346 gmx::KeyValueTreeTransformer transform;
2347 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2348 mdModules->initMdpTransform(transform.rules());
2349 for (const auto& path : transform.mappedPaths())
2351 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2352 mark_einp_set(inp, path[0].c_str());
2354 MdpErrorHandler errorHandler(wi);
2355 auto result = transform.transform(convertedValues, &errorHandler);
2356 ir->params = new gmx::KeyValueTreeObject(result.object());
2357 mdModules->adjustInputrecBasedOnModules(ir);
2358 errorHandler.setBackMapping(result.backMapping());
2359 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2362 /* Ion/water position swapping ("computational electrophysiology") */
2363 printStringNewline(&inp,
2364 "Ion/water position swapping for computational electrophysiology setups");
2365 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2366 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2367 if (ir->eSwapCoords != eswapNO)
2374 printStringNoNewline(&inp, "Swap attempt frequency");
2375 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2376 printStringNoNewline(&inp, "Number of ion types to be controlled");
2377 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2380 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2382 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2383 snew(ir->swap->grp, ir->swap->ngrp);
2384 for (i = 0; i < ir->swap->ngrp; i++)
2386 snew(ir->swap->grp[i].molname, STRLEN);
2388 printStringNoNewline(&inp,
2389 "Two index groups that contain the compartment-partitioning atoms");
2390 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2391 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2392 printStringNoNewline(&inp,
2393 "Use center of mass of split groups (yes/no), otherwise center of "
2394 "geometry is used");
2395 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2396 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2398 printStringNoNewline(&inp, "Name of solvent molecules");
2399 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2401 printStringNoNewline(&inp,
2402 "Split cylinder: radius, upper and lower extension (nm) (this will "
2403 "define the channels)");
2404 printStringNoNewline(&inp,
2405 "Note that the split cylinder settings do not have an influence on "
2406 "the swapping protocol,");
2407 printStringNoNewline(
2409 "however, if correctly defined, the permeation events are recorded per channel");
2410 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2411 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2412 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2413 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2414 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2415 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2417 printStringNoNewline(
2419 "Average the number of ions per compartment over these many swap attempt steps");
2420 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2422 printStringNoNewline(
2423 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2424 printStringNoNewline(
2425 &inp, "and the requested number of ions of this type in compartments A and B");
2426 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2427 for (i = 0; i < nIonTypes; i++)
2429 int ig = eSwapFixedGrpNR + i;
2431 sprintf(buf, "iontype%d-name", i);
2432 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2433 sprintf(buf, "iontype%d-in-A", i);
2434 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2435 sprintf(buf, "iontype%d-in-B", i);
2436 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2439 printStringNoNewline(
2441 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2442 printStringNoNewline(
2444 "at maximum distance (= bulk concentration) to the split group layers. However,");
2445 printStringNoNewline(&inp,
2446 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2447 "layer from the middle at 0.0");
2448 printStringNoNewline(&inp,
2449 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2450 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2451 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2452 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2453 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2455 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2458 printStringNoNewline(
2459 &inp, "Start to swap ions if threshold difference to requested count is reached");
2460 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2463 /* AdResS is no longer supported, but we need grompp to be able to
2464 refuse to process old .mdp files that used it. */
2465 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2467 /* User defined thingies */
2468 printStringNewline(&inp, "User defined thingies");
2469 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2470 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2471 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2472 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2473 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2474 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2475 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2476 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2477 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2478 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2482 gmx::TextOutputFile stream(mdparout);
2483 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2485 // Transform module data into a flat key-value tree for output.
2486 gmx::KeyValueTreeBuilder builder;
2487 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2488 mdModules->buildMdpOutput(&builderObject);
2490 gmx::TextWriter writer(&stream);
2491 writeKeyValueTreeAsMdp(&writer, builder.build());
2496 /* Process options if necessary */
2497 for (m = 0; m < 2; m++)
2499 for (i = 0; i < 2 * DIM; i++)
2508 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2512 "Pressure coupling incorrect number of values (I need exactly 1)");
2514 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2516 case epctSEMIISOTROPIC:
2517 case epctSURFACETENSION:
2518 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2522 "Pressure coupling incorrect number of values (I need exactly 2)");
2524 dumdub[m][YY] = dumdub[m][XX];
2526 case epctANISOTROPIC:
2527 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf", &(dumdub[m][XX]), &(dumdub[m][YY]),
2528 &(dumdub[m][ZZ]), &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5]))
2533 "Pressure coupling incorrect number of values (I need exactly 6)");
2537 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2538 epcoupltype_names[ir->epct]);
2542 clear_mat(ir->ref_p);
2543 clear_mat(ir->compress);
2544 for (i = 0; i < DIM; i++)
2546 ir->ref_p[i][i] = dumdub[1][i];
2547 ir->compress[i][i] = dumdub[0][i];
2549 if (ir->epct == epctANISOTROPIC)
2551 ir->ref_p[XX][YY] = dumdub[1][3];
2552 ir->ref_p[XX][ZZ] = dumdub[1][4];
2553 ir->ref_p[YY][ZZ] = dumdub[1][5];
2554 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2557 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2558 "apply a threefold shear stress?\n");
2560 ir->compress[XX][YY] = dumdub[0][3];
2561 ir->compress[XX][ZZ] = dumdub[0][4];
2562 ir->compress[YY][ZZ] = dumdub[0][5];
2563 for (i = 0; i < DIM; i++)
2565 for (m = 0; m < i; m++)
2567 ir->ref_p[i][m] = ir->ref_p[m][i];
2568 ir->compress[i][m] = ir->compress[m][i];
2573 if (ir->comm_mode == ecmNO)
2578 opts->couple_moltype = nullptr;
2579 if (strlen(inputrecStrings->couple_moltype) > 0)
2581 if (ir->efep != efepNO)
2583 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2584 if (opts->couple_lam0 == opts->couple_lam1)
2586 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2588 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2592 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2599 "Free energy is turned off, so we will not decouple the molecule listed "
2603 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2604 if (ir->efep != efepNO)
2606 if (fep->delta_lambda != 0)
2608 ir->efep = efepSLOWGROWTH;
2612 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2614 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2616 "Old option for dhdl-print-energy given: "
2617 "changing \"yes\" to \"total\"\n");
2620 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2622 /* always print out the energy to dhdl if we are doing
2623 expanded ensemble, since we need the total energy for
2624 analysis if the temperature is changing. In some
2625 conditions one may only want the potential energy, so
2626 we will allow that if the appropriate mdp setting has
2627 been enabled. Otherwise, total it is:
2629 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2632 if ((ir->efep != efepNO) || ir->bSimTemp)
2634 ir->bExpanded = FALSE;
2635 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2637 ir->bExpanded = TRUE;
2639 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2640 if (ir->bSimTemp) /* done after fep params */
2642 do_simtemp_params(ir);
2645 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2646 * setup and not on the old way of specifying the free-energy setup,
2647 * we should check for using soft-core when not needed, since that
2648 * can complicate the sampling significantly.
2649 * Note that we only check for the automated coupling setup.
2650 * If the (advanced) user does FEP through manual topology changes,
2651 * this check will not be triggered.
2653 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 && ir->fepvals->sc_alpha != 0
2654 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2657 "You are using soft-core interactions while the Van der Waals interactions are "
2658 "not decoupled (note that the sc-coul option is only active when using lambda "
2659 "states). Although this will not lead to errors, you will need much more "
2660 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2665 ir->fepvals->n_lambda = 0;
2668 /* WALL PARAMETERS */
2670 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2672 /* ORIENTATION RESTRAINT PARAMETERS */
2674 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2676 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2679 /* DEFORMATION PARAMETERS */
2681 clear_mat(ir->deform);
2682 for (i = 0; i < 6; i++)
2687 double gmx_unused canary;
2688 int ndeform = sscanf(inputrecStrings->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]),
2689 &(dumdub[0][1]), &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]),
2690 &(dumdub[0][5]), &canary);
2692 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2696 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2697 inputrecStrings->deform)
2700 for (i = 0; i < 3; i++)
2702 ir->deform[i][i] = dumdub[0][i];
2704 ir->deform[YY][XX] = dumdub[0][3];
2705 ir->deform[ZZ][XX] = dumdub[0][4];
2706 ir->deform[ZZ][YY] = dumdub[0][5];
2707 if (ir->epc != epcNO)
2709 for (i = 0; i < 3; i++)
2711 for (j = 0; j <= i; j++)
2713 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2715 warning_error(wi, "A box element has deform set and compressibility > 0");
2719 for (i = 0; i < 3; i++)
2721 for (j = 0; j < i; j++)
2723 if (ir->deform[i][j] != 0)
2725 for (m = j; m < DIM; m++)
2727 if (ir->compress[m][j] != 0)
2730 "An off-diagonal box element has deform set while "
2731 "compressibility > 0 for the same component of another box "
2732 "vector, this might lead to spurious periodicity effects.");
2733 warning(wi, warn_buf);
2741 /* Ion/water position swapping checks */
2742 if (ir->eSwapCoords != eswapNO)
2744 if (ir->swap->nstswap < 1)
2746 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2748 if (ir->swap->nAverage < 1)
2750 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2752 if (ir->swap->threshold < 1.0)
2754 warning_error(wi, "Ion count threshold must be at least 1.\n");
2758 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2761 setupMtsLevels(ir->mtsLevels, *ir, *opts, wi);
2766 gmx::checkAwhParams(ir->awhParams, ir, wi);
2773 /* We would like gn to be const as well, but C doesn't allow this */
2774 /* TODO this is utility functionality (search for the index of a
2775 string in a collection), so should be refactored and located more
2777 int search_string(const char* s, int ng, char* gn[])
2781 for (i = 0; (i < ng); i++)
2783 if (gmx_strcasecmp(s, gn[i]) == 0)
2790 "Group %s referenced in the .mdp file was not found in the index file.\n"
2791 "Group names must match either [moleculetype] names or custom index group\n"
2792 "names, in which case you must supply an index file to the '-n' option\n"
2797 static void do_numbering(int natoms,
2798 SimulationGroups* groups,
2799 gmx::ArrayRef<std::string> groupsFromMdpFile,
2802 SimulationAtomGroupType gtype,
2808 unsigned short* cbuf;
2809 AtomGroupIndices* grps = &(groups->groups[gtype]);
2810 int j, gid, aj, ognr, ntot = 0;
2812 char warn_buf[STRLEN];
2814 title = shortName(gtype);
2817 /* Mark all id's as not set */
2818 for (int i = 0; (i < natoms); i++)
2823 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2825 /* Lookup the group name in the block structure */
2826 gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2827 if ((grptp != egrptpONE) || (i == 0))
2829 grps->emplace_back(gid);
2832 /* Now go over the atoms in the group */
2833 for (j = block->index[gid]; (j < block->index[gid + 1]); j++)
2838 /* Range checking */
2839 if ((aj < 0) || (aj >= natoms))
2841 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2843 /* Lookup up the old group number */
2847 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
2852 /* Store the group number in buffer */
2853 if (grptp == egrptpONE)
2866 /* Now check whether we have done all atoms */
2869 if (grptp == egrptpALL)
2871 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2873 else if (grptp == egrptpPART)
2875 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2876 warning_note(wi, warn_buf);
2878 /* Assign all atoms currently unassigned to a rest group */
2879 for (j = 0; (j < natoms); j++)
2881 if (cbuf[j] == NOGID)
2883 cbuf[j] = grps->size();
2886 if (grptp != egrptpPART)
2890 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title,
2893 /* Add group name "rest" */
2894 grps->emplace_back(restnm);
2896 /* Assign the rest name to all atoms not currently assigned to a group */
2897 for (j = 0; (j < natoms); j++)
2899 if (cbuf[j] == NOGID)
2901 // group size was not updated before this here, so need to use -1.
2902 cbuf[j] = grps->size() - 1;
2908 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2910 /* All atoms are part of one (or no) group, no index required */
2911 groups->groupNumbers[gtype].clear();
2915 for (int j = 0; (j < natoms); j++)
2917 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2924 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2927 pull_params_t* pull;
2928 int natoms, imin, jmin;
2929 int * nrdf2, *na_vcm, na_tot;
2930 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2935 * First calc 3xnr-atoms for each group
2936 * then subtract half a degree of freedom for each constraint
2938 * Only atoms and nuclei contribute to the degrees of freedom...
2943 const SimulationGroups& groups = mtop->groups;
2944 natoms = mtop->natoms;
2946 /* Allocate one more for a possible rest group */
2947 /* We need to sum degrees of freedom into doubles,
2948 * since floats give too low nrdf's above 3 million atoms.
2950 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
2951 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2952 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2953 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2954 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2956 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2960 for (gmx::index i = 0;
2961 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; i++)
2964 clear_ivec(dof_vcm[i]);
2966 nrdf_vcm_sub[i] = 0;
2968 snew(nrdf2, natoms);
2969 for (const AtomProxy atomP : AtomRange(*mtop))
2971 const t_atom& local = atomP.atom();
2972 int i = atomP.globalAtomNumber();
2974 if (local.ptype == eptAtom || local.ptype == eptNucleus)
2976 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
2977 for (int d = 0; d < DIM; d++)
2979 if (opts->nFreeze[g][d] == 0)
2981 /* Add one DOF for particle i (counted as 2*1) */
2983 /* VCM group i has dim d as a DOF */
2984 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
2988 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
2990 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
2996 for (const gmx_molblock_t& molb : mtop->molblock)
2998 const gmx_moltype_t& molt = mtop->moltype[molb.type];
2999 const t_atom* atom = molt.atoms.atom;
3000 for (int mol = 0; mol < molb.nmol; mol++)
3002 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
3004 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
3005 for (int i = 0; i < molt.ilist[ftype].size();)
3007 /* Subtract degrees of freedom for the constraints,
3008 * if the particles still have degrees of freedom left.
3009 * If one of the particles is a vsite or a shell, then all
3010 * constraint motion will go there, but since they do not
3011 * contribute to the constraints the degrees of freedom do not
3014 int ai = as + ia[i + 1];
3015 int aj = as + ia[i + 2];
3016 if (((atom[ia[i + 1]].ptype == eptNucleus) || (atom[ia[i + 1]].ptype == eptAtom))
3017 && ((atom[ia[i + 2]].ptype == eptNucleus) || (atom[ia[i + 2]].ptype == eptAtom)))
3035 imin = std::min(imin, nrdf2[ai]);
3036 jmin = std::min(jmin, nrdf2[aj]);
3039 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3041 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3043 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3045 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3048 i += interaction_function[ftype].nratoms + 1;
3051 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3052 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3054 /* Subtract 1 dof from every atom in the SETTLE */
3055 for (int j = 0; j < 3; j++)
3057 int ai = as + ia[i + 1 + j];
3058 imin = std::min(2, nrdf2[ai]);
3060 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3062 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3067 as += molt.atoms.nr;
3073 /* Correct nrdf for the COM constraints.
3074 * We correct using the TC and VCM group of the first atom
3075 * in the reference and pull group. If atoms in one pull group
3076 * belong to different TC or VCM groups it is anyhow difficult
3077 * to determine the optimal nrdf assignment.
3079 pull = ir->pull.get();
3081 for (int i = 0; i < pull->ncoord; i++)
3083 if (pull->coord[i].eType != epullCONSTRAINT)
3090 for (int j = 0; j < 2; j++)
3092 const t_pull_group* pgrp;
3094 pgrp = &pull->group[pull->coord[i].group[j]];
3096 if (!pgrp->ind.empty())
3098 /* Subtract 1/2 dof from each group */
3099 int ai = pgrp->ind[0];
3100 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3102 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3104 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3107 "Center of mass pulling constraints caused the number of degrees "
3108 "of freedom for temperature coupling group %s to be negative",
3109 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3110 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3115 /* We need to subtract the whole DOF from group j=1 */
3122 if (ir->nstcomm != 0)
3124 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3125 "Expect at least one group when removing COM motion");
3127 /* We remove COM motion up to dim ndof_com() */
3128 const int ndim_rm_vcm = ndof_com(ir);
3130 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3131 * the number of degrees of freedom in each vcm group when COM
3132 * translation is removed and 6 when rotation is removed as well.
3133 * Note that we do not and should not include the rest group here.
3135 for (gmx::index j = 0;
3136 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]); j++)
3138 switch (ir->comm_mode)
3141 case ecmLINEAR_ACCELERATION_CORRECTION:
3142 nrdf_vcm_sub[j] = 0;
3143 for (int d = 0; d < ndim_rm_vcm; d++)
3151 case ecmANGULAR: nrdf_vcm_sub[j] = 6; break;
3152 default: gmx_incons("Checking comm_mode");
3156 for (gmx::index i = 0;
3157 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
3159 /* Count the number of atoms of TC group i for every VCM group */
3160 for (gmx::index j = 0;
3161 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3166 for (int ai = 0; ai < natoms; ai++)
3168 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3170 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3174 /* Correct for VCM removal according to the fraction of each VCM
3175 * group present in this TC group.
3177 nrdf_uc = nrdf_tc[i];
3179 for (gmx::index j = 0;
3180 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3182 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3184 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3185 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3190 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3192 opts->nrdf[i] = nrdf_tc[i];
3193 if (opts->nrdf[i] < 0)
3197 fprintf(stderr, "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3198 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
3206 sfree(nrdf_vcm_sub);
3209 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3211 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3212 * But since this is much larger than STRLEN, such a line can not be parsed.
3213 * The real maximum is the number of names that fit in a string: STRLEN/2.
3215 #define EGP_MAX (STRLEN / 2)
3219 auto names = gmx::splitString(val);
3220 if (names.size() % 2 != 0)
3222 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3224 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3226 for (size_t i = 0; i < names.size() / 2; i++)
3228 // TODO this needs to be replaced by a solution using std::find_if
3232 names[2 * i].c_str(),
3233 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3239 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3244 names[2 * i + 1].c_str(),
3245 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3251 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3253 if ((j < nr) && (k < nr))
3255 ir->opts.egp_flags[nr * j + k] |= flag;
3256 ir->opts.egp_flags[nr * k + j] |= flag;
3265 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3267 int ig = -1, i = 0, gind;
3271 /* Just a quick check here, more thorough checks are in mdrun */
3272 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3274 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3277 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3278 for (ig = 0; ig < swap->ngrp; ig++)
3280 swapg = &swap->grp[ig];
3281 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3282 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3286 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3287 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap", swap->grp[ig].molname, swapg->nat);
3288 snew(swapg->ind, swapg->nat);
3289 for (i = 0; i < swapg->nat; i++)
3291 swapg->ind[i] = grps->a[grps->index[gind] + i];
3296 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3302 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3307 ig = search_string(IMDgname, grps->nr, gnames);
3308 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3310 if (IMDgroup->nat > 0)
3313 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3315 IMDgname, IMDgroup->nat);
3316 snew(IMDgroup->ind, IMDgroup->nat);
3317 for (i = 0; i < IMDgroup->nat; i++)
3319 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3324 /* Checks whether atoms are both part of a COM removal group and frozen.
3325 * If a fully frozen atom is part of a COM removal group, it is removed
3326 * from the COM removal group. A note is issued if such atoms are present.
3327 * A warning is issued for atom with one or two dimensions frozen that
3328 * are part of a COM removal group (mdrun would need to compute COM mass
3329 * per dimension to handle this correctly).
3330 * Also issues a warning when non-frozen atoms are not part of a COM
3331 * removal group while COM removal is active.
3333 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3335 const t_grpopts& opts,
3338 const int vcmRestGroup =
3339 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3341 int numFullyFrozenVcmAtoms = 0;
3342 int numPartiallyFrozenVcmAtoms = 0;
3343 int numNonVcmAtoms = 0;
3344 for (int a = 0; a < numAtoms; a++)
3346 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3347 int numFrozenDims = 0;
3348 for (int d = 0; d < DIM; d++)
3350 numFrozenDims += opts.nFreeze[freezeGroup][d];
3353 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3354 if (vcmGroup < vcmRestGroup)
3356 if (numFrozenDims == DIM)
3358 /* Do not remove COM motion for this fully frozen atom */
3359 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3361 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3364 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3365 numFullyFrozenVcmAtoms++;
3367 else if (numFrozenDims > 0)
3369 numPartiallyFrozenVcmAtoms++;
3372 else if (numFrozenDims < DIM)
3378 if (numFullyFrozenVcmAtoms > 0)
3380 std::string warningText = gmx::formatString(
3381 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3382 "removing these atoms from the COMM removal group(s)",
3383 numFullyFrozenVcmAtoms);
3384 warning_note(wi, warningText.c_str());
3386 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3388 std::string warningText = gmx::formatString(
3389 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3390 "removal group(s), due to limitations in the code these still contribute to the "
3391 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3393 numPartiallyFrozenVcmAtoms, DIM);
3394 warning(wi, warningText.c_str());
3396 if (numNonVcmAtoms > 0)
3398 std::string warningText = gmx::formatString(
3399 "%d atoms are not part of any center of mass motion removal group.\n"
3400 "This may lead to artifacts.\n"
3401 "In most cases one should use one group for the whole system.",
3403 warning(wi, warningText.c_str());
3407 void do_index(const char* mdparin,
3411 const gmx::MdModulesNotifier& notifier,
3415 t_blocka* defaultIndexGroups;
3423 int i, j, k, restnm;
3424 bool bExcl, bTable, bAnneal;
3425 char warn_buf[STRLEN];
3429 fprintf(stderr, "processing index file...\n");
3433 snew(defaultIndexGroups, 1);
3434 snew(defaultIndexGroups->index, 1);
3436 atoms_all = gmx_mtop_global_atoms(mtop);
3437 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3438 done_atom(&atoms_all);
3442 defaultIndexGroups = init_index(ndx, &gnames);
3445 SimulationGroups* groups = &mtop->groups;
3446 natoms = mtop->natoms;
3447 symtab = &mtop->symtab;
3449 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3451 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3453 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3454 restnm = groups->groupNames.size() - 1;
3455 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3456 srenew(gnames, defaultIndexGroups->nr + 1);
3457 gnames[restnm] = *(groups->groupNames.back());
3459 set_warning_line(wi, mdparin, -1);
3461 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3462 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3463 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3464 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3465 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3468 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3470 temperatureCouplingGroupNames.size(), temperatureCouplingReferenceValues.size(),
3471 temperatureCouplingTauValues.size());
3474 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3475 do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
3476 SimulationAtomGroupType::TemperatureCoupling, restnm,
3477 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3478 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3480 snew(ir->opts.nrdf, nr);
3481 snew(ir->opts.tau_t, nr);
3482 snew(ir->opts.ref_t, nr);
3483 if (ir->eI == eiBD && ir->bd_fric == 0)
3485 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3488 if (useReferenceTemperature)
3490 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3492 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3496 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3497 for (i = 0; (i < nr); i++)
3499 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3501 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3502 warning_error(wi, warn_buf);
3505 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3509 "tau-t = -1 is the value to signal that a group should not have "
3510 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3513 if (ir->opts.tau_t[i] >= 0)
3515 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3518 if (ir->etc != etcNO && ir->nsttcouple == -1)
3520 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3525 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3528 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3529 "md-vv; use either vrescale temperature with berendsen pressure or "
3530 "Nose-Hoover temperature with MTTK pressure");
3532 if (ir->epc == epcMTTK)
3534 if (ir->etc != etcNOSEHOOVER)
3537 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3542 if (ir->nstpcouple != ir->nsttcouple)
3544 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3545 ir->nstpcouple = ir->nsttcouple = mincouple;
3547 "for current Trotter decomposition methods with vv, nsttcouple and "
3548 "nstpcouple must be equal. Both have been reset to "
3549 "min(nsttcouple,nstpcouple) = %d",
3551 warning_note(wi, warn_buf);
3556 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3557 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3559 if (ETC_ANDERSEN(ir->etc))
3561 if (ir->nsttcouple != 1)
3565 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3566 "need for larger nsttcouple > 1, since no global parameters are computed. "
3567 "nsttcouple has been reset to 1");
3568 warning_note(wi, warn_buf);
3571 nstcmin = tcouple_min_integration_steps(ir->etc);
3574 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3577 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3578 "least %d times larger than nsttcouple*dt (%g)",
3579 ETCOUPLTYPE(ir->etc), tau_min, nstcmin, ir->nsttcouple * ir->delta_t);
3580 warning(wi, warn_buf);
3583 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3584 for (i = 0; (i < nr); i++)
3586 if (ir->opts.ref_t[i] < 0)
3588 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3591 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3592 if we are in this conditional) if mc_temp is negative */
3593 if (ir->expandedvals->mc_temp < 0)
3595 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3599 /* Simulated annealing for each group. There are nr groups */
3600 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3601 if (simulatedAnnealingGroupNames.size() == 1
3602 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3604 simulatedAnnealingGroupNames.resize(0);
3606 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3608 gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
3609 simulatedAnnealingGroupNames.size(), nr);
3613 snew(ir->opts.annealing, nr);
3614 snew(ir->opts.anneal_npoints, nr);
3615 snew(ir->opts.anneal_time, nr);
3616 snew(ir->opts.anneal_temp, nr);
3617 for (i = 0; i < nr; i++)
3619 ir->opts.annealing[i] = eannNO;
3620 ir->opts.anneal_npoints[i] = 0;
3621 ir->opts.anneal_time[i] = nullptr;
3622 ir->opts.anneal_temp[i] = nullptr;
3624 if (!simulatedAnnealingGroupNames.empty())
3627 for (i = 0; i < nr; i++)
3629 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3631 ir->opts.annealing[i] = eannNO;
3633 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3635 ir->opts.annealing[i] = eannSINGLE;
3638 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3640 ir->opts.annealing[i] = eannPERIODIC;
3646 /* Read the other fields too */
3647 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3648 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3650 gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
3651 simulatedAnnealingPoints.size(), simulatedAnnealingGroupNames.size());
3653 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3654 size_t numSimulatedAnnealingFields = 0;
3655 for (i = 0; i < nr; i++)
3657 if (ir->opts.anneal_npoints[i] == 1)
3661 "Please specify at least a start and an end point for annealing\n");
3663 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3664 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3665 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3668 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3670 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3672 gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
3673 simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
3675 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3676 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3678 gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
3679 simulatedAnnealingTemperatures.size(), numSimulatedAnnealingFields);
3682 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3683 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3684 convertReals(wi, simulatedAnnealingTimes, "anneal-time",
3685 allSimulatedAnnealingTimes.data());
3686 convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp",
3687 allSimulatedAnnealingTemperatures.data());
3688 for (i = 0, k = 0; i < nr; i++)
3690 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3692 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3693 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3696 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3698 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3704 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3707 "Annealing timepoints out of order: t=%f comes after "
3709 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j - 1]);
3712 if (ir->opts.anneal_temp[i][j] < 0)
3714 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n",
3715 ir->opts.anneal_temp[i][j]);
3720 /* Print out some summary information, to make sure we got it right */
3721 for (i = 0; i < nr; i++)
3723 if (ir->opts.annealing[i] != eannNO)
3725 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3726 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3727 *(groups->groupNames[j]), eann_names[ir->opts.annealing[i]],
3728 ir->opts.anneal_npoints[i]);
3729 fprintf(stderr, "Time (ps) Temperature (K)\n");
3730 /* All terms except the last one */
3731 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3733 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3734 ir->opts.anneal_temp[i][j]);
3737 /* Finally the last one */
3738 j = ir->opts.anneal_npoints[i] - 1;
3739 if (ir->opts.annealing[i] == eannSINGLE)
3741 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j],
3742 ir->opts.anneal_temp[i][j]);
3746 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3747 ir->opts.anneal_temp[i][j]);
3748 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3751 "There is a temperature jump when your annealing "
3763 process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3765 checkPullCoords(ir->pull->group, ir->pull->coord);
3770 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3773 if (ir->eSwapCoords != eswapNO)
3775 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3778 /* Make indices for IMD session */
3781 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3784 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3785 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3786 notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
3788 auto accelerations = gmx::splitString(inputrecStrings->acc);
3789 auto accelerationGroupNames = gmx::splitString(inputrecStrings->accgrps);
3790 if (accelerationGroupNames.size() * DIM != accelerations.size())
3792 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3793 accelerationGroupNames.size(), accelerations.size());
3795 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3796 SimulationAtomGroupType::Acceleration, restnm, egrptpALL_GENREST, bVerbose, wi);
3797 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3798 snew(ir->opts.acc, nr);
3799 ir->opts.ngacc = nr;
3801 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3803 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3804 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3805 if (freezeDims.size() != DIM * freezeGroupNames.size())
3807 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3808 freezeGroupNames.size(), freezeDims.size());
3810 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3811 SimulationAtomGroupType::Freeze, restnm, egrptpALL_GENREST, bVerbose, wi);
3812 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3813 ir->opts.ngfrz = nr;
3814 snew(ir->opts.nFreeze, nr);
3815 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3817 for (j = 0; (j < DIM); j++, k++)
3819 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3820 if (!ir->opts.nFreeze[i][j])
3822 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3825 "Please use Y(ES) or N(O) for freezedim only "
3827 freezeDims[k].c_str());
3828 warning(wi, warn_buf);
3833 for (; (i < nr); i++)
3835 for (j = 0; (j < DIM); j++)
3837 ir->opts.nFreeze[i][j] = 0;
3841 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3842 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3843 SimulationAtomGroupType::EnergyOutput, restnm, egrptpALL_GENREST, bVerbose, wi);
3844 add_wall_energrps(groups, ir->nwall, symtab);
3845 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3846 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3847 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3848 SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
3849 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3851 if (ir->comm_mode != ecmNO)
3853 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3856 /* Now we have filled the freeze struct, so we can calculate NRDF */
3857 calc_nrdf(mtop, ir, gnames);
3859 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3860 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3861 SimulationAtomGroupType::User1, restnm, egrptpALL_GENREST, bVerbose, wi);
3862 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3863 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3864 SimulationAtomGroupType::User2, restnm, egrptpALL_GENREST, bVerbose, wi);
3865 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3866 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3867 SimulationAtomGroupType::CompressedPositionOutput, restnm, egrptpONE, bVerbose, wi);
3868 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3869 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3870 SimulationAtomGroupType::OrientationRestraintsFit, restnm, egrptpALL_GENREST,
3873 /* MiMiC QMMM input processing */
3874 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
3875 if (qmGroupNames.size() > 1)
3877 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3879 /* group rest, if any, is always MM! */
3880 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3881 SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
3882 ir->opts.ngQM = qmGroupNames.size();
3884 /* end of MiMiC QMMM input */
3888 for (auto group : gmx::keysOf(groups->groups))
3890 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3891 for (const auto& entry : groups->groups[group])
3893 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3895 fprintf(stderr, "\n");
3899 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3900 snew(ir->opts.egp_flags, nr * nr);
3902 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
3903 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3905 warning_error(wi, "Energy group exclusions are currently not supported");
3907 if (bExcl && EEL_FULL(ir->coulombtype))
3909 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3912 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
3913 if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
3914 && !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
3917 "Can only have energy group pair tables in combination with user tables for VdW "
3921 /* final check before going out of scope if simulated tempering variables
3922 * need to be set to default values.
3924 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3926 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
3927 warning(wi, gmx::formatString(
3928 "the value for nstexpanded was not specified for "
3929 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3930 "by default, but it is recommended to set it to an explicit value!",
3931 ir->expandedvals->nstexpanded));
3933 for (i = 0; (i < defaultIndexGroups->nr); i++)
3938 done_blocka(defaultIndexGroups);
3939 sfree(defaultIndexGroups);
3943 static void check_disre(const gmx_mtop_t* mtop)
3945 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3947 const gmx_ffparams_t& ffparams = mtop->ffparams;
3950 for (int i = 0; i < ffparams.numTypes(); i++)
3952 int ftype = ffparams.functype[i];
3953 if (ftype == F_DISRES)
3955 int label = ffparams.iparams[i].disres.label;
3956 if (label == old_label)
3958 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3967 "Found %d double distance restraint indices,\n"
3968 "probably the parameters for multiple pairs in one restraint "
3969 "are not identical\n",
3975 static bool absolute_reference(const t_inputrec* ir, const gmx_mtop_t* sys, const bool posres_only, ivec AbsRef)
3978 gmx_mtop_ilistloop_t iloop;
3980 const t_iparams* pr;
3987 for (d = 0; d < DIM; d++)
3989 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3991 /* Check for freeze groups */
3992 for (g = 0; g < ir->opts.ngfrz; g++)
3994 for (d = 0; d < DIM; d++)
3996 if (ir->opts.nFreeze[g][d] != 0)
4004 /* Check for position restraints */
4005 iloop = gmx_mtop_ilistloop_init(sys);
4006 while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
4008 if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
4010 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
4012 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
4013 for (d = 0; d < DIM; d++)
4015 if (pr->posres.fcA[d] != 0)
4021 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
4023 /* Check for flat-bottom posres */
4024 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
4025 if (pr->fbposres.k != 0)
4027 switch (pr->fbposres.geom)
4029 case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
4030 case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
4031 case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
4032 case efbposresCYLINDER:
4033 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4034 case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
4035 case efbposresX: /* d=XX */
4036 case efbposresY: /* d=YY */
4037 case efbposresZ: /* d=ZZ */
4038 d = pr->fbposres.geom - efbposresX;
4043 " Invalid geometry for flat-bottom position restraint.\n"
4044 "Expected nr between 1 and %d. Found %d\n",
4045 efbposresNR - 1, pr->fbposres.geom);
4052 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
4055 static void check_combination_rule_differences(const gmx_mtop_t* mtop,
4057 bool* bC6ParametersWorkWithGeometricRules,
4058 bool* bC6ParametersWorkWithLBRules,
4059 bool* bLBRulesPossible)
4061 int ntypes, tpi, tpj;
4064 double c6i, c6j, c12i, c12j;
4065 double c6, c6_geometric, c6_LB;
4066 double sigmai, sigmaj, epsi, epsj;
4067 bool bCanDoLBRules, bCanDoGeometricRules;
4070 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4071 * force-field floating point parameters.
4074 ptr = getenv("GMX_LJCOMB_TOL");
4078 double gmx_unused canary;
4080 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4083 "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4088 *bC6ParametersWorkWithLBRules = TRUE;
4089 *bC6ParametersWorkWithGeometricRules = TRUE;
4090 bCanDoLBRules = TRUE;
4091 ntypes = mtop->ffparams.atnr;
4092 snew(typecount, ntypes);
4093 gmx_mtop_count_atomtypes(mtop, state, typecount);
4094 *bLBRulesPossible = TRUE;
4095 for (tpi = 0; tpi < ntypes; ++tpi)
4097 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4098 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4099 for (tpj = tpi; tpj < ntypes; ++tpj)
4101 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4102 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4103 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4104 c6_geometric = std::sqrt(c6i * c6j);
4105 if (!gmx_numzero(c6_geometric))
4107 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4109 sigmai = gmx::sixthroot(c12i / c6i);
4110 sigmaj = gmx::sixthroot(c12j / c6j);
4111 epsi = c6i * c6i / (4.0 * c12i);
4112 epsj = c6j * c6j / (4.0 * c12j);
4113 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4117 *bLBRulesPossible = FALSE;
4118 c6_LB = c6_geometric;
4120 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4125 *bC6ParametersWorkWithLBRules = FALSE;
4128 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4130 if (!bCanDoGeometricRules)
4132 *bC6ParametersWorkWithGeometricRules = FALSE;
4139 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop, warninp_t wi)
4141 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4143 check_combination_rule_differences(mtop, 0, &bC6ParametersWorkWithGeometricRules,
4144 &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4145 if (ir->ljpme_combination_rule == eljpmeLB)
4147 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4150 "You are using arithmetic-geometric combination rules "
4151 "in LJ-PME, but your non-bonded C6 parameters do not "
4152 "follow these rules.");
4157 if (!bC6ParametersWorkWithGeometricRules)
4159 if (ir->eDispCorr != edispcNO)
4162 "You are using geometric combination rules in "
4163 "LJ-PME, but your non-bonded C6 parameters do "
4164 "not follow these rules. "
4165 "This will introduce very small errors in the forces and energies in "
4166 "your simulations. Dispersion correction will correct total energy "
4167 "and/or pressure for isotropic systems, but not forces or surface "
4173 "You are using geometric combination rules in "
4174 "LJ-PME, but your non-bonded C6 parameters do "
4175 "not follow these rules. "
4176 "This will introduce very small errors in the forces and energies in "
4177 "your simulations. If your system is homogeneous, consider using "
4178 "dispersion correction "
4179 "for the total energy and pressure.");
4185 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4187 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4188 gmx::assertMtsRequirements(*ir);
4190 char err_buf[STRLEN];
4195 gmx_mtop_atomloop_block_t aloopb;
4197 char warn_buf[STRLEN];
4199 set_warning_line(wi, mdparin, -1);
4201 if (absolute_reference(ir, sys, false, AbsRef))
4204 "Removing center of mass motion in the presence of position restraints might "
4205 "cause artifacts. When you are using position restraints to equilibrate a "
4206 "macro-molecule, the artifacts are usually negligible.");
4209 if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 && ir->nstlist > 1
4210 && ((EI_MD(ir->eI) || EI_SD(ir->eI)) && (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4212 /* Check if a too small Verlet buffer might potentially
4213 * cause more drift than the thermostat can couple off.
4215 /* Temperature error fraction for warning and suggestion */
4216 const real T_error_warn = 0.002;
4217 const real T_error_suggest = 0.001;
4218 /* For safety: 2 DOF per atom (typical with constraints) */
4219 const real nrdf_at = 2;
4220 real T, tau, max_T_error;
4225 for (i = 0; i < ir->opts.ngtc; i++)
4227 T = std::max(T, ir->opts.ref_t[i]);
4228 tau = std::max(tau, ir->opts.tau_t[i]);
4232 /* This is a worst case estimate of the temperature error,
4233 * assuming perfect buffer estimation and no cancelation
4234 * of errors. The factor 0.5 is because energy distributes
4235 * equally over Ekin and Epot.
4237 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * BOLTZ * T);
4238 if (max_T_error > T_error_warn)
4241 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4242 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4243 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4244 "%.0e or decrease tau_t.",
4245 ir->verletbuf_tol, T, tau, 100 * max_T_error, 100 * T_error_suggest,
4246 ir->verletbuf_tol * T_error_suggest / max_T_error);
4247 warning(wi, warn_buf);
4252 if (ETC_ANDERSEN(ir->etc))
4256 for (i = 0; i < ir->opts.ngtc; i++)
4259 "all tau_t must currently be equal using Andersen temperature control, "
4260 "violated for group %d",
4262 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4264 "all tau_t must be positive using Andersen temperature control, "
4266 i, ir->opts.tau_t[i]);
4267 CHECK(ir->opts.tau_t[i] < 0);
4270 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
4272 for (i = 0; i < ir->opts.ngtc; i++)
4274 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4276 "tau_t/delta_t for group %d for temperature control method %s must be a "
4277 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4278 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4280 i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
4281 CHECK(nsteps % ir->nstcomm != 0);
4286 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD && ir->comm_mode == ecmNO
4287 && !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
4290 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4291 "rounding errors can lead to build up of kinetic energy of the center of mass");
4294 if (ir->epc == epcPARRINELLORAHMAN && ir->etc == etcNOSEHOOVER)
4297 for (int g = 0; g < ir->opts.ngtc; g++)
4299 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4301 if (ir->tau_p < 1.9 * tau_t_max)
4303 std::string message = gmx::formatString(
4304 "With %s T-coupling and %s p-coupling, "
4305 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4306 etcoupl_names[ir->etc], epcoupl_names[ir->epc], "tau-p", ir->tau_p, "tau-t",
4308 warning(wi, message.c_str());
4312 /* Check for pressure coupling with absolute position restraints */
4313 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4315 absolute_reference(ir, sys, TRUE, AbsRef);
4317 for (m = 0; m < DIM; m++)
4319 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4322 "You are using pressure coupling with absolute position restraints, "
4323 "this will give artifacts. Use the refcoord_scaling option.");
4331 aloopb = gmx_mtop_atomloop_block_init(sys);
4333 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4335 if (atom->q != 0 || atom->qB != 0)
4343 if (EEL_FULL(ir->coulombtype))
4346 "You are using full electrostatics treatment %s for a system without charges.\n"
4347 "This costs a lot of performance for just processing zeros, consider using %s "
4349 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4350 warning(wi, err_buf);
4355 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4358 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4359 "You might want to consider using %s electrostatics.\n",
4361 warning_note(wi, err_buf);
4365 /* Check if combination rules used in LJ-PME are the same as in the force field */
4366 if (EVDW_PME(ir->vdwtype))
4368 check_combination_rules(ir, sys, wi);
4371 /* Generalized reaction field */
4372 if (ir->coulombtype == eelGRF_NOTUSED)
4375 "Generalized reaction-field electrostatics is no longer supported. "
4376 "You can use normal reaction-field instead and compute the reaction-field "
4377 "constant by hand.");
4381 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4383 for (m = 0; (m < DIM); m++)
4385 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4394 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4395 for (const AtomProxy atomP : AtomRange(*sys))
4397 const t_atom& local = atomP.atom();
4398 int i = atomP.globalAtomNumber();
4399 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4402 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4404 for (m = 0; (m < DIM); m++)
4406 acc[m] += ir->opts.acc[i][m] * mgrp[i];
4410 for (m = 0; (m < DIM); m++)
4412 if (fabs(acc[m]) > 1e-6)
4414 const char* dim[DIM] = { "X", "Y", "Z" };
4415 fprintf(stderr, "Net Acceleration in %s direction, will %s be corrected\n", dim[m],
4416 ir->nstcomm != 0 ? "" : "not");
4417 if (ir->nstcomm != 0 && m < ndof_com(ir))
4421 (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4423 ir->opts.acc[i][m] -= acc[m];
4431 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0
4432 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4434 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4442 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4444 if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
4446 absolute_reference(ir, sys, FALSE, AbsRef);
4447 for (m = 0; m < DIM; m++)
4449 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4452 "You are using an absolute reference for pulling, but the rest of "
4453 "the system does not have an absolute reference. This will lead to "
4462 for (i = 0; i < 3; i++)
4464 for (m = 0; m <= i; m++)
4466 if ((ir->epc != epcNO && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4468 for (c = 0; c < ir->pull->ncoord; c++)
4470 if (ir->pull->coord[c].eGeom == epullgDIRPBC && ir->pull->coord[c].vec[m] != 0)
4473 "Can not have dynamic box while using pull geometry '%s' "
4475 EPULLGEOM(ir->pull->coord[c].eGeom), 'x' + m);
4486 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4488 char warn_buf[STRLEN];
4491 ptr = check_box(ir->pbcType, box);
4494 warning_error(wi, ptr);
4497 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4499 if (ir->shake_tol <= 0.0)
4501 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4502 warning_error(wi, warn_buf);
4506 if ((ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4508 /* If we have Lincs constraints: */
4509 if (ir->eI == eiMD && ir->etc == etcNO && ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4512 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4513 warning_note(wi, warn_buf);
4516 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4519 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4521 warning_note(wi, warn_buf);
4523 if (ir->epc == epcMTTK)
4525 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4529 if (bHasAnyConstraints && ir->epc == epcMTTK)
4531 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4534 if (ir->LincsWarnAngle > 90.0)
4536 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4537 warning(wi, warn_buf);
4538 ir->LincsWarnAngle = 90.0;
4541 if (ir->pbcType != PbcType::No)
4543 if (ir->nstlist == 0)
4546 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4547 "atoms might cause the simulation to crash.");
4549 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4552 "ERROR: The cut-off length is longer than half the shortest box vector or "
4553 "longer than the smallest box diagonal element. Increase the box size or "
4554 "decrease rlist.\n");
4555 warning_error(wi, warn_buf);