2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, The GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
52 #include "gromacs/applied_forces/awh/read_params.h"
53 #include "gromacs/fileio/readinp.h"
54 #include "gromacs/fileio/warninp.h"
55 #include "gromacs/gmxlib/network.h"
56 #include "gromacs/gmxpreprocess/toputil.h"
57 #include "gromacs/math/functions.h"
58 #include "gromacs/math/units.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/mdlib/calc_verletbuf.h"
61 #include "gromacs/mdrun/mdmodules.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/mdtypes/multipletimestepping.h"
65 #include "gromacs/mdtypes/pull_params.h"
66 #include "gromacs/options/options.h"
67 #include "gromacs/options/treesupport.h"
68 #include "gromacs/pbcutil/pbc.h"
69 #include "gromacs/selection/indexutil.h"
70 #include "gromacs/topology/block.h"
71 #include "gromacs/topology/ifunc.h"
72 #include "gromacs/topology/index.h"
73 #include "gromacs/topology/mtop_util.h"
74 #include "gromacs/topology/symtab.h"
75 #include "gromacs/topology/topology.h"
76 #include "gromacs/utility/arrayref.h"
77 #include "gromacs/utility/cstringutil.h"
78 #include "gromacs/utility/exceptions.h"
79 #include "gromacs/utility/fatalerror.h"
80 #include "gromacs/utility/filestream.h"
81 #include "gromacs/utility/gmxassert.h"
82 #include "gromacs/utility/ikeyvaluetreeerror.h"
83 #include "gromacs/utility/keyvaluetree.h"
84 #include "gromacs/utility/keyvaluetreebuilder.h"
85 #include "gromacs/utility/keyvaluetreemdpwriter.h"
86 #include "gromacs/utility/keyvaluetreetransform.h"
87 #include "gromacs/utility/mdmodulenotification.h"
88 #include "gromacs/utility/smalloc.h"
89 #include "gromacs/utility/strconvert.h"
90 #include "gromacs/utility/stringcompare.h"
91 #include "gromacs/utility/stringutil.h"
92 #include "gromacs/utility/textwriter.h"
97 /* Resource parameters
98 * Do not change any of these until you read the instruction
99 * in readinp.h. Some cpp's do not take spaces after the backslash
100 * (like the c-shell), which will give you a very weird compiler
104 struct gmx_inputrec_strings
106 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], freeze[STRLEN], frdim[STRLEN],
107 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
108 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
109 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
110 char fep_lambda[efptNR][STRLEN];
111 char lambda_weights[STRLEN];
112 std::vector<std::string> pullGroupNames;
113 std::vector<std::string> rotateGroupNames;
114 char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
117 static gmx_inputrec_strings* inputrecStrings = nullptr;
119 void init_inputrec_strings()
124 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
125 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
127 inputrecStrings = new gmx_inputrec_strings();
130 void done_inputrec_strings()
132 delete inputrecStrings;
133 inputrecStrings = nullptr;
139 egrptpALL, /* All particles have to be a member of a group. */
140 egrptpALL_GENREST, /* A rest group with name is generated for particles *
141 * that are not part of any group. */
142 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
143 * for the rest group. */
144 egrptpONE /* Merge all selected groups into one group, *
145 * make a rest group for the remaining particles. */
148 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
149 "h-angles", "all-angles", nullptr };
151 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
153 static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lambdas)
158 for (i = 0; i < ntemps; i++)
160 /* simple linear scaling -- allows more control */
161 if (simtemp->eSimTempScale == esimtempLINEAR)
163 simtemp->temperatures[i] =
165 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
167 else if (simtemp->eSimTempScale
168 == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
170 simtemp->temperatures[i] = simtemp->simtemp_low
171 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
172 static_cast<real>((1.0 * i) / (ntemps - 1)));
174 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
176 simtemp->temperatures[i] = simtemp->simtemp_low
177 + (simtemp->simtemp_high - simtemp->simtemp_low)
178 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
183 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
184 gmx_fatal(FARGS, "%s", errorstr);
190 static void _low_check(bool b, const char* s, warninp_t wi)
194 warning_error(wi, s);
198 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
202 if (*p > 0 && *p % nst != 0)
204 /* Round up to the next multiple of nst */
205 *p = ((*p) / nst + 1) * nst;
206 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
211 //! Convert legacy mdp entries to modern ones.
212 static void process_interaction_modifier(int* eintmod)
214 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
216 *eintmod = eintmodPOTSHIFT;
220 void check_ir(const char* mdparin,
221 const gmx::MdModulesNotifier& mdModulesNotifier,
225 /* Check internal consistency.
226 * NOTE: index groups are not set here yet, don't check things
227 * like temperature coupling group options here, but in triple_check
230 /* Strange macro: first one fills the err_buf, and then one can check
231 * the condition, which will print the message and increase the error
234 #define CHECK(b) _low_check(b, err_buf, wi)
235 char err_buf[256], warn_buf[STRLEN];
238 t_lambda* fep = ir->fepvals;
239 t_expanded* expand = ir->expandedvals;
241 set_warning_line(wi, mdparin, -1);
243 /* We cannot check MTS requirements with an invalid MTS setup
244 * and we will already have generated errors with an invalid MTS setup.
246 if (gmx::haveValidMtsSetup(*ir))
248 std::vector<std::string> errorMessages = gmx::checkMtsRequirements(*ir);
250 for (const auto& errorMessage : errorMessages)
252 warning_error(wi, errorMessage.c_str());
256 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
258 sprintf(warn_buf, "%s electrostatics is no longer supported", eel_names[eelRF_NEC_UNSUPPORTED]);
259 warning_error(wi, warn_buf);
262 /* BASIC CUT-OFF STUFF */
263 if (ir->rcoulomb < 0)
265 warning_error(wi, "rcoulomb should be >= 0");
269 warning_error(wi, "rvdw should be >= 0");
271 if (ir->rlist < 0 && !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
273 warning_error(wi, "rlist should be >= 0");
276 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
277 "neighbour-list update scheme for efficient buffering for improved energy "
278 "conservation, please use the Verlet cut-off scheme instead.)");
279 CHECK(ir->nstlist < 0);
281 process_interaction_modifier(&ir->coulomb_modifier);
282 process_interaction_modifier(&ir->vdw_modifier);
284 if (ir->cutoff_scheme == ecutsGROUP)
287 "The group cutoff scheme has been removed since GROMACS 2020. "
288 "Please use the Verlet cutoff scheme.");
290 if (ir->cutoff_scheme == ecutsVERLET)
294 /* Normal Verlet type neighbor-list, currently only limited feature support */
295 if (inputrec2nboundeddim(ir) < 3)
297 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
300 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
301 if (ir->rcoulomb != ir->rvdw)
303 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
304 // for PME load balancing, we can support this exception.
305 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == evdwCUT
306 && ir->rcoulomb > ir->rvdw);
307 if (!bUsesPmeTwinRangeKernel)
310 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
311 "rcoulomb>rvdw with PME electrostatics)");
315 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
317 if (ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT)
319 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
322 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
324 evdw_names[ir->vdwtype],
326 eintmod_names[ir->vdw_modifier]);
327 warning_note(wi, warn_buf);
329 ir->vdwtype = evdwCUT;
334 "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
335 evdw_names[ir->vdwtype],
336 eintmod_names[ir->vdw_modifier]);
337 warning_error(wi, warn_buf);
341 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
344 "With Verlet lists only cut-off and PME LJ interactions are supported");
346 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) || EEL_PME(ir->coulombtype)
347 || ir->coulombtype == eelEWALD))
350 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
351 "electrostatics are supported");
353 if (!(ir->coulomb_modifier == eintmodNONE || ir->coulomb_modifier == eintmodPOTSHIFT))
355 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
356 warning_error(wi, warn_buf);
359 if (EEL_USER(ir->coulombtype))
362 "Coulomb type %s is not supported with the verlet scheme",
363 eel_names[ir->coulombtype]);
364 warning_error(wi, warn_buf);
367 if (ir->nstlist <= 0)
369 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
372 if (ir->nstlist < 10)
375 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
376 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
380 rc_max = std::max(ir->rvdw, ir->rcoulomb);
384 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
385 ir->verletbuf_tol = 0;
388 else if (ir->verletbuf_tol <= 0)
390 if (ir->verletbuf_tol == 0)
392 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
395 if (ir->rlist < rc_max)
398 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
401 if (ir->rlist == rc_max && ir->nstlist > 1)
405 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
406 "buffer. The cluster pair list does have a buffering effect, but choosing "
407 "a larger rlist might be necessary for good energy conservation.");
412 if (ir->rlist > rc_max)
415 "You have set rlist larger than the interaction cut-off, but you also "
416 "have verlet-buffer-tolerance > 0. Will set rlist using "
417 "verlet-buffer-tolerance.");
420 if (ir->nstlist == 1)
422 /* No buffer required */
427 if (EI_DYNAMICS(ir->eI))
429 if (inputrec2nboundeddim(ir) < 3)
432 "The box volume is required for calculating rlist from the "
433 "energy drift with verlet-buffer-tolerance > 0. You are "
434 "using at least one unbounded dimension, so no volume can be "
435 "computed. Either use a finite box, or set rlist yourself "
436 "together with verlet-buffer-tolerance = -1.");
438 /* Set rlist temporarily so we can continue processing */
443 /* Set the buffer to 5% of the cut-off */
444 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
450 /* GENERAL INTEGRATOR STUFF */
453 if (ir->etc != TemperatureCoupling::No)
455 if (EI_RANDOM(ir->eI))
458 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
459 "implicitly. See the documentation for more information on which "
460 "parameters affect temperature for %s.",
461 enumValueToString(ir->etc),
468 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
470 enumValueToString(ir->etc),
473 warning_note(wi, warn_buf);
475 ir->etc = TemperatureCoupling::No;
477 if (ir->eI == eiVVAK)
480 "Integrator method %s is implemented primarily for validation purposes; for "
481 "molecular dynamics, you should probably be using %s or %s",
485 warning_note(wi, warn_buf);
487 if (!EI_DYNAMICS(ir->eI))
489 if (ir->epc != PressureCoupling::No)
492 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
493 enumValueToString(ir->epc),
495 warning_note(wi, warn_buf);
497 ir->epc = PressureCoupling::No;
499 if (EI_DYNAMICS(ir->eI))
501 if (ir->nstcalcenergy < 0)
503 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
504 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
506 /* nstcalcenergy larger than nstener does not make sense.
507 * We ideally want nstcalcenergy=nstener.
511 ir->nstcalcenergy = std::gcd(ir->nstenergy, ir->nstlist);
515 ir->nstcalcenergy = ir->nstenergy;
519 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
520 || (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
521 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
524 const char* nsten = "nstenergy";
525 const char* nstdh = "nstdhdl";
526 const char* min_name = nsten;
527 int min_nst = ir->nstenergy;
529 /* find the smallest of ( nstenergy, nstdhdl ) */
530 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
531 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
533 min_nst = ir->fepvals->nstdhdl;
536 /* If the user sets nstenergy small, we should respect that */
537 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy, min_name, min_nst);
538 warning_note(wi, warn_buf);
539 ir->nstcalcenergy = min_nst;
542 if (ir->epc != PressureCoupling::No)
544 if (ir->nstpcouple < 0)
546 ir->nstpcouple = ir_optimal_nstpcouple(ir);
548 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
551 "With multiple time stepping, nstpcouple should be a mutiple of "
556 if (ir->nstcalcenergy > 0)
558 if (ir->efep != efepNO)
560 /* nstdhdl should be a multiple of nstcalcenergy */
561 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
565 /* nstexpanded should be a multiple of nstcalcenergy */
566 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded", &ir->expandedvals->nstexpanded, wi);
568 /* for storing exact averages nstenergy should be
569 * a multiple of nstcalcenergy
571 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
574 // Inquire all MdModules, if their parameters match with the energy
575 // calculation frequency
576 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
577 mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
579 // Emit all errors from the energy calculation frequency checks
580 for (const std::string& energyFrequencyErrorMessage :
581 energyCalculationFrequencyErrors.errorMessages())
583 warning_error(wi, energyFrequencyErrorMessage);
587 if (ir->nsteps == 0 && !ir->bContinuation)
590 "For a correct single-point energy evaluation with nsteps = 0, use "
591 "continuation = yes to avoid constraining the input coordinates.");
595 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
598 "You are doing a continuation with SD or BD, make sure that ld_seed is "
599 "different from the previous run (using ld_seed=-1 will ensure this)");
605 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
606 CHECK(ir->pbcType != PbcType::Xyz);
607 sprintf(err_buf, "with TPI nstlist should be larger than zero");
608 CHECK(ir->nstlist <= 0);
609 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
610 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
614 if ((opts->nshake > 0) && (opts->bMorse))
616 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
617 warning(wi, warn_buf);
620 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
623 "You are doing a continuation with SD or BD, make sure that ld_seed is "
624 "different from the previous run (using ld_seed=-1 will ensure this)");
626 /* verify simulated tempering options */
630 bool bAllTempZero = TRUE;
631 for (i = 0; i < fep->n_lambda; i++)
634 "Entry %d for %s must be between 0 and 1, instead is %g",
636 efpt_names[efptTEMPERATURE],
637 fep->all_lambda[efptTEMPERATURE][i]);
638 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
639 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
641 bAllTempZero = FALSE;
644 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
645 CHECK(bAllTempZero == TRUE);
647 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
648 CHECK(ir->eI != eiVV);
650 /* check compatability of the temperature coupling with simulated tempering */
652 if (ir->etc == TemperatureCoupling::NoseHoover)
655 "Nose-Hoover based temperature control such as [%s] my not be "
656 "entirelyconsistent with simulated tempering",
657 enumValueToString(ir->etc));
658 warning_note(wi, warn_buf);
661 /* check that the temperatures make sense */
664 "Higher simulated tempering temperature (%g) must be >= than the simulated "
665 "tempering lower temperature (%g)",
666 ir->simtempvals->simtemp_high,
667 ir->simtempvals->simtemp_low);
668 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
671 "Higher simulated tempering temperature (%g) must be >= zero",
672 ir->simtempvals->simtemp_high);
673 CHECK(ir->simtempvals->simtemp_high <= 0);
676 "Lower simulated tempering temperature (%g) must be >= zero",
677 ir->simtempvals->simtemp_low);
678 CHECK(ir->simtempvals->simtemp_low <= 0);
681 /* verify free energy options */
683 if (ir->efep != efepNO)
686 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
687 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
690 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
692 static_cast<int>(fep->sc_r_power));
693 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
696 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
699 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
702 "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
704 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
706 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
707 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
709 sprintf(err_buf, "Free-energy not implemented for Ewald");
710 CHECK(ir->coulombtype == eelEWALD);
712 /* check validty of lambda inputs */
713 if (fep->n_lambda == 0)
715 /* Clear output in case of no states:*/
716 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
717 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
722 "initial thermodynamic state %d does not exist, only goes to %d",
725 CHECK((fep->init_fep_state >= fep->n_lambda));
729 "Lambda state must be set, either with init-lambda-state or with init-lambda");
730 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
733 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
734 "init-lambda-state or with init-lambda, but not both",
736 fep->init_fep_state);
737 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
740 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
744 for (i = 0; i < efptNR; i++)
746 if (fep->separate_dvdl[i])
751 if (n_lambda_terms > 1)
754 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
755 "use init-lambda to set lambda state (except for slow growth). Use "
756 "init-lambda-state instead.");
757 warning(wi, warn_buf);
760 if (n_lambda_terms < 2 && fep->n_lambda > 0)
763 "init-lambda is deprecated for setting lambda state (except for slow "
764 "growth). Use init-lambda-state instead.");
768 for (j = 0; j < efptNR; j++)
770 for (i = 0; i < fep->n_lambda; i++)
773 "Entry %d for %s must be between 0 and 1, instead is %g",
776 fep->all_lambda[j][i]);
777 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
781 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
783 for (i = 0; i < fep->n_lambda; i++)
786 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
787 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
788 "crashes, and is not supported.",
790 fep->all_lambda[efptVDW][i],
791 fep->all_lambda[efptCOUL][i]);
792 CHECK((fep->sc_alpha > 0)
793 && (((fep->all_lambda[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL][i] < 1.0))
794 && ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0))));
798 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
800 real sigma, lambda, r_sc;
803 /* Maximum estimate for A and B charges equal with lambda power 1 */
805 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
806 1.0 / fep->sc_r_power);
808 "With PME there is a minor soft core effect present at the cut-off, "
809 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
810 "energy conservation, but usually other effects dominate. With a common sigma "
811 "value of %g nm the fraction of the particle-particle potential at the cut-off "
812 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
818 warning_note(wi, warn_buf);
821 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
822 be treated differently, but that's the next step */
824 for (i = 0; i < efptNR; i++)
826 for (j = 0; j < fep->n_lambda; j++)
828 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
829 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
834 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
838 /* checking equilibration of weights inputs for validity */
841 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
843 expand->equil_n_at_lam,
844 elmceq_names[elmceqNUMATLAM]);
845 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
848 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
850 expand->equil_samples,
851 elmceq_names[elmceqSAMPLES]);
852 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
855 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
857 elmceq_names[elmceqSTEPS]);
858 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
861 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
862 expand->equil_samples,
863 elmceq_names[elmceqWLDELTA]);
864 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
867 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
869 elmceq_names[elmceqRATIO]);
870 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
873 "weight-equil-number-all-lambda (%d) must be a positive integer if "
874 "lmc-weights-equil=%s",
875 expand->equil_n_at_lam,
876 elmceq_names[elmceqNUMATLAM]);
877 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
880 "weight-equil-number-samples (%d) must be a positive integer if "
881 "lmc-weights-equil=%s",
882 expand->equil_samples,
883 elmceq_names[elmceqSAMPLES]);
884 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
887 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
889 elmceq_names[elmceqSTEPS]);
890 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
893 "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
894 expand->equil_wl_delta,
895 elmceq_names[elmceqWLDELTA]);
896 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
899 "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
901 elmceq_names[elmceqRATIO]);
902 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
905 "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
906 elmceq_names[elmceqWLDELTA],
907 elamstats_names[elamstatsWL],
908 elamstats_names[elamstatsWWL]);
909 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
911 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
912 CHECK((expand->lmc_repeats <= 0));
913 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
914 CHECK((expand->minvarmin <= 0));
915 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
916 CHECK((expand->c_range < 0));
918 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
921 expand->lmc_forced_nstart);
922 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
923 && (expand->elmcmove != elmcmoveNO));
924 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
925 CHECK((expand->lmc_forced_nstart < 0));
927 "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
928 fep->init_fep_state);
929 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
931 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
932 CHECK((expand->init_wl_delta < 0));
933 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
934 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
935 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
936 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
938 /* if there is no temperature control, we need to specify an MC temperature */
939 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO)
940 && (expand->mc_temp <= 0.0))
943 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
944 "to a positive number");
945 warning_error(wi, err_buf);
947 if (expand->nstTij > 0)
949 sprintf(err_buf, "nstlog must be non-zero");
950 CHECK(ir->nstlog == 0);
951 // Avoid modulus by zero in the case that already triggered an error exit.
955 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
958 CHECK((expand->nstTij % ir->nstlog) != 0);
964 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
965 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
968 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
970 if (ir->pbcType == PbcType::No)
972 if (ir->epc != PressureCoupling::No)
974 warning(wi, "Turning off pressure coupling for vacuum system");
975 ir->epc = PressureCoupling::No;
981 "Can not have pressure coupling with pbc=%s",
982 c_pbcTypeNames[ir->pbcType].c_str());
983 CHECK(ir->epc != PressureCoupling::No);
985 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
986 CHECK(EEL_FULL(ir->coulombtype));
989 "Can not have dispersion correction with pbc=%s",
990 c_pbcTypeNames[ir->pbcType].c_str());
991 CHECK(ir->eDispCorr != edispcNO);
994 if (ir->rlist == 0.0)
997 "can only have neighborlist cut-off zero (=infinite)\n"
998 "with coulombtype = %s or coulombtype = %s\n"
999 "without periodic boundary conditions (pbc = %s) and\n"
1000 "rcoulomb and rvdw set to zero",
1003 c_pbcTypeNames[PbcType::No].c_str());
1004 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
1005 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1007 if (ir->nstlist > 0)
1010 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1011 "nstype=simple and only one MPI rank");
1016 if (ir->nstcomm == 0)
1018 // TODO Change this behaviour. There should be exactly one way
1019 // to turn off an algorithm.
1020 ir->comm_mode = ecmNO;
1022 if (ir->comm_mode != ecmNO)
1024 if (ir->nstcomm < 0)
1026 // TODO Such input was once valid. Now that we've been
1027 // helpful for a few years, we should reject such input,
1028 // lest we have to support every historical decision
1031 "If you want to remove the rotation around the center of mass, you should set "
1032 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1033 "its absolute value");
1034 ir->nstcomm = abs(ir->nstcomm);
1037 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
1040 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
1041 "nstcomm to nstcalcenergy");
1042 ir->nstcomm = ir->nstcalcenergy;
1045 if (ir->comm_mode == ecmANGULAR)
1048 "Can not remove the rotation around the center of mass with periodic "
1050 CHECK(ir->bPeriodicMols);
1051 if (ir->pbcType != PbcType::No)
1054 "Removing the rotation around the center of mass in a periodic system, "
1055 "this can lead to artifacts. Only use this on a single (cluster of) "
1056 "molecules. This cluster should not cross periodic boundaries.");
1061 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
1064 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1065 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1068 warning_note(wi, warn_buf);
1071 /* TEMPERATURE COUPLING */
1072 if (ir->etc == TemperatureCoupling::Yes)
1074 ir->etc = TemperatureCoupling::Berendsen;
1076 "Old option for temperature coupling given: "
1077 "changing \"yes\" to \"Berendsen\"\n");
1080 if ((ir->etc == TemperatureCoupling::NoseHoover) || (ir->epc == PressureCoupling::Mttk))
1082 if (ir->opts.nhchainlength < 1)
1085 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1087 ir->opts.nhchainlength);
1088 ir->opts.nhchainlength = 1;
1089 warning(wi, warn_buf);
1092 if (ir->etc == TemperatureCoupling::NoseHoover && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1096 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1097 ir->opts.nhchainlength = 1;
1102 ir->opts.nhchainlength = 0;
1105 if (ir->eI == eiVVAK)
1108 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1111 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1114 if (ETC_ANDERSEN(ir->etc))
1117 "%s temperature control not supported for integrator %s.",
1118 enumValueToString(ir->etc),
1120 CHECK(!(EI_VV(ir->eI)));
1122 if (ir->nstcomm > 0 && (ir->etc == TemperatureCoupling::Andersen))
1125 "Center of mass removal not necessary for %s. All velocities of coupled "
1126 "groups are rerandomized periodically, so flying ice cube errors will not "
1128 enumValueToString(ir->etc));
1129 warning_note(wi, warn_buf);
1133 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1134 "randomized every time step",
1136 enumValueToString(ir->etc));
1137 CHECK(ir->nstcomm > 1 && (ir->etc == TemperatureCoupling::Andersen));
1140 if (ir->etc == TemperatureCoupling::Berendsen)
1143 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1144 "might want to consider using the %s thermostat.",
1145 enumValueToString(ir->etc),
1146 enumValueToString(TemperatureCoupling::VRescale));
1147 warning_note(wi, warn_buf);
1150 if ((ir->etc == TemperatureCoupling::NoseHoover || ETC_ANDERSEN(ir->etc))
1151 && ir->epc == PressureCoupling::Berendsen)
1154 "Using Berendsen pressure coupling invalidates the "
1155 "true ensemble for the thermostat");
1156 warning(wi, warn_buf);
1159 /* PRESSURE COUPLING */
1160 if (ir->epc == PressureCoupling::Isotropic)
1162 ir->epc = PressureCoupling::Berendsen;
1164 "Old option for pressure coupling given: "
1165 "changing \"Isotropic\" to \"Berendsen\"\n");
1168 if (ir->epc != PressureCoupling::No)
1170 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1172 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1173 CHECK(ir->tau_p <= 0);
1175 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1178 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1179 "times larger than nstpcouple*dt (%g)",
1180 enumValueToString(ir->epc),
1182 pcouple_min_integration_steps(ir->epc),
1184 warning(wi, warn_buf);
1188 "compressibility must be > 0 when using pressure"
1190 enumValueToString(ir->epc));
1191 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1192 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1193 && ir->compress[ZZ][YY] <= 0));
1195 if (PressureCoupling::ParrinelloRahman == ir->epc && opts->bGenVel)
1198 "You are generating velocities so I am assuming you "
1199 "are equilibrating a system. You are using "
1200 "%s pressure coupling, but this can be "
1201 "unstable for equilibration. If your system crashes, try "
1202 "equilibrating first with Berendsen pressure coupling. If "
1203 "you are not equilibrating the system, you can probably "
1204 "ignore this warning.",
1205 enumValueToString(ir->epc));
1206 warning(wi, warn_buf);
1212 if (ir->epc == PressureCoupling::Mttk)
1214 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1218 /* ELECTROSTATICS */
1219 /* More checks are in triple check (grompp.c) */
1221 if (ir->coulombtype == eelSWITCH)
1224 "coulombtype = %s is only for testing purposes and can lead to serious "
1225 "artifacts, advice: use coulombtype = %s",
1226 eel_names[ir->coulombtype],
1227 eel_names[eelRF_ZERO]);
1228 warning(wi, warn_buf);
1231 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1234 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1235 "format and exchanging epsilon-r and epsilon-rf",
1237 warning(wi, warn_buf);
1238 ir->epsilon_rf = ir->epsilon_r;
1239 ir->epsilon_r = 1.0;
1242 if (ir->epsilon_r == 0)
1245 "It is pointless to use long-range electrostatics with infinite relative "
1247 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1249 CHECK(EEL_FULL(ir->coulombtype));
1252 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1254 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1255 CHECK(ir->epsilon_r < 0);
1258 if (EEL_RF(ir->coulombtype))
1260 /* reaction field (at the cut-off) */
1262 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1265 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1266 eel_names[ir->coulombtype]);
1267 warning(wi, warn_buf);
1268 ir->epsilon_rf = 0.0;
1271 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1272 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1273 if (ir->epsilon_rf == ir->epsilon_r)
1276 "Using epsilon-rf = epsilon-r with %s does not make sense",
1277 eel_names[ir->coulombtype]);
1278 warning(wi, warn_buf);
1281 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1282 * means the interaction is zero outside rcoulomb, but it helps to
1283 * provide accurate energy conservation.
1285 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1287 if (ir_coulomb_switched(ir))
1290 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1291 "potential modifier options!",
1292 eel_names[ir->coulombtype]);
1293 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1297 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1300 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1301 "secondary coulomb-modifier.");
1302 CHECK(ir->coulomb_modifier != eintmodNONE);
1304 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1307 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1308 "secondary vdw-modifier.");
1309 CHECK(ir->vdw_modifier != eintmodNONE);
1312 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->vdwtype == evdwSWITCH
1313 || ir->vdwtype == evdwSHIFT)
1316 "The switch/shift interaction settings are just for compatibility; you will get "
1318 "performance from applying potential modifiers to your interactions!\n");
1319 warning_note(wi, warn_buf);
1322 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1324 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1326 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1328 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1329 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1330 "will be good regardless, since ewald_rtol = %g.",
1332 ir->rcoulomb_switch,
1335 warning(wi, warn_buf);
1339 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1341 if (ir->rvdw_switch == 0)
1344 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1345 "potential. This suggests it was not set in the mdp, which can lead to large "
1346 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1347 "switching range.");
1348 warning(wi, warn_buf);
1352 if (EEL_FULL(ir->coulombtype))
1354 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER
1355 || ir->coulombtype == eelPMEUSERSWITCH)
1357 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist", eel_names[ir->coulombtype]);
1358 CHECK(ir->rcoulomb > ir->rlist);
1362 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1364 // TODO: Move these checks into the ewald module with the options class
1366 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1368 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1371 "With coulombtype = %s, you should have %d <= pme-order <= %d",
1372 eel_names[ir->coulombtype],
1375 warning_error(wi, warn_buf);
1379 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1381 if (ir->ewald_geometry == eewg3D)
1384 "With pbc=%s you should use ewald-geometry=%s",
1385 c_pbcTypeNames[ir->pbcType].c_str(),
1386 eewg_names[eewg3DC]);
1387 warning(wi, warn_buf);
1389 /* This check avoids extra pbc coding for exclusion corrections */
1390 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1391 CHECK(ir->wall_ewald_zfac < 2);
1393 if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
1396 "With %s and ewald_geometry = %s you should use pbc = %s",
1397 eel_names[ir->coulombtype],
1398 eewg_names[eewg3DC],
1399 c_pbcTypeNames[PbcType::XY].c_str());
1400 warning(wi, warn_buf);
1402 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1404 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1405 CHECK(ir->bPeriodicMols);
1406 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1407 warning_note(wi, warn_buf);
1409 "With epsilon_surface > 0 you can only use domain decomposition "
1410 "when there are only small molecules with all bonds constrained (mdrun will check "
1412 warning_note(wi, warn_buf);
1415 if (ir_vdw_switched(ir))
1417 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1418 CHECK(ir->rvdw_switch >= ir->rvdw);
1420 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1423 "You are applying a switch function to vdw forces or potentials from %g to %g "
1424 "nm, which is more than half the interaction range, whereas switch functions "
1425 "are intended to act only close to the cut-off.",
1428 warning_note(wi, warn_buf);
1432 if (ir->vdwtype == evdwPME)
1434 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1437 "With vdwtype = %s, the only supported modifiers are %s and %s",
1438 evdw_names[ir->vdwtype],
1439 eintmod_names[eintmodPOTSHIFT],
1440 eintmod_names[eintmodNONE]);
1441 warning_error(wi, err_buf);
1445 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1448 "You have selected user tables with dispersion correction, the dispersion "
1449 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1450 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1451 "really want dispersion correction to -C6/r^6.");
1454 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT) && ir->rvdw != 0)
1456 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1459 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1461 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1464 /* IMPLICIT SOLVENT */
1465 if (ir->coulombtype == eelGB_NOTUSED)
1467 sprintf(warn_buf, "Invalid option %s for coulombtype", eel_names[ir->coulombtype]);
1468 warning_error(wi, warn_buf);
1473 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1478 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1481 // cosine acceleration is only supported in leap-frog
1482 if (ir->cos_accel != 0.0 && ir->eI != eiMD)
1484 warning_error(wi, "cos-acceleration is only supported by integrator = md");
1488 /* interpret a number of doubles from a string and put them in an array,
1489 after allocating space for them.
1490 str = the input string
1491 n = the (pre-allocated) number of doubles read
1492 r = the output array of doubles. */
1493 static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
1495 auto values = gmx::splitString(str);
1499 for (int i = 0; i < *n; i++)
1503 (*r)[i] = gmx::fromString<real>(values[i]);
1505 catch (gmx::GromacsException&)
1508 "Invalid value " + values[i]
1509 + " in string in mdp file. Expected a real number.");
1515 static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1518 int i, j, max_n_lambda, nweights, nfep[efptNR];
1519 t_lambda* fep = ir->fepvals;
1520 t_expanded* expand = ir->expandedvals;
1521 real** count_fep_lambdas;
1522 bool bOneLambda = TRUE;
1524 snew(count_fep_lambdas, efptNR);
1526 /* FEP input processing */
1527 /* first, identify the number of lambda values for each type.
1528 All that are nonzero must have the same number */
1530 for (i = 0; i < efptNR; i++)
1532 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1535 /* now, determine the number of components. All must be either zero, or equal. */
1538 for (i = 0; i < efptNR; i++)
1540 if (nfep[i] > max_n_lambda)
1542 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1543 must have the same number if its not zero.*/
1548 for (i = 0; i < efptNR; i++)
1552 ir->fepvals->separate_dvdl[i] = FALSE;
1554 else if (nfep[i] == max_n_lambda)
1556 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute
1557 the derivative with respect to the temperature currently */
1559 ir->fepvals->separate_dvdl[i] = TRUE;
1565 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1572 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1573 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1575 /* the number of lambdas is the number we've read in, which is either zero
1576 or the same for all */
1577 fep->n_lambda = max_n_lambda;
1579 /* allocate space for the array of lambda values */
1580 snew(fep->all_lambda, efptNR);
1581 /* if init_lambda is defined, we need to set lambda */
1582 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1584 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1586 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1587 for (i = 0; i < efptNR; i++)
1589 snew(fep->all_lambda[i], fep->n_lambda);
1590 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1593 for (j = 0; j < fep->n_lambda; j++)
1595 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1597 sfree(count_fep_lambdas[i]);
1600 sfree(count_fep_lambdas);
1602 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1603 internal bookkeeping -- for now, init_lambda */
1605 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1607 for (i = 0; i < fep->n_lambda; i++)
1609 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1613 /* check to see if only a single component lambda is defined, and soft core is defined.
1614 In this case, turn on coulomb soft core */
1616 if (max_n_lambda == 0)
1622 for (i = 0; i < efptNR; i++)
1624 if ((nfep[i] != 0) && (i != efptFEP))
1630 if ((bOneLambda) && (fep->sc_alpha > 0))
1632 fep->bScCoul = TRUE;
1635 /* Fill in the others with the efptFEP if they are not explicitly
1636 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1637 they are all zero. */
1639 for (i = 0; i < efptNR; i++)
1641 if ((nfep[i] == 0) && (i != efptFEP))
1643 for (j = 0; j < fep->n_lambda; j++)
1645 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1651 /* now read in the weights */
1652 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1655 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1657 else if (nweights != fep->n_lambda)
1660 "Number of weights (%d) is not equal to number of lambda values (%d)",
1664 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1666 expand->nstexpanded = fep->nstdhdl;
1667 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1672 static void do_simtemp_params(t_inputrec* ir)
1675 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1676 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1679 template<typename T>
1680 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1683 for (const auto& input : inputs)
1687 outputs[i] = gmx::fromStdString<T>(input);
1689 catch (gmx::GromacsException&)
1691 auto message = gmx::formatString(
1692 "Invalid value for mdp option %s. %s should only consist of integers separated "
1696 warning_error(wi, message);
1702 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1705 for (const auto& input : inputs)
1709 outputs[i] = gmx::fromString<real>(input);
1711 catch (gmx::GromacsException&)
1713 auto message = gmx::formatString(
1714 "Invalid value for mdp option %s. %s should only consist of real numbers "
1715 "separated by spaces.",
1718 warning_error(wi, message);
1724 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1726 opts->wall_atomtype[0] = nullptr;
1727 opts->wall_atomtype[1] = nullptr;
1729 ir->wall_atomtype[0] = -1;
1730 ir->wall_atomtype[1] = -1;
1731 ir->wall_density[0] = 0;
1732 ir->wall_density[1] = 0;
1736 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1737 if (wallAtomTypes.size() != size_t(ir->nwall))
1740 "Expected %d elements for wall_atomtype, found %zu",
1742 wallAtomTypes.size());
1744 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1745 for (int i = 0; i < ir->nwall; i++)
1747 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1750 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1752 auto wallDensity = gmx::splitString(wall_density);
1753 if (wallDensity.size() != size_t(ir->nwall))
1756 "Expected %d elements for wall-density, found %zu",
1758 wallDensity.size());
1760 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1761 for (int i = 0; i < ir->nwall; i++)
1763 if (ir->wall_density[i] <= 0)
1765 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1772 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1776 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1777 for (int i = 0; i < nwall; i++)
1779 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1780 grps->emplace_back(groups->groupNames.size() - 1);
1785 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1787 /* read expanded ensemble parameters */
1788 printStringNewline(inp, "expanded ensemble variables");
1789 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1790 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1791 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1792 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1793 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1794 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1795 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1796 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1797 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1798 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1799 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1800 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1801 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1802 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1803 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1804 expand->bSymmetrizedTMatrix =
1805 (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1806 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1807 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1808 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1809 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1810 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1811 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1812 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1815 template<typename EnumType>
1816 EnumType getEnum(std::vector<t_inpfile>* inp, const char* name, warninp_t wi)
1818 // If we there's no valid option, we'll use the first enum entry as default.
1819 // Note, this assumes the enum is zero based, which is also assumed by
1820 // EnumerationWrapper and EnumerationArray.
1821 const auto defaultEnumValue = EnumType::Default;
1822 const auto& defaultName = enumValueToString(defaultEnumValue);
1823 // Get index of option in input
1824 const auto ii = get_einp(inp, name);
1827 // If the option wasn't set, we use the first enum entry as default
1828 inp->back().value_.assign(defaultName);
1829 return defaultEnumValue;
1832 // Check if option string can be mapped to a valid enum value
1833 const auto* optionString = (*inp)[ii].value_.c_str();
1834 for (auto enumValue : gmx::EnumerationWrapper<EnumType>{})
1836 if (gmx_strcasecmp_min(enumValueToString(enumValue), optionString) == 0)
1842 // If we get here, the option set is invalid. Print error.
1843 std::string errorMessage = gmx::formatString(
1844 "Invalid enum '%s' for variable %s, using '%s'\n", optionString, name, defaultName);
1845 errorMessage += gmx::formatString("Next time, use one of:");
1846 for (auto enumValue : gmx::EnumerationWrapper<EnumType>{})
1848 errorMessage += gmx::formatString(" '%s'", enumValueToString(enumValue));
1852 warning_error(wi, errorMessage);
1856 fprintf(stderr, "%s\n", errorMessage.c_str());
1858 (*inp)[ii].value_.assign(defaultName);
1859 return defaultEnumValue;
1862 /*! \brief Return whether an end state with the given coupling-lambda
1863 * value describes fully-interacting VDW.
1865 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1866 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1868 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1870 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1876 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1879 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1881 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1883 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1886 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1887 std::string message = gmx::formatExceptionMessageToString(*ex);
1888 warning_error(wi_, message.c_str());
1893 std::string getOptionName(const gmx::KeyValueTreePath& context)
1895 if (mapping_ != nullptr)
1897 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1898 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1901 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1906 const gmx::IKeyValueTreeBackMapping* mapping_;
1911 void get_ir(const char* mdparin,
1912 const char* mdparout,
1913 gmx::MDModules* mdModules,
1916 WriteMdpHeader writeMdpHeader,
1920 double dumdub[2][6];
1922 char warn_buf[STRLEN];
1923 t_lambda* fep = ir->fepvals;
1924 t_expanded* expand = ir->expandedvals;
1926 const char* no_names[] = { "no", nullptr };
1928 init_inputrec_strings();
1929 gmx::TextInputFile stream(mdparin);
1930 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1932 snew(dumstr[0], STRLEN);
1933 snew(dumstr[1], STRLEN);
1935 /* ignore the following deprecated commands */
1936 replace_inp_entry(inp, "title", nullptr);
1937 replace_inp_entry(inp, "cpp", nullptr);
1938 replace_inp_entry(inp, "domain-decomposition", nullptr);
1939 replace_inp_entry(inp, "andersen-seed", nullptr);
1940 replace_inp_entry(inp, "dihre", nullptr);
1941 replace_inp_entry(inp, "dihre-fc", nullptr);
1942 replace_inp_entry(inp, "dihre-tau", nullptr);
1943 replace_inp_entry(inp, "nstdihreout", nullptr);
1944 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1945 replace_inp_entry(inp, "optimize-fft", nullptr);
1946 replace_inp_entry(inp, "adress_type", nullptr);
1947 replace_inp_entry(inp, "adress_const_wf", nullptr);
1948 replace_inp_entry(inp, "adress_ex_width", nullptr);
1949 replace_inp_entry(inp, "adress_hy_width", nullptr);
1950 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1951 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1952 replace_inp_entry(inp, "adress_site", nullptr);
1953 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1954 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1955 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1956 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1957 replace_inp_entry(inp, "rlistlong", nullptr);
1958 replace_inp_entry(inp, "nstcalclr", nullptr);
1959 replace_inp_entry(inp, "pull-print-com2", nullptr);
1960 replace_inp_entry(inp, "gb-algorithm", nullptr);
1961 replace_inp_entry(inp, "nstgbradii", nullptr);
1962 replace_inp_entry(inp, "rgbradii", nullptr);
1963 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1964 replace_inp_entry(inp, "gb-saltconc", nullptr);
1965 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1966 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1967 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1968 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1969 replace_inp_entry(inp, "sa-algorithm", nullptr);
1970 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1971 replace_inp_entry(inp, "ns-type", nullptr);
1973 /* replace the following commands with the clearer new versions*/
1974 replace_inp_entry(inp, "unconstrained-start", "continuation");
1975 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1976 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1977 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1978 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1979 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1980 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1982 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1983 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1984 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1985 setStringEntry(&inp, "include", opts->include, nullptr);
1986 printStringNoNewline(
1987 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1988 setStringEntry(&inp, "define", opts->define, nullptr);
1990 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1991 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1992 printStringNoNewline(&inp, "Start time and timestep in ps");
1993 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1994 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1995 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1996 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1997 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1998 printStringNoNewline(
1999 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
2000 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
2001 printStringNoNewline(&inp, "Multiple time-stepping");
2002 ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
2005 gmx::GromppMtsOpts& mtsOpts = opts->mtsOpts;
2006 mtsOpts.numLevels = get_eint(&inp, "mts-levels", 2, wi);
2007 mtsOpts.level2Forces = setStringEntry(&inp, "mts-level2-forces", "longrange-nonbonded");
2008 mtsOpts.level2Factor = get_eint(&inp, "mts-level2-factor", 2, wi);
2010 // We clear after reading without dynamics to not force the user to remove MTS mdp options
2011 if (!EI_DYNAMICS(ir->eI))
2016 printStringNoNewline(&inp, "mode for center of mass motion removal");
2017 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
2018 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
2019 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
2020 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
2021 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
2023 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2024 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2025 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2026 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2029 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2030 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2031 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2032 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2033 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2034 ir->niter = get_eint(&inp, "niter", 20, wi);
2035 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2036 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2037 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2038 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2039 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2041 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2042 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2044 /* Output options */
2045 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2046 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2047 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2048 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2049 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2050 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2051 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2052 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2053 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2054 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2055 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2056 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2057 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2058 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2059 printStringNoNewline(&inp, "default, all atoms will be written.");
2060 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2061 printStringNoNewline(&inp, "Selection of energy groups");
2062 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2064 /* Neighbor searching */
2065 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2066 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2067 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
2068 printStringNoNewline(&inp, "nblist update frequency");
2069 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2070 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2071 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2072 std::vector<const char*> pbcTypesNamesChar;
2073 for (const auto& pbcTypeName : c_pbcTypeNames)
2075 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2077 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2078 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
2079 printStringNoNewline(&inp,
2080 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2081 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2082 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2083 printStringNoNewline(&inp, "nblist cut-off");
2084 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2085 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2087 /* Electrostatics */
2088 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2089 printStringNoNewline(&inp, "Method for doing electrostatics");
2090 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
2091 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
2092 printStringNoNewline(&inp, "cut-off lengths");
2093 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2094 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2095 printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
2096 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2097 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2098 printStringNoNewline(&inp, "Method for doing Van der Waals");
2099 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
2100 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
2101 printStringNoNewline(&inp, "cut-off lengths");
2102 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2103 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2104 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2105 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
2106 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2107 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2108 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2109 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2110 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2111 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2112 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2113 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2114 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2115 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2116 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2117 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2118 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2119 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2120 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
2121 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
2122 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2124 /* Implicit solvation is no longer supported, but we need grompp
2125 to be able to refuse old .mdp files that would have built a tpr
2126 to run it. Thus, only "no" is accepted. */
2127 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2129 /* Coupling stuff */
2130 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2131 printStringNoNewline(&inp, "Temperature coupling");
2132 ir->etc = getEnum<TemperatureCoupling>(&inp, "tcoupl", wi);
2133 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2134 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2135 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
2136 printStringNoNewline(&inp, "Groups to couple separately");
2137 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2138 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2139 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2140 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2141 printStringNoNewline(&inp, "pressure coupling");
2142 ir->epc = getEnum<PressureCoupling>(&inp, "pcoupl", wi);
2143 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
2144 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2145 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2146 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2147 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2148 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2149 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2150 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
2153 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2154 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
2155 printStringNoNewline(&inp, "Groups treated with MiMiC");
2156 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2158 /* Simulated annealing */
2159 printStringNewline(&inp, "SIMULATED ANNEALING");
2160 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2161 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2162 printStringNoNewline(&inp,
2163 "Number of time points to use for specifying annealing in each group");
2164 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2165 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2166 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2167 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2168 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2171 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2172 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
2173 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2174 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2177 printStringNewline(&inp, "OPTIONS FOR BONDS");
2178 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2179 printStringNoNewline(&inp, "Type of constraint algorithm");
2180 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
2181 printStringNoNewline(&inp, "Do not constrain the start configuration");
2182 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
2183 printStringNoNewline(&inp,
2184 "Use successive overrelaxation to reduce the number of shake iterations");
2185 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
2186 printStringNoNewline(&inp, "Relative tolerance of shake");
2187 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2188 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2189 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2190 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2191 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2192 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2193 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2194 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2195 printStringNoNewline(&inp, "rotates over more degrees than");
2196 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2197 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2198 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
2200 /* Energy group exclusions */
2201 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2202 printStringNoNewline(
2203 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2204 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2207 printStringNewline(&inp, "WALLS");
2208 printStringNoNewline(
2209 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2210 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2211 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
2212 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2213 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2214 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2215 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2218 printStringNewline(&inp, "COM PULLING");
2219 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2222 ir->pull = std::make_unique<pull_params_t>();
2223 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
2227 for (int c = 0; c < ir->pull->ncoord; c++)
2229 if (ir->pull->coord[c].eType == epullCONSTRAINT)
2232 "Constraint COM pulling is not supported in combination with "
2233 "multiple time stepping");
2241 NOTE: needs COM pulling or free energy input */
2242 printStringNewline(&inp, "AWH biasing");
2243 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2246 ir->awhParams = gmx::readAwhParams(&inp, wi);
2249 /* Enforced rotation */
2250 printStringNewline(&inp, "ENFORCED ROTATION");
2251 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2252 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2256 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2259 /* Interactive MD */
2261 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2262 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2263 if (inputrecStrings->imd_grp[0] != '\0')
2270 printStringNewline(&inp, "NMR refinement stuff");
2271 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2272 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2273 printStringNoNewline(
2274 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2275 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2276 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2277 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2278 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2279 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2280 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2281 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2282 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2283 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2284 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2285 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2286 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2287 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2288 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2289 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2291 /* free energy variables */
2292 printStringNewline(&inp, "Free energy variables");
2293 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2294 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2295 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2296 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2297 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2299 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2301 it was not entered */
2302 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2303 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2304 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2305 setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
2306 setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
2307 setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
2308 setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
2309 setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
2310 setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
2311 setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
2312 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2313 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2314 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2315 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2316 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2317 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2318 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2319 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2320 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2321 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2322 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2323 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2324 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2325 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2327 /* Non-equilibrium MD stuff */
2328 printStringNewline(&inp, "Non-equilibrium MD stuff");
2329 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2330 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2331 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2332 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2334 /* simulated tempering variables */
2335 printStringNewline(&inp, "simulated tempering variables");
2336 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2337 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2338 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2339 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2341 /* expanded ensemble variables */
2342 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2344 read_expandedparams(&inp, expand, wi);
2347 /* Electric fields */
2349 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2350 gmx::KeyValueTreeTransformer transform;
2351 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2352 mdModules->initMdpTransform(transform.rules());
2353 for (const auto& path : transform.mappedPaths())
2355 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2356 mark_einp_set(inp, path[0].c_str());
2358 MdpErrorHandler errorHandler(wi);
2359 auto result = transform.transform(convertedValues, &errorHandler);
2360 ir->params = new gmx::KeyValueTreeObject(result.object());
2361 mdModules->adjustInputrecBasedOnModules(ir);
2362 errorHandler.setBackMapping(result.backMapping());
2363 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2366 /* Ion/water position swapping ("computational electrophysiology") */
2367 printStringNewline(&inp,
2368 "Ion/water position swapping for computational electrophysiology setups");
2369 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2370 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2371 if (ir->eSwapCoords != eswapNO)
2378 printStringNoNewline(&inp, "Swap attempt frequency");
2379 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2380 printStringNoNewline(&inp, "Number of ion types to be controlled");
2381 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2384 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2386 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2387 snew(ir->swap->grp, ir->swap->ngrp);
2388 for (i = 0; i < ir->swap->ngrp; i++)
2390 snew(ir->swap->grp[i].molname, STRLEN);
2392 printStringNoNewline(&inp,
2393 "Two index groups that contain the compartment-partitioning atoms");
2394 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2395 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2396 printStringNoNewline(&inp,
2397 "Use center of mass of split groups (yes/no), otherwise center of "
2398 "geometry is used");
2399 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2400 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2402 printStringNoNewline(&inp, "Name of solvent molecules");
2403 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2405 printStringNoNewline(&inp,
2406 "Split cylinder: radius, upper and lower extension (nm) (this will "
2407 "define the channels)");
2408 printStringNoNewline(&inp,
2409 "Note that the split cylinder settings do not have an influence on "
2410 "the swapping protocol,");
2411 printStringNoNewline(
2413 "however, if correctly defined, the permeation events are recorded per channel");
2414 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2415 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2416 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2417 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2418 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2419 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2421 printStringNoNewline(
2423 "Average the number of ions per compartment over these many swap attempt steps");
2424 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2426 printStringNoNewline(
2427 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2428 printStringNoNewline(
2429 &inp, "and the requested number of ions of this type in compartments A and B");
2430 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2431 for (i = 0; i < nIonTypes; i++)
2433 int ig = eSwapFixedGrpNR + i;
2435 sprintf(buf, "iontype%d-name", i);
2436 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2437 sprintf(buf, "iontype%d-in-A", i);
2438 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2439 sprintf(buf, "iontype%d-in-B", i);
2440 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2443 printStringNoNewline(
2445 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2446 printStringNoNewline(
2448 "at maximum distance (= bulk concentration) to the split group layers. However,");
2449 printStringNoNewline(&inp,
2450 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2451 "layer from the middle at 0.0");
2452 printStringNoNewline(&inp,
2453 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2454 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2455 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2456 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2457 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2459 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2462 printStringNoNewline(
2463 &inp, "Start to swap ions if threshold difference to requested count is reached");
2464 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2467 /* AdResS is no longer supported, but we need grompp to be able to
2468 refuse to process old .mdp files that used it. */
2469 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2471 /* User defined thingies */
2472 printStringNewline(&inp, "User defined thingies");
2473 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2474 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2475 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2476 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2477 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2478 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2479 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2480 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2481 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2482 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2486 gmx::TextOutputFile stream(mdparout);
2487 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2489 // Transform module data into a flat key-value tree for output.
2490 gmx::KeyValueTreeBuilder builder;
2491 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2492 mdModules->buildMdpOutput(&builderObject);
2494 gmx::TextWriter writer(&stream);
2495 writeKeyValueTreeAsMdp(&writer, builder.build());
2500 /* Process options if necessary */
2501 for (m = 0; m < 2; m++)
2503 for (i = 0; i < 2 * DIM; i++)
2507 if (ir->epc != PressureCoupling::No)
2512 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2516 "Pressure coupling incorrect number of values (I need exactly 1)");
2518 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2520 case epctSEMIISOTROPIC:
2521 case epctSURFACETENSION:
2522 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2526 "Pressure coupling incorrect number of values (I need exactly 2)");
2528 dumdub[m][YY] = dumdub[m][XX];
2530 case epctANISOTROPIC:
2531 if (sscanf(dumstr[m],
2532 "%lf%lf%lf%lf%lf%lf",
2543 "Pressure coupling incorrect number of values (I need exactly 6)");
2548 "Pressure coupling type %s not implemented yet",
2549 epcoupltype_names[ir->epct]);
2553 clear_mat(ir->ref_p);
2554 clear_mat(ir->compress);
2555 for (i = 0; i < DIM; i++)
2557 ir->ref_p[i][i] = dumdub[1][i];
2558 ir->compress[i][i] = dumdub[0][i];
2560 if (ir->epct == epctANISOTROPIC)
2562 ir->ref_p[XX][YY] = dumdub[1][3];
2563 ir->ref_p[XX][ZZ] = dumdub[1][4];
2564 ir->ref_p[YY][ZZ] = dumdub[1][5];
2565 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2568 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2569 "apply a threefold shear stress?\n");
2571 ir->compress[XX][YY] = dumdub[0][3];
2572 ir->compress[XX][ZZ] = dumdub[0][4];
2573 ir->compress[YY][ZZ] = dumdub[0][5];
2574 for (i = 0; i < DIM; i++)
2576 for (m = 0; m < i; m++)
2578 ir->ref_p[i][m] = ir->ref_p[m][i];
2579 ir->compress[i][m] = ir->compress[m][i];
2584 if (ir->comm_mode == ecmNO)
2589 opts->couple_moltype = nullptr;
2590 if (strlen(inputrecStrings->couple_moltype) > 0)
2592 if (ir->efep != efepNO)
2594 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2595 if (opts->couple_lam0 == opts->couple_lam1)
2597 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2599 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2603 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2610 "Free energy is turned off, so we will not decouple the molecule listed "
2614 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2615 if (ir->efep != efepNO)
2617 if (fep->delta_lambda != 0)
2619 ir->efep = efepSLOWGROWTH;
2623 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2625 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2627 "Old option for dhdl-print-energy given: "
2628 "changing \"yes\" to \"total\"\n");
2631 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2633 /* always print out the energy to dhdl if we are doing
2634 expanded ensemble, since we need the total energy for
2635 analysis if the temperature is changing. In some
2636 conditions one may only want the potential energy, so
2637 we will allow that if the appropriate mdp setting has
2638 been enabled. Otherwise, total it is:
2640 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2643 if ((ir->efep != efepNO) || ir->bSimTemp)
2645 ir->bExpanded = FALSE;
2646 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2648 ir->bExpanded = TRUE;
2650 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2651 if (ir->bSimTemp) /* done after fep params */
2653 do_simtemp_params(ir);
2656 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2657 * setup and not on the old way of specifying the free-energy setup,
2658 * we should check for using soft-core when not needed, since that
2659 * can complicate the sampling significantly.
2660 * Note that we only check for the automated coupling setup.
2661 * If the (advanced) user does FEP through manual topology changes,
2662 * this check will not be triggered.
2664 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 && ir->fepvals->sc_alpha != 0
2665 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2668 "You are using soft-core interactions while the Van der Waals interactions are "
2669 "not decoupled (note that the sc-coul option is only active when using lambda "
2670 "states). Although this will not lead to errors, you will need much more "
2671 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2676 ir->fepvals->n_lambda = 0;
2679 /* WALL PARAMETERS */
2681 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2683 /* ORIENTATION RESTRAINT PARAMETERS */
2685 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2687 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2690 /* DEFORMATION PARAMETERS */
2692 clear_mat(ir->deform);
2693 for (i = 0; i < 6; i++)
2698 double gmx_unused canary;
2699 int ndeform = sscanf(inputrecStrings->deform,
2700 "%lf %lf %lf %lf %lf %lf %lf",
2709 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2713 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2714 inputrecStrings->deform)
2717 for (i = 0; i < 3; i++)
2719 ir->deform[i][i] = dumdub[0][i];
2721 ir->deform[YY][XX] = dumdub[0][3];
2722 ir->deform[ZZ][XX] = dumdub[0][4];
2723 ir->deform[ZZ][YY] = dumdub[0][5];
2724 if (ir->epc != PressureCoupling::No)
2726 for (i = 0; i < 3; i++)
2728 for (j = 0; j <= i; j++)
2730 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2732 warning_error(wi, "A box element has deform set and compressibility > 0");
2736 for (i = 0; i < 3; i++)
2738 for (j = 0; j < i; j++)
2740 if (ir->deform[i][j] != 0)
2742 for (m = j; m < DIM; m++)
2744 if (ir->compress[m][j] != 0)
2747 "An off-diagonal box element has deform set while "
2748 "compressibility > 0 for the same component of another box "
2749 "vector, this might lead to spurious periodicity effects.");
2750 warning(wi, warn_buf);
2758 /* Ion/water position swapping checks */
2759 if (ir->eSwapCoords != eswapNO)
2761 if (ir->swap->nstswap < 1)
2763 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2765 if (ir->swap->nAverage < 1)
2767 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2769 if (ir->swap->threshold < 1.0)
2771 warning_error(wi, "Ion count threshold must be at least 1.\n");
2775 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2778 std::vector<std::string> errorMessages;
2779 ir->mtsLevels = gmx::setupMtsLevels(opts->mtsOpts, &errorMessages);
2781 for (const auto& errorMessage : errorMessages)
2783 warning_error(wi, errorMessage.c_str());
2789 gmx::checkAwhParams(ir->awhParams, ir, wi);
2796 /* We would like gn to be const as well, but C doesn't allow this */
2797 /* TODO this is utility functionality (search for the index of a
2798 string in a collection), so should be refactored and located more
2800 int search_string(const char* s, int ng, char* gn[])
2804 for (i = 0; (i < ng); i++)
2806 if (gmx_strcasecmp(s, gn[i]) == 0)
2813 "Group %s referenced in the .mdp file was not found in the index file.\n"
2814 "Group names must match either [moleculetype] names or custom index group\n"
2815 "names, in which case you must supply an index file to the '-n' option\n"
2820 static void atomGroupRangeValidation(int natoms, int groupIndex, const t_blocka& block)
2822 /* Now go over the atoms in the group */
2823 for (int j = block.index[groupIndex]; (j < block.index[groupIndex + 1]); j++)
2825 int aj = block.a[j];
2827 /* Range checking */
2828 if ((aj < 0) || (aj >= natoms))
2830 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2835 static void do_numbering(int natoms,
2836 SimulationGroups* groups,
2837 gmx::ArrayRef<std::string> groupsFromMdpFile,
2840 SimulationAtomGroupType gtype,
2846 unsigned short* cbuf;
2847 AtomGroupIndices* grps = &(groups->groups[gtype]);
2850 char warn_buf[STRLEN];
2852 title = shortName(gtype);
2855 /* Mark all id's as not set */
2856 for (int i = 0; (i < natoms); i++)
2861 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2863 /* Lookup the group name in the block structure */
2864 const int gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2865 if ((grptp != egrptpONE) || (i == 0))
2867 grps->emplace_back(gid);
2869 GMX_ASSERT(block, "Can't have a nullptr block");
2870 atomGroupRangeValidation(natoms, gid, *block);
2871 /* Now go over the atoms in the group */
2872 for (int j = block->index[gid]; (j < block->index[gid + 1]); j++)
2874 const int aj = block->a[j];
2875 /* Lookup up the old group number */
2876 const int ognr = cbuf[aj];
2879 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title, ognr + 1, i + 1);
2883 /* Store the group number in buffer */
2884 if (grptp == egrptpONE)
2897 /* Now check whether we have done all atoms */
2900 if (grptp == egrptpALL)
2902 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2904 else if (grptp == egrptpPART)
2906 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2907 warning_note(wi, warn_buf);
2909 /* Assign all atoms currently unassigned to a rest group */
2910 for (int j = 0; (j < natoms); j++)
2912 if (cbuf[j] == NOGID)
2914 cbuf[j] = grps->size();
2917 if (grptp != egrptpPART)
2921 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title, natoms - ntot);
2923 /* Add group name "rest" */
2924 grps->emplace_back(restnm);
2926 /* Assign the rest name to all atoms not currently assigned to a group */
2927 for (int j = 0; (j < natoms); j++)
2929 if (cbuf[j] == NOGID)
2931 // group size was not updated before this here, so need to use -1.
2932 cbuf[j] = grps->size() - 1;
2938 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2940 /* All atoms are part of one (or no) group, no index required */
2941 groups->groupNumbers[gtype].clear();
2945 for (int j = 0; (j < natoms); j++)
2947 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2954 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2957 pull_params_t* pull;
2958 int natoms, imin, jmin;
2959 int * nrdf2, *na_vcm, na_tot;
2960 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2965 * First calc 3xnr-atoms for each group
2966 * then subtract half a degree of freedom for each constraint
2968 * Only atoms and nuclei contribute to the degrees of freedom...
2973 const SimulationGroups& groups = mtop->groups;
2974 natoms = mtop->natoms;
2976 /* Allocate one more for a possible rest group */
2977 /* We need to sum degrees of freedom into doubles,
2978 * since floats give too low nrdf's above 3 million atoms.
2980 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
2981 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2982 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2983 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2984 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2986 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2990 for (gmx::index i = 0;
2991 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
2995 clear_ivec(dof_vcm[i]);
2997 nrdf_vcm_sub[i] = 0;
2999 snew(nrdf2, natoms);
3000 for (const AtomProxy atomP : AtomRange(*mtop))
3002 const t_atom& local = atomP.atom();
3003 int i = atomP.globalAtomNumber();
3005 if (local.ptype == eptAtom || local.ptype == eptNucleus)
3007 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
3008 for (int d = 0; d < DIM; d++)
3010 if (opts->nFreeze[g][d] == 0)
3012 /* Add one DOF for particle i (counted as 2*1) */
3014 /* VCM group i has dim d as a DOF */
3015 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
3019 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
3021 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
3027 for (const gmx_molblock_t& molb : mtop->molblock)
3029 const gmx_moltype_t& molt = mtop->moltype[molb.type];
3030 const t_atom* atom = molt.atoms.atom;
3031 for (int mol = 0; mol < molb.nmol; mol++)
3033 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
3035 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
3036 for (int i = 0; i < molt.ilist[ftype].size();)
3038 /* Subtract degrees of freedom for the constraints,
3039 * if the particles still have degrees of freedom left.
3040 * If one of the particles is a vsite or a shell, then all
3041 * constraint motion will go there, but since they do not
3042 * contribute to the constraints the degrees of freedom do not
3045 int ai = as + ia[i + 1];
3046 int aj = as + ia[i + 2];
3047 if (((atom[ia[i + 1]].ptype == eptNucleus) || (atom[ia[i + 1]].ptype == eptAtom))
3048 && ((atom[ia[i + 2]].ptype == eptNucleus) || (atom[ia[i + 2]].ptype == eptAtom)))
3066 imin = std::min(imin, nrdf2[ai]);
3067 jmin = std::min(jmin, nrdf2[aj]);
3070 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3072 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3074 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3076 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3079 i += interaction_function[ftype].nratoms + 1;
3082 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3083 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3085 /* Subtract 1 dof from every atom in the SETTLE */
3086 for (int j = 0; j < 3; j++)
3088 int ai = as + ia[i + 1 + j];
3089 imin = std::min(2, nrdf2[ai]);
3091 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3093 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3098 as += molt.atoms.nr;
3104 /* Correct nrdf for the COM constraints.
3105 * We correct using the TC and VCM group of the first atom
3106 * in the reference and pull group. If atoms in one pull group
3107 * belong to different TC or VCM groups it is anyhow difficult
3108 * to determine the optimal nrdf assignment.
3110 pull = ir->pull.get();
3112 for (int i = 0; i < pull->ncoord; i++)
3114 if (pull->coord[i].eType != epullCONSTRAINT)
3121 for (int j = 0; j < 2; j++)
3123 const t_pull_group* pgrp;
3125 pgrp = &pull->group[pull->coord[i].group[j]];
3127 if (!pgrp->ind.empty())
3129 /* Subtract 1/2 dof from each group */
3130 int ai = pgrp->ind[0];
3131 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3133 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3135 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3138 "Center of mass pulling constraints caused the number of degrees "
3139 "of freedom for temperature coupling group %s to be negative",
3140 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3141 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3146 /* We need to subtract the whole DOF from group j=1 */
3153 if (ir->nstcomm != 0)
3155 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3156 "Expect at least one group when removing COM motion");
3158 /* We remove COM motion up to dim ndof_com() */
3159 const int ndim_rm_vcm = ndof_com(ir);
3161 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3162 * the number of degrees of freedom in each vcm group when COM
3163 * translation is removed and 6 when rotation is removed as well.
3164 * Note that we do not and should not include the rest group here.
3166 for (gmx::index j = 0;
3167 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]);
3170 switch (ir->comm_mode)
3173 case ecmLINEAR_ACCELERATION_CORRECTION:
3174 nrdf_vcm_sub[j] = 0;
3175 for (int d = 0; d < ndim_rm_vcm; d++)
3183 case ecmANGULAR: nrdf_vcm_sub[j] = 6; break;
3184 default: gmx_incons("Checking comm_mode");
3188 for (gmx::index i = 0;
3189 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]);
3192 /* Count the number of atoms of TC group i for every VCM group */
3193 for (gmx::index j = 0;
3194 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3200 for (int ai = 0; ai < natoms; ai++)
3202 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3204 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3208 /* Correct for VCM removal according to the fraction of each VCM
3209 * group present in this TC group.
3211 nrdf_uc = nrdf_tc[i];
3213 for (gmx::index j = 0;
3214 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1;
3217 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3219 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3220 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3225 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3227 opts->nrdf[i] = nrdf_tc[i];
3228 if (opts->nrdf[i] < 0)
3233 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3234 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]],
3243 sfree(nrdf_vcm_sub);
3246 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3248 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3249 * But since this is much larger than STRLEN, such a line can not be parsed.
3250 * The real maximum is the number of names that fit in a string: STRLEN/2.
3252 #define EGP_MAX (STRLEN / 2)
3256 auto names = gmx::splitString(val);
3257 if (names.size() % 2 != 0)
3259 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3261 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3263 for (size_t i = 0; i < names.size() / 2; i++)
3265 // TODO this needs to be replaced by a solution using std::find_if
3269 names[2 * i].c_str(),
3270 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3276 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3281 names[2 * i + 1].c_str(),
3282 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3288 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3290 if ((j < nr) && (k < nr))
3292 ir->opts.egp_flags[nr * j + k] |= flag;
3293 ir->opts.egp_flags[nr * k + j] |= flag;
3302 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3304 int ig = -1, i = 0, gind;
3308 /* Just a quick check here, more thorough checks are in mdrun */
3309 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3311 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3314 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3315 for (ig = 0; ig < swap->ngrp; ig++)
3317 swapg = &swap->grp[ig];
3318 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3319 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3324 "%s group '%s' contains %d atoms.\n",
3325 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3326 swap->grp[ig].molname,
3328 snew(swapg->ind, swapg->nat);
3329 for (i = 0; i < swapg->nat; i++)
3331 swapg->ind[i] = grps->a[grps->index[gind] + i];
3336 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3342 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3347 ig = search_string(IMDgname, grps->nr, gnames);
3348 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3350 if (IMDgroup->nat > 0)
3353 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3357 snew(IMDgroup->ind, IMDgroup->nat);
3358 for (i = 0; i < IMDgroup->nat; i++)
3360 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3365 /* Checks whether atoms are both part of a COM removal group and frozen.
3366 * If a fully frozen atom is part of a COM removal group, it is removed
3367 * from the COM removal group. A note is issued if such atoms are present.
3368 * A warning is issued for atom with one or two dimensions frozen that
3369 * are part of a COM removal group (mdrun would need to compute COM mass
3370 * per dimension to handle this correctly).
3371 * Also issues a warning when non-frozen atoms are not part of a COM
3372 * removal group while COM removal is active.
3374 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3376 const t_grpopts& opts,
3379 const int vcmRestGroup =
3380 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3382 int numFullyFrozenVcmAtoms = 0;
3383 int numPartiallyFrozenVcmAtoms = 0;
3384 int numNonVcmAtoms = 0;
3385 for (int a = 0; a < numAtoms; a++)
3387 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3388 int numFrozenDims = 0;
3389 for (int d = 0; d < DIM; d++)
3391 numFrozenDims += opts.nFreeze[freezeGroup][d];
3394 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3395 if (vcmGroup < vcmRestGroup)
3397 if (numFrozenDims == DIM)
3399 /* Do not remove COM motion for this fully frozen atom */
3400 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3402 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3405 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3406 numFullyFrozenVcmAtoms++;
3408 else if (numFrozenDims > 0)
3410 numPartiallyFrozenVcmAtoms++;
3413 else if (numFrozenDims < DIM)
3419 if (numFullyFrozenVcmAtoms > 0)
3421 std::string warningText = gmx::formatString(
3422 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3423 "removing these atoms from the COMM removal group(s)",
3424 numFullyFrozenVcmAtoms);
3425 warning_note(wi, warningText.c_str());
3427 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3429 std::string warningText = gmx::formatString(
3430 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3431 "removal group(s), due to limitations in the code these still contribute to the "
3432 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3434 numPartiallyFrozenVcmAtoms,
3436 warning(wi, warningText.c_str());
3438 if (numNonVcmAtoms > 0)
3440 std::string warningText = gmx::formatString(
3441 "%d atoms are not part of any center of mass motion removal group.\n"
3442 "This may lead to artifacts.\n"
3443 "In most cases one should use one group for the whole system.",
3445 warning(wi, warningText.c_str());
3449 void do_index(const char* mdparin,
3453 const gmx::MdModulesNotifier& notifier,
3457 t_blocka* defaultIndexGroups;
3465 int i, j, k, restnm;
3466 bool bExcl, bTable, bAnneal;
3467 char warn_buf[STRLEN];
3471 fprintf(stderr, "processing index file...\n");
3475 snew(defaultIndexGroups, 1);
3476 snew(defaultIndexGroups->index, 1);
3478 atoms_all = gmx_mtop_global_atoms(mtop);
3479 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3480 done_atom(&atoms_all);
3484 defaultIndexGroups = init_index(ndx, &gnames);
3487 SimulationGroups* groups = &mtop->groups;
3488 natoms = mtop->natoms;
3489 symtab = &mtop->symtab;
3491 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3493 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3495 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3496 restnm = groups->groupNames.size() - 1;
3497 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3498 srenew(gnames, defaultIndexGroups->nr + 1);
3499 gnames[restnm] = *(groups->groupNames.back());
3501 set_warning_line(wi, mdparin, -1);
3503 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3504 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3505 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3506 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3507 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3510 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3512 temperatureCouplingGroupNames.size(),
3513 temperatureCouplingReferenceValues.size(),
3514 temperatureCouplingTauValues.size());
3517 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3518 do_numbering(natoms,
3520 temperatureCouplingGroupNames,
3523 SimulationAtomGroupType::TemperatureCoupling,
3525 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST,
3528 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3530 snew(ir->opts.nrdf, nr);
3531 snew(ir->opts.tau_t, nr);
3532 snew(ir->opts.ref_t, nr);
3533 if (ir->eI == eiBD && ir->bd_fric == 0)
3535 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3538 if (useReferenceTemperature)
3540 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3542 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3546 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3547 for (i = 0; (i < nr); i++)
3549 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3551 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3552 warning_error(wi, warn_buf);
3555 if (ir->etc != TemperatureCoupling::VRescale && ir->opts.tau_t[i] == 0)
3559 "tau-t = -1 is the value to signal that a group should not have "
3560 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3563 if (ir->opts.tau_t[i] >= 0)
3565 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3568 if (ir->etc != TemperatureCoupling::No && ir->nsttcouple == -1)
3570 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3575 if ((ir->etc == TemperatureCoupling::NoseHoover) && (ir->epc == PressureCoupling::Berendsen))
3578 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3579 "md-vv; use either vrescale temperature with berendsen pressure or "
3580 "Nose-Hoover temperature with MTTK pressure");
3582 if (ir->epc == PressureCoupling::Mttk)
3584 if (ir->etc != TemperatureCoupling::NoseHoover)
3587 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3592 if (ir->nstpcouple != ir->nsttcouple)
3594 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3595 ir->nstpcouple = ir->nsttcouple = mincouple;
3597 "for current Trotter decomposition methods with vv, nsttcouple and "
3598 "nstpcouple must be equal. Both have been reset to "
3599 "min(nsttcouple,nstpcouple) = %d",
3601 warning_note(wi, warn_buf);
3606 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3607 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3609 if (ETC_ANDERSEN(ir->etc))
3611 if (ir->nsttcouple != 1)
3615 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3616 "need for larger nsttcouple > 1, since no global parameters are computed. "
3617 "nsttcouple has been reset to 1");
3618 warning_note(wi, warn_buf);
3621 nstcmin = tcouple_min_integration_steps(ir->etc);
3624 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3627 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3628 "least %d times larger than nsttcouple*dt (%g)",
3629 enumValueToString(ir->etc),
3632 ir->nsttcouple * ir->delta_t);
3633 warning(wi, warn_buf);
3636 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3637 for (i = 0; (i < nr); i++)
3639 if (ir->opts.ref_t[i] < 0)
3641 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3644 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3645 if we are in this conditional) if mc_temp is negative */
3646 if (ir->expandedvals->mc_temp < 0)
3648 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3652 /* Simulated annealing for each group. There are nr groups */
3653 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3654 if (simulatedAnnealingGroupNames.size() == 1
3655 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3657 simulatedAnnealingGroupNames.resize(0);
3659 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3662 "Wrong number of annealing values: %zu (for %d groups)\n",
3663 simulatedAnnealingGroupNames.size(),
3668 snew(ir->opts.annealing, nr);
3669 snew(ir->opts.anneal_npoints, nr);
3670 snew(ir->opts.anneal_time, nr);
3671 snew(ir->opts.anneal_temp, nr);
3672 for (i = 0; i < nr; i++)
3674 ir->opts.annealing[i] = eannNO;
3675 ir->opts.anneal_npoints[i] = 0;
3676 ir->opts.anneal_time[i] = nullptr;
3677 ir->opts.anneal_temp[i] = nullptr;
3679 if (!simulatedAnnealingGroupNames.empty())
3682 for (i = 0; i < nr; i++)
3684 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3686 ir->opts.annealing[i] = eannNO;
3688 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3690 ir->opts.annealing[i] = eannSINGLE;
3693 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3695 ir->opts.annealing[i] = eannPERIODIC;
3701 /* Read the other fields too */
3702 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3703 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3706 "Found %zu annealing-npoints values for %zu groups\n",
3707 simulatedAnnealingPoints.size(),
3708 simulatedAnnealingGroupNames.size());
3710 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3711 size_t numSimulatedAnnealingFields = 0;
3712 for (i = 0; i < nr; i++)
3714 if (ir->opts.anneal_npoints[i] == 1)
3718 "Please specify at least a start and an end point for annealing\n");
3720 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3721 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3722 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3725 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3727 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3730 "Found %zu annealing-time values, wanted %zu\n",
3731 simulatedAnnealingTimes.size(),
3732 numSimulatedAnnealingFields);
3734 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3735 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3738 "Found %zu annealing-temp values, wanted %zu\n",
3739 simulatedAnnealingTemperatures.size(),
3740 numSimulatedAnnealingFields);
3743 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3744 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3745 convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
3747 simulatedAnnealingTemperatures,
3749 allSimulatedAnnealingTemperatures.data());
3750 for (i = 0, k = 0; i < nr; i++)
3752 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3754 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3755 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3758 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3760 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3766 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3769 "Annealing timepoints out of order: t=%f comes after "
3771 ir->opts.anneal_time[i][j],
3772 ir->opts.anneal_time[i][j - 1]);
3775 if (ir->opts.anneal_temp[i][j] < 0)
3778 "Found negative temperature in annealing: %f\n",
3779 ir->opts.anneal_temp[i][j]);
3784 /* Print out some summary information, to make sure we got it right */
3785 for (i = 0; i < nr; i++)
3787 if (ir->opts.annealing[i] != eannNO)
3789 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3791 "Simulated annealing for group %s: %s, %d timepoints\n",
3792 *(groups->groupNames[j]),
3793 eann_names[ir->opts.annealing[i]],
3794 ir->opts.anneal_npoints[i]);
3795 fprintf(stderr, "Time (ps) Temperature (K)\n");
3796 /* All terms except the last one */
3797 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3801 ir->opts.anneal_time[i][j],
3802 ir->opts.anneal_temp[i][j]);
3805 /* Finally the last one */
3806 j = ir->opts.anneal_npoints[i] - 1;
3807 if (ir->opts.annealing[i] == eannSINGLE)
3811 ir->opts.anneal_time[i][j],
3812 ir->opts.anneal_temp[i][j]);
3818 ir->opts.anneal_time[i][j],
3819 ir->opts.anneal_temp[i][j]);
3820 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3823 "There is a temperature jump when your annealing "
3835 for (int i = 1; i < ir->pull->ngroup; i++)
3837 const int gid = search_string(
3838 inputrecStrings->pullGroupNames[i].c_str(), defaultIndexGroups->nr, gnames);
3839 GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
3840 atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
3843 process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3845 checkPullCoords(ir->pull->group, ir->pull->coord);
3850 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3853 if (ir->eSwapCoords != eswapNO)
3855 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3858 /* Make indices for IMD session */
3861 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3864 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3865 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3866 notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
3868 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3869 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3870 if (freezeDims.size() != DIM * freezeGroupNames.size())
3873 "Invalid Freezing input: %zu groups and %zu freeze values",
3874 freezeGroupNames.size(),
3877 do_numbering(natoms,
3882 SimulationAtomGroupType::Freeze,
3887 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3888 ir->opts.ngfrz = nr;
3889 snew(ir->opts.nFreeze, nr);
3890 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3892 for (j = 0; (j < DIM); j++, k++)
3894 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3895 if (!ir->opts.nFreeze[i][j])
3897 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3900 "Please use Y(ES) or N(O) for freezedim only "
3902 freezeDims[k].c_str());
3903 warning(wi, warn_buf);
3908 for (; (i < nr); i++)
3910 for (j = 0; (j < DIM); j++)
3912 ir->opts.nFreeze[i][j] = 0;
3916 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3917 do_numbering(natoms,
3922 SimulationAtomGroupType::EnergyOutput,
3927 add_wall_energrps(groups, ir->nwall, symtab);
3928 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3929 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3930 do_numbering(natoms,
3935 SimulationAtomGroupType::MassCenterVelocityRemoval,
3937 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART,
3941 if (ir->comm_mode != ecmNO)
3943 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3946 /* Now we have filled the freeze struct, so we can calculate NRDF */
3947 calc_nrdf(mtop, ir, gnames);
3949 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3950 do_numbering(natoms,
3955 SimulationAtomGroupType::User1,
3960 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3961 do_numbering(natoms,
3966 SimulationAtomGroupType::User2,
3971 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3972 do_numbering(natoms,
3974 compressedXGroupNames,
3977 SimulationAtomGroupType::CompressedPositionOutput,
3982 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3983 do_numbering(natoms,
3985 orirefFitGroupNames,
3988 SimulationAtomGroupType::OrientationRestraintsFit,
3994 /* MiMiC QMMM input processing */
3995 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
3996 if (qmGroupNames.size() > 1)
3998 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
4000 /* group rest, if any, is always MM! */
4001 do_numbering(natoms,
4006 SimulationAtomGroupType::QuantumMechanics,
4011 ir->opts.ngQM = qmGroupNames.size();
4013 /* end of MiMiC QMMM input */
4017 for (auto group : gmx::keysOf(groups->groups))
4019 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
4020 for (const auto& entry : groups->groups[group])
4022 fprintf(stderr, " %s", *(groups->groupNames[entry]));
4024 fprintf(stderr, "\n");
4028 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
4029 snew(ir->opts.egp_flags, nr * nr);
4031 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
4032 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
4034 warning_error(wi, "Energy group exclusions are currently not supported");
4036 if (bExcl && EEL_FULL(ir->coulombtype))
4038 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
4041 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
4042 if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
4043 && !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
4046 "Can only have energy group pair tables in combination with user tables for VdW "
4050 /* final check before going out of scope if simulated tempering variables
4051 * need to be set to default values.
4053 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
4055 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
4058 "the value for nstexpanded was not specified for "
4059 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
4060 "by default, but it is recommended to set it to an explicit value!",
4061 ir->expandedvals->nstexpanded));
4063 for (i = 0; (i < defaultIndexGroups->nr); i++)
4068 done_blocka(defaultIndexGroups);
4069 sfree(defaultIndexGroups);
4073 static void check_disre(const gmx_mtop_t* mtop)
4075 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
4077 const gmx_ffparams_t& ffparams = mtop->ffparams;
4080 for (int i = 0; i < ffparams.numTypes(); i++)
4082 int ftype = ffparams.functype[i];
4083 if (ftype == F_DISRES)
4085 int label = ffparams.iparams[i].disres.label;
4086 if (label == old_label)
4088 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
4097 "Found %d double distance restraint indices,\n"
4098 "probably the parameters for multiple pairs in one restraint "
4099 "are not identical\n",
4105 static bool absolute_reference(const t_inputrec* ir, const gmx_mtop_t* sys, const bool posres_only, ivec AbsRef)
4108 gmx_mtop_ilistloop_t iloop;
4110 const t_iparams* pr;
4117 for (d = 0; d < DIM; d++)
4119 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
4121 /* Check for freeze groups */
4122 for (g = 0; g < ir->opts.ngfrz; g++)
4124 for (d = 0; d < DIM; d++)
4126 if (ir->opts.nFreeze[g][d] != 0)
4134 /* Check for position restraints */
4135 iloop = gmx_mtop_ilistloop_init(sys);
4136 while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
4138 if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
4140 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
4142 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
4143 for (d = 0; d < DIM; d++)
4145 if (pr->posres.fcA[d] != 0)
4151 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
4153 /* Check for flat-bottom posres */
4154 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
4155 if (pr->fbposres.k != 0)
4157 switch (pr->fbposres.geom)
4159 case efbposresSPHERE: AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1; break;
4160 case efbposresCYLINDERX: AbsRef[YY] = AbsRef[ZZ] = 1; break;
4161 case efbposresCYLINDERY: AbsRef[XX] = AbsRef[ZZ] = 1; break;
4162 case efbposresCYLINDER:
4163 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4164 case efbposresCYLINDERZ: AbsRef[XX] = AbsRef[YY] = 1; break;
4165 case efbposresX: /* d=XX */
4166 case efbposresY: /* d=YY */
4167 case efbposresZ: /* d=ZZ */
4168 d = pr->fbposres.geom - efbposresX;
4173 " Invalid geometry for flat-bottom position restraint.\n"
4174 "Expected nr between 1 and %d. Found %d\n",
4183 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
4186 static void check_combination_rule_differences(const gmx_mtop_t* mtop,
4188 bool* bC6ParametersWorkWithGeometricRules,
4189 bool* bC6ParametersWorkWithLBRules,
4190 bool* bLBRulesPossible)
4192 int ntypes, tpi, tpj;
4195 double c6i, c6j, c12i, c12j;
4196 double c6, c6_geometric, c6_LB;
4197 double sigmai, sigmaj, epsi, epsj;
4198 bool bCanDoLBRules, bCanDoGeometricRules;
4201 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4202 * force-field floating point parameters.
4205 ptr = getenv("GMX_LJCOMB_TOL");
4209 double gmx_unused canary;
4211 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4214 FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4219 *bC6ParametersWorkWithLBRules = TRUE;
4220 *bC6ParametersWorkWithGeometricRules = TRUE;
4221 bCanDoLBRules = TRUE;
4222 ntypes = mtop->ffparams.atnr;
4223 snew(typecount, ntypes);
4224 gmx_mtop_count_atomtypes(mtop, state, typecount);
4225 *bLBRulesPossible = TRUE;
4226 for (tpi = 0; tpi < ntypes; ++tpi)
4228 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4229 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4230 for (tpj = tpi; tpj < ntypes; ++tpj)
4232 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4233 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4234 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4235 c6_geometric = std::sqrt(c6i * c6j);
4236 if (!gmx_numzero(c6_geometric))
4238 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4240 sigmai = gmx::sixthroot(c12i / c6i);
4241 sigmaj = gmx::sixthroot(c12j / c6j);
4242 epsi = c6i * c6i / (4.0 * c12i);
4243 epsj = c6j * c6j / (4.0 * c12j);
4244 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4248 *bLBRulesPossible = FALSE;
4249 c6_LB = c6_geometric;
4251 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4256 *bC6ParametersWorkWithLBRules = FALSE;
4259 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4261 if (!bCanDoGeometricRules)
4263 *bC6ParametersWorkWithGeometricRules = FALSE;
4270 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop, warninp_t wi)
4272 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4274 check_combination_rule_differences(
4275 mtop, 0, &bC6ParametersWorkWithGeometricRules, &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4276 if (ir->ljpme_combination_rule == eljpmeLB)
4278 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4281 "You are using arithmetic-geometric combination rules "
4282 "in LJ-PME, but your non-bonded C6 parameters do not "
4283 "follow these rules.");
4288 if (!bC6ParametersWorkWithGeometricRules)
4290 if (ir->eDispCorr != edispcNO)
4293 "You are using geometric combination rules in "
4294 "LJ-PME, but your non-bonded C6 parameters do "
4295 "not follow these rules. "
4296 "This will introduce very small errors in the forces and energies in "
4297 "your simulations. Dispersion correction will correct total energy "
4298 "and/or pressure for isotropic systems, but not forces or surface "
4304 "You are using geometric combination rules in "
4305 "LJ-PME, but your non-bonded C6 parameters do "
4306 "not follow these rules. "
4307 "This will introduce very small errors in the forces and energies in "
4308 "your simulations. If your system is homogeneous, consider using "
4309 "dispersion correction "
4310 "for the total energy and pressure.");
4316 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4318 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4319 GMX_ASSERT(gmx::checkMtsRequirements(*ir).empty(), "All MTS requirements should be met here");
4321 char err_buf[STRLEN];
4324 gmx_mtop_atomloop_block_t aloopb;
4326 char warn_buf[STRLEN];
4328 set_warning_line(wi, mdparin, -1);
4330 if (absolute_reference(ir, sys, false, AbsRef))
4333 "Removing center of mass motion in the presence of position restraints might "
4334 "cause artifacts. When you are using position restraints to equilibrate a "
4335 "macro-molecule, the artifacts are usually negligible.");
4338 if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 && ir->nstlist > 1
4339 && ((EI_MD(ir->eI) || EI_SD(ir->eI))
4340 && (ir->etc == TemperatureCoupling::VRescale || ir->etc == TemperatureCoupling::Berendsen)))
4342 /* Check if a too small Verlet buffer might potentially
4343 * cause more drift than the thermostat can couple off.
4345 /* Temperature error fraction for warning and suggestion */
4346 const real T_error_warn = 0.002;
4347 const real T_error_suggest = 0.001;
4348 /* For safety: 2 DOF per atom (typical with constraints) */
4349 const real nrdf_at = 2;
4350 real T, tau, max_T_error;
4355 for (i = 0; i < ir->opts.ngtc; i++)
4357 T = std::max(T, ir->opts.ref_t[i]);
4358 tau = std::max(tau, ir->opts.tau_t[i]);
4362 /* This is a worst case estimate of the temperature error,
4363 * assuming perfect buffer estimation and no cancelation
4364 * of errors. The factor 0.5 is because energy distributes
4365 * equally over Ekin and Epot.
4367 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * BOLTZ * T);
4368 if (max_T_error > T_error_warn)
4371 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4372 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4373 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4374 "%.0e or decrease tau_t.",
4379 100 * T_error_suggest,
4380 ir->verletbuf_tol * T_error_suggest / max_T_error);
4381 warning(wi, warn_buf);
4386 if (ETC_ANDERSEN(ir->etc))
4390 for (i = 0; i < ir->opts.ngtc; i++)
4393 "all tau_t must currently be equal using Andersen temperature control, "
4394 "violated for group %d",
4396 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4398 "all tau_t must be positive using Andersen temperature control, "
4402 CHECK(ir->opts.tau_t[i] < 0);
4405 if (ir->etc == TemperatureCoupling::AndersenMassive && ir->comm_mode != ecmNO)
4407 for (i = 0; i < ir->opts.ngtc; i++)
4409 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4411 "tau_t/delta_t for group %d for temperature control method %s must be a "
4412 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4413 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4416 enumValueToString(ir->etc),
4420 CHECK(nsteps % ir->nstcomm != 0);
4425 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD && ir->comm_mode == ecmNO
4426 && !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) && !ETC_ANDERSEN(ir->etc))
4429 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4430 "rounding errors can lead to build up of kinetic energy of the center of mass");
4433 if (ir->epc == PressureCoupling::ParrinelloRahman && ir->etc == TemperatureCoupling::NoseHoover)
4436 for (int g = 0; g < ir->opts.ngtc; g++)
4438 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4440 if (ir->tau_p < 1.9 * tau_t_max)
4442 std::string message = gmx::formatString(
4443 "With %s T-coupling and %s p-coupling, "
4444 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4445 enumValueToString(ir->etc),
4446 enumValueToString(ir->epc),
4451 warning(wi, message.c_str());
4455 /* Check for pressure coupling with absolute position restraints */
4456 if (ir->epc != PressureCoupling::No && ir->refcoord_scaling == erscNO)
4458 absolute_reference(ir, sys, TRUE, AbsRef);
4460 for (m = 0; m < DIM; m++)
4462 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4465 "You are using pressure coupling with absolute position restraints, "
4466 "this will give artifacts. Use the refcoord_scaling option.");
4474 aloopb = gmx_mtop_atomloop_block_init(sys);
4476 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4478 if (atom->q != 0 || atom->qB != 0)
4486 if (EEL_FULL(ir->coulombtype))
4489 "You are using full electrostatics treatment %s for a system without charges.\n"
4490 "This costs a lot of performance for just processing zeros, consider using %s "
4492 EELTYPE(ir->coulombtype),
4494 warning(wi, err_buf);
4499 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4502 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4503 "You might want to consider using %s electrostatics.\n",
4505 warning_note(wi, err_buf);
4509 /* Check if combination rules used in LJ-PME are the same as in the force field */
4510 if (EVDW_PME(ir->vdwtype))
4512 check_combination_rules(ir, sys, wi);
4515 /* Generalized reaction field */
4516 if (ir->coulombtype == eelGRF_NOTUSED)
4519 "Generalized reaction-field electrostatics is no longer supported. "
4520 "You can use normal reaction-field instead and compute the reaction-field "
4521 "constant by hand.");
4524 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0
4525 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4527 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4535 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4537 if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
4539 absolute_reference(ir, sys, FALSE, AbsRef);
4540 for (m = 0; m < DIM; m++)
4542 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4545 "You are using an absolute reference for pulling, but the rest of "
4546 "the system does not have an absolute reference. This will lead to "
4555 for (i = 0; i < 3; i++)
4557 for (m = 0; m <= i; m++)
4559 if ((ir->epc != PressureCoupling::No && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4561 for (c = 0; c < ir->pull->ncoord; c++)
4563 if (ir->pull->coord[c].eGeom == epullgDIRPBC && ir->pull->coord[c].vec[m] != 0)
4566 "Can not have dynamic box while using pull geometry '%s' "
4568 EPULLGEOM(ir->pull->coord[c].eGeom),
4580 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4582 char warn_buf[STRLEN];
4585 ptr = check_box(ir->pbcType, box);
4588 warning_error(wi, ptr);
4591 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4593 if (ir->shake_tol <= 0.0)
4595 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4596 warning_error(wi, warn_buf);
4600 if ((ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4602 /* If we have Lincs constraints: */
4603 if (ir->eI == eiMD && ir->etc == TemperatureCoupling::No && ir->eConstrAlg == econtLINCS
4604 && ir->nLincsIter == 1)
4607 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4608 warning_note(wi, warn_buf);
4611 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4614 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4616 warning_note(wi, warn_buf);
4618 if (ir->epc == PressureCoupling::Mttk)
4620 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4624 if (bHasAnyConstraints && ir->epc == PressureCoupling::Mttk)
4626 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4629 if (ir->LincsWarnAngle > 90.0)
4631 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4632 warning(wi, warn_buf);
4633 ir->LincsWarnAngle = 90.0;
4636 if (ir->pbcType != PbcType::No)
4638 if (ir->nstlist == 0)
4641 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4642 "atoms might cause the simulation to crash.");
4644 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4647 "ERROR: The cut-off length is longer than half the shortest box vector or "
4648 "longer than the smallest box diagonal element. Increase the box size or "
4649 "decrease rlist.\n");
4650 warning_error(wi, warn_buf);