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.
51 #include "gromacs/applied_forces/awh/read_params.h"
52 #include "gromacs/fileio/readinp.h"
53 #include "gromacs/fileio/warninp.h"
54 #include "gromacs/gmxlib/network.h"
55 #include "gromacs/gmxpreprocess/toputil.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/calc_verletbuf.h"
60 #include "gromacs/mdrun/mdmodules.h"
61 #include "gromacs/mdtypes/inputrec.h"
62 #include "gromacs/mdtypes/md_enums.h"
63 #include "gromacs/mdtypes/multipletimestepping.h"
64 #include "gromacs/mdtypes/pull_params.h"
65 #include "gromacs/options/options.h"
66 #include "gromacs/options/treesupport.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/selection/indexutil.h"
69 #include "gromacs/topology/block.h"
70 #include "gromacs/topology/ifunc.h"
71 #include "gromacs/topology/index.h"
72 #include "gromacs/topology/mtop_util.h"
73 #include "gromacs/topology/symtab.h"
74 #include "gromacs/topology/topology.h"
75 #include "gromacs/utility/arrayref.h"
76 #include "gromacs/utility/cstringutil.h"
77 #include "gromacs/utility/exceptions.h"
78 #include "gromacs/utility/fatalerror.h"
79 #include "gromacs/utility/filestream.h"
80 #include "gromacs/utility/gmxassert.h"
81 #include "gromacs/utility/ikeyvaluetreeerror.h"
82 #include "gromacs/utility/keyvaluetree.h"
83 #include "gromacs/utility/keyvaluetreebuilder.h"
84 #include "gromacs/utility/keyvaluetreemdpwriter.h"
85 #include "gromacs/utility/keyvaluetreetransform.h"
86 #include "gromacs/utility/mdmodulenotification.h"
87 #include "gromacs/utility/smalloc.h"
88 #include "gromacs/utility/strconvert.h"
89 #include "gromacs/utility/stringcompare.h"
90 #include "gromacs/utility/stringutil.h"
91 #include "gromacs/utility/textwriter.h"
96 using gmx::BasicVector;
98 /* Resource parameters
99 * Do not change any of these until you read the instruction
100 * in readinp.h. Some cpp's do not take spaces after the backslash
101 * (like the c-shell), which will give you a very weird compiler
105 struct gmx_inputrec_strings
107 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN], acc[STRLEN], accgrps[STRLEN], freeze[STRLEN],
108 frdim[STRLEN], energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN],
109 x_compressed_groups[STRLEN], couple_moltype[STRLEN], orirefitgrp[STRLEN],
110 egptable[STRLEN], egpexcl[STRLEN], wall_atomtype[STRLEN], wall_density[STRLEN],
111 deform[STRLEN], QMMM[STRLEN], imd_grp[STRLEN];
112 char fep_lambda[efptNR][STRLEN];
113 char lambda_weights[STRLEN];
114 std::vector<std::string> pullGroupNames;
115 std::vector<std::string> rotateGroupNames;
116 char anneal[STRLEN], anneal_npoints[STRLEN], anneal_time[STRLEN], anneal_temp[STRLEN];
119 static gmx_inputrec_strings* inputrecStrings = nullptr;
121 void init_inputrec_strings()
126 "Attempted to call init_inputrec_strings before calling done_inputrec_strings. "
127 "Only one inputrec (i.e. .mdp file) can be parsed at a time.");
129 inputrecStrings = new gmx_inputrec_strings();
132 void done_inputrec_strings()
134 delete inputrecStrings;
135 inputrecStrings = nullptr;
141 egrptpALL, /* All particles have to be a member of a group. */
142 egrptpALL_GENREST, /* A rest group with name is generated for particles *
143 * that are not part of any group. */
144 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
145 * for the rest group. */
146 egrptpONE /* Merge all selected groups into one group, *
147 * make a rest group for the remaining particles. */
150 static const char* constraints[eshNR + 1] = { "none", "h-bonds", "all-bonds",
151 "h-angles", "all-angles", nullptr };
153 static const char* couple_lam[ecouplamNR + 1] = { "vdw-q", "vdw", "q", "none", nullptr };
155 static void GetSimTemps(int ntemps, t_simtemp* simtemp, double* temperature_lambdas)
160 for (i = 0; i < ntemps; i++)
162 /* simple linear scaling -- allows more control */
163 if (simtemp->eSimTempScale == esimtempLINEAR)
165 simtemp->temperatures[i] =
167 + (simtemp->simtemp_high - simtemp->simtemp_low) * temperature_lambdas[i];
169 else if (simtemp->eSimTempScale
170 == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
172 simtemp->temperatures[i] = simtemp->simtemp_low
173 * std::pow(simtemp->simtemp_high / simtemp->simtemp_low,
174 static_cast<real>((1.0 * i) / (ntemps - 1)));
176 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
178 simtemp->temperatures[i] = simtemp->simtemp_low
179 + (simtemp->simtemp_high - simtemp->simtemp_low)
180 * (std::expm1(temperature_lambdas[i]) / std::expm1(1.0));
185 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
186 gmx_fatal(FARGS, "%s", errorstr);
192 static void _low_check(bool b, const char* s, warninp_t wi)
196 warning_error(wi, s);
200 static void check_nst(const char* desc_nst, int nst, const char* desc_p, int* p, warninp_t wi)
204 if (*p > 0 && *p % nst != 0)
206 /* Round up to the next multiple of nst */
207 *p = ((*p) / nst + 1) * nst;
208 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n", desc_p, desc_nst, desc_p, *p);
213 static int lcd(int n1, int n2)
218 for (i = 2; (i <= n1 && i <= n2); i++)
220 if (n1 % i == 0 && n2 % i == 0)
229 //! Convert legacy mdp entries to modern ones.
230 static void process_interaction_modifier(int* eintmod)
232 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
234 *eintmod = eintmodPOTSHIFT;
238 static void checkMtsRequirement(const t_inputrec& ir, const char* param, const int nstValue, warninp_t wi)
240 GMX_RELEASE_ASSERT(ir.mtsLevels.size() >= 2, "Need at least two levels for MTS");
241 const int mtsFactor = ir.mtsLevels.back().stepFactor;
242 if (nstValue % mtsFactor != 0)
244 auto message = gmx::formatString(
245 "With MTS, %s = %d should be a multiple of mts-factor = %d", param, nstValue, mtsFactor);
246 warning_error(wi, message.c_str());
250 static void setupMtsLevels(gmx::ArrayRef<gmx::MtsLevel> mtsLevels,
251 const t_inputrec& ir,
252 const t_gromppopts& opts,
255 /* MD-VV has no MTS support yet.
256 * SD1 needs different scaling coefficients for the different MTS forces
257 * and the different forces are currently not available in ForceBuffers.
261 auto message = gmx::formatString(
262 "Multiple time stepping is only supported with integrator %s", ei_names[eiMD]);
263 warning_error(wi, message.c_str());
265 if (opts.numMtsLevels != 2)
267 warning_error(wi, "Only mts-levels = 2 is supported");
271 const std::vector<std::string> inputForceGroups = gmx::splitString(opts.mtsLevel2Forces);
272 auto& forceGroups = mtsLevels[1].forceGroups;
273 for (const auto& inputForceGroup : inputForceGroups)
277 for (const auto& forceGroupName : gmx::mtsForceGroupNames)
279 if (gmx::equalCaseInsensitive(inputForceGroup, forceGroupName))
281 forceGroups.set(nameIndex);
289 gmx::formatString("Unknown MTS force group '%s'", inputForceGroup.c_str());
290 warning_error(wi, message.c_str());
294 if (mtsLevels[1].stepFactor <= 1)
296 gmx_fatal(FARGS, "mts-factor should be larger than 1");
299 // Make the level 0 use the complement of the force groups of group 1
300 mtsLevels[0].forceGroups = ~mtsLevels[1].forceGroups;
301 mtsLevels[0].stepFactor = 1;
303 if ((EEL_FULL(ir.coulombtype) || EVDW_PME(ir.vdwtype))
304 && !mtsLevels[1].forceGroups[static_cast<int>(gmx::MtsForceGroups::LongrangeNonbonded)])
307 "With long-range electrostatics and/or LJ treatment, the long-range part "
308 "has to be part of the mts-level2-forces");
311 if (ir.nstcalcenergy > 0)
313 checkMtsRequirement(ir, "nstcalcenergy", ir.nstcalcenergy, wi);
315 checkMtsRequirement(ir, "nstenergy", ir.nstenergy, wi);
316 checkMtsRequirement(ir, "nstlog", ir.nstlog, wi);
317 if (ir.efep != efepNO)
319 checkMtsRequirement(ir, "nstdhdl", ir.fepvals->nstdhdl, wi);
324 const int pullMtsLevel = gmx::forceGroupMtsLevel(ir.mtsLevels, gmx::MtsForceGroups::Pull);
325 if (ir.pull->nstxout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
327 warning_error(wi, "pull-nstxout should be a multiple of mts-factor");
329 if (ir.pull->nstfout % ir.mtsLevels[pullMtsLevel].stepFactor != 0)
331 warning_error(wi, "pull-nstfout should be a multiple of mts-factor");
337 void check_ir(const char* mdparin,
338 const gmx::MdModulesNotifier& mdModulesNotifier,
342 /* Check internal consistency.
343 * NOTE: index groups are not set here yet, don't check things
344 * like temperature coupling group options here, but in triple_check
347 /* Strange macro: first one fills the err_buf, and then one can check
348 * the condition, which will print the message and increase the error
351 #define CHECK(b) _low_check(b, err_buf, wi)
352 char err_buf[256], warn_buf[STRLEN];
355 t_lambda* fep = ir->fepvals;
356 t_expanded* expand = ir->expandedvals;
358 set_warning_line(wi, mdparin, -1);
360 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
362 sprintf(warn_buf, "%s electrostatics is no longer supported", eel_names[eelRF_NEC_UNSUPPORTED]);
363 warning_error(wi, warn_buf);
366 /* BASIC CUT-OFF STUFF */
367 if (ir->rcoulomb < 0)
369 warning_error(wi, "rcoulomb should be >= 0");
373 warning_error(wi, "rvdw should be >= 0");
375 if (ir->rlist < 0 && !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
377 warning_error(wi, "rlist should be >= 0");
380 "nstlist can not be smaller than 0. (If you were trying to use the heuristic "
381 "neighbour-list update scheme for efficient buffering for improved energy "
382 "conservation, please use the Verlet cut-off scheme instead.)");
383 CHECK(ir->nstlist < 0);
385 process_interaction_modifier(&ir->coulomb_modifier);
386 process_interaction_modifier(&ir->vdw_modifier);
388 if (ir->cutoff_scheme == ecutsGROUP)
391 "The group cutoff scheme has been removed since GROMACS 2020. "
392 "Please use the Verlet cutoff scheme.");
394 if (ir->cutoff_scheme == ecutsVERLET)
398 /* Normal Verlet type neighbor-list, currently only limited feature support */
399 if (inputrec2nboundeddim(ir) < 3)
401 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
404 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
405 if (ir->rcoulomb != ir->rvdw)
407 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
408 // for PME load balancing, we can support this exception.
409 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) && ir->vdwtype == evdwCUT
410 && ir->rcoulomb > ir->rvdw);
411 if (!bUsesPmeTwinRangeKernel)
414 "With Verlet lists rcoulomb!=rvdw is not supported (except for "
415 "rcoulomb>rvdw with PME electrostatics)");
419 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
421 if (ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT)
423 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
426 "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and "
428 evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
429 warning_note(wi, warn_buf);
431 ir->vdwtype = evdwCUT;
435 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s",
436 evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
437 warning_error(wi, warn_buf);
441 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
444 "With Verlet lists only cut-off and PME LJ interactions are supported");
446 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) || EEL_PME(ir->coulombtype)
447 || ir->coulombtype == eelEWALD))
450 "With Verlet lists only cut-off, reaction-field, PME and Ewald "
451 "electrostatics are supported");
453 if (!(ir->coulomb_modifier == eintmodNONE || ir->coulomb_modifier == eintmodPOTSHIFT))
455 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
456 warning_error(wi, warn_buf);
459 if (EEL_USER(ir->coulombtype))
461 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme",
462 eel_names[ir->coulombtype]);
463 warning_error(wi, warn_buf);
466 if (ir->nstlist <= 0)
468 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
471 if (ir->nstlist < 10)
474 "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note "
475 "that with the Verlet scheme, nstlist has no effect on the accuracy of "
479 rc_max = std::max(ir->rvdw, ir->rcoulomb);
483 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
484 ir->verletbuf_tol = 0;
487 else if (ir->verletbuf_tol <= 0)
489 if (ir->verletbuf_tol == 0)
491 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
494 if (ir->rlist < rc_max)
497 "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
500 if (ir->rlist == rc_max && ir->nstlist > 1)
504 "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet "
505 "buffer. The cluster pair list does have a buffering effect, but choosing "
506 "a larger rlist might be necessary for good energy conservation.");
511 if (ir->rlist > rc_max)
514 "You have set rlist larger than the interaction cut-off, but you also "
515 "have verlet-buffer-tolerance > 0. Will set rlist using "
516 "verlet-buffer-tolerance.");
519 if (ir->nstlist == 1)
521 /* No buffer required */
526 if (EI_DYNAMICS(ir->eI))
528 if (inputrec2nboundeddim(ir) < 3)
531 "The box volume is required for calculating rlist from the "
532 "energy drift with verlet-buffer-tolerance > 0. You are "
533 "using at least one unbounded dimension, so no volume can be "
534 "computed. Either use a finite box, or set rlist yourself "
535 "together with verlet-buffer-tolerance = -1.");
537 /* Set rlist temporarily so we can continue processing */
542 /* Set the buffer to 5% of the cut-off */
543 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics) * rc_max;
549 /* GENERAL INTEGRATOR STUFF */
552 if (ir->etc != etcNO)
554 if (EI_RANDOM(ir->eI))
557 "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling "
558 "implicitly. See the documentation for more information on which "
559 "parameters affect temperature for %s.",
560 etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
565 "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to "
567 etcoupl_names[ir->etc], ei_names[ir->eI]);
569 warning_note(wi, warn_buf);
573 if (ir->eI == eiVVAK)
576 "Integrator method %s is implemented primarily for validation purposes; for "
577 "molecular dynamics, you should probably be using %s or %s",
578 ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
579 warning_note(wi, warn_buf);
581 if (!EI_DYNAMICS(ir->eI))
583 if (ir->epc != epcNO)
586 "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.",
587 epcoupl_names[ir->epc], ei_names[ir->eI]);
588 warning_note(wi, warn_buf);
592 if (EI_DYNAMICS(ir->eI))
594 if (ir->nstcalcenergy < 0)
596 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
597 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
599 /* nstcalcenergy larger than nstener does not make sense.
600 * We ideally want nstcalcenergy=nstener.
604 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
608 ir->nstcalcenergy = ir->nstenergy;
612 else if ((ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy)
613 || (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
614 && (ir->nstcalcenergy > ir->fepvals->nstdhdl)))
617 const char* nsten = "nstenergy";
618 const char* nstdh = "nstdhdl";
619 const char* min_name = nsten;
620 int min_nst = ir->nstenergy;
622 /* find the smallest of ( nstenergy, nstdhdl ) */
623 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0
624 && (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
626 min_nst = ir->fepvals->nstdhdl;
629 /* If the user sets nstenergy small, we should respect that */
630 sprintf(warn_buf, "Setting nstcalcenergy (%d) equal to %s (%d)", ir->nstcalcenergy,
632 warning_note(wi, warn_buf);
633 ir->nstcalcenergy = min_nst;
636 if (ir->epc != epcNO)
638 if (ir->nstpcouple < 0)
640 ir->nstpcouple = ir_optimal_nstpcouple(ir);
642 if (ir->useMts && ir->nstpcouple % ir->mtsLevels.back().stepFactor != 0)
645 "With multiple time stepping, nstpcouple should be a mutiple of "
650 if (ir->nstcalcenergy > 0)
652 if (ir->efep != efepNO)
654 /* nstdhdl should be a multiple of nstcalcenergy */
655 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstdhdl", &ir->fepvals->nstdhdl, wi);
659 /* nstexpanded should be a multiple of nstcalcenergy */
660 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstexpanded",
661 &ir->expandedvals->nstexpanded, wi);
663 /* for storing exact averages nstenergy should be
664 * a multiple of nstcalcenergy
666 check_nst("nstcalcenergy", ir->nstcalcenergy, "nstenergy", &ir->nstenergy, wi);
669 // Inquire all MdModules, if their parameters match with the energy
670 // calculation frequency
671 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
672 mdModulesNotifier.preProcessingNotifications_.notify(&energyCalculationFrequencyErrors);
674 // Emit all errors from the energy calculation frequency checks
675 for (const std::string& energyFrequencyErrorMessage :
676 energyCalculationFrequencyErrors.errorMessages())
678 warning_error(wi, energyFrequencyErrorMessage);
682 if (ir->nsteps == 0 && !ir->bContinuation)
685 "For a correct single-point energy evaluation with nsteps = 0, use "
686 "continuation = yes to avoid constraining the input coordinates.");
690 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
693 "You are doing a continuation with SD or BD, make sure that ld_seed is "
694 "different from the previous run (using ld_seed=-1 will ensure this)");
700 sprintf(err_buf, "TPI only works with pbc = %s", c_pbcTypeNames[PbcType::Xyz].c_str());
701 CHECK(ir->pbcType != PbcType::Xyz);
702 sprintf(err_buf, "with TPI nstlist should be larger than zero");
703 CHECK(ir->nstlist <= 0);
704 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
705 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
709 if ((opts->nshake > 0) && (opts->bMorse))
711 sprintf(warn_buf, "Using morse bond-potentials while constraining bonds is useless");
712 warning(wi, warn_buf);
715 if ((EI_SD(ir->eI) || ir->eI == eiBD) && ir->bContinuation && ir->ld_seed != -1)
718 "You are doing a continuation with SD or BD, make sure that ld_seed is "
719 "different from the previous run (using ld_seed=-1 will ensure this)");
721 /* verify simulated tempering options */
725 bool bAllTempZero = TRUE;
726 for (i = 0; i < fep->n_lambda; i++)
728 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
729 efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
730 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
731 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
733 bAllTempZero = FALSE;
736 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
737 CHECK(bAllTempZero == TRUE);
739 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
740 CHECK(ir->eI != eiVV);
742 /* check compatability of the temperature coupling with simulated tempering */
744 if (ir->etc == etcNOSEHOOVER)
747 "Nose-Hoover based temperature control such as [%s] my not be "
748 "entirelyconsistent with simulated tempering",
749 etcoupl_names[ir->etc]);
750 warning_note(wi, warn_buf);
753 /* check that the temperatures make sense */
756 "Higher simulated tempering temperature (%g) must be >= than the simulated "
757 "tempering lower temperature (%g)",
758 ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
759 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
761 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero",
762 ir->simtempvals->simtemp_high);
763 CHECK(ir->simtempvals->simtemp_high <= 0);
765 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero",
766 ir->simtempvals->simtemp_low);
767 CHECK(ir->simtempvals->simtemp_low <= 0);
770 /* verify free energy options */
772 if (ir->efep != efepNO)
775 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2", fep->sc_power);
776 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
779 "The soft-core sc-r-power is %d and can only be 6. (sc-r-power 48 is no longer "
781 static_cast<int>(fep->sc_r_power));
782 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0);
785 "Can't use positive delta-lambda (%g) if initial state/lambda does not start at "
788 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
790 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations",
792 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
794 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
795 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
797 sprintf(err_buf, "Free-energy not implemented for Ewald");
798 CHECK(ir->coulombtype == eelEWALD);
800 /* check validty of lambda inputs */
801 if (fep->n_lambda == 0)
803 /* Clear output in case of no states:*/
804 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.",
805 fep->init_fep_state);
806 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
810 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d",
811 fep->init_fep_state, fep->n_lambda - 1);
812 CHECK((fep->init_fep_state >= fep->n_lambda));
816 "Lambda state must be set, either with init-lambda-state or with init-lambda");
817 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
820 "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with "
821 "init-lambda-state or with init-lambda, but not both",
822 fep->init_lambda, fep->init_fep_state);
823 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
826 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
830 for (i = 0; i < efptNR; i++)
832 if (fep->separate_dvdl[i])
837 if (n_lambda_terms > 1)
840 "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't "
841 "use init-lambda to set lambda state (except for slow growth). Use "
842 "init-lambda-state instead.");
843 warning(wi, warn_buf);
846 if (n_lambda_terms < 2 && fep->n_lambda > 0)
849 "init-lambda is deprecated for setting lambda state (except for slow "
850 "growth). Use init-lambda-state instead.");
854 for (j = 0; j < efptNR; j++)
856 for (i = 0; i < fep->n_lambda; i++)
858 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i,
859 efpt_names[j], fep->all_lambda[j][i]);
860 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
864 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
866 for (i = 0; i < fep->n_lambda; i++)
869 "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while "
870 "coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to "
871 "crashes, and is not supported.",
872 i, fep->all_lambda[efptVDW][i], fep->all_lambda[efptCOUL][i]);
873 CHECK((fep->sc_alpha > 0)
874 && (((fep->all_lambda[efptCOUL][i] > 0.0) && (fep->all_lambda[efptCOUL][i] < 1.0))
875 && ((fep->all_lambda[efptVDW][i] > 0.0) && (fep->all_lambda[efptVDW][i] < 1.0))));
879 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
881 real sigma, lambda, r_sc;
884 /* Maximum estimate for A and B charges equal with lambda power 1 */
886 r_sc = std::pow(lambda * fep->sc_alpha * std::pow(sigma / ir->rcoulomb, fep->sc_r_power) + 1.0,
887 1.0 / fep->sc_r_power);
889 "With PME there is a minor soft core effect present at the cut-off, "
890 "proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on "
891 "energy conservation, but usually other effects dominate. With a common sigma "
892 "value of %g nm the fraction of the particle-particle potential at the cut-off "
893 "at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
894 fep->sc_r_power, sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
895 warning_note(wi, warn_buf);
898 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
899 be treated differently, but that's the next step */
901 for (i = 0; i < efptNR; i++)
903 for (j = 0; j < fep->n_lambda; j++)
905 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
906 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
911 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
915 /* checking equilibration of weights inputs for validity */
918 "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal "
920 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
921 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
924 "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to "
926 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
927 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
930 "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
931 expand->equil_steps, elmceq_names[elmceqSTEPS]);
932 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
935 "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
936 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
937 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
940 "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
941 expand->equil_ratio, elmceq_names[elmceqRATIO]);
942 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
945 "weight-equil-number-all-lambda (%d) must be a positive integer if "
946 "lmc-weights-equil=%s",
947 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
948 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
951 "weight-equil-number-samples (%d) must be a positive integer if "
952 "lmc-weights-equil=%s",
953 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
954 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
957 "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
958 expand->equil_steps, elmceq_names[elmceqSTEPS]);
959 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
961 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
962 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
963 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
965 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
966 expand->equil_ratio, elmceq_names[elmceqRATIO]);
967 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
969 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
970 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
971 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
973 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
974 CHECK((expand->lmc_repeats <= 0));
975 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
976 CHECK((expand->minvarmin <= 0));
977 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
978 CHECK((expand->c_range < 0));
980 "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != "
982 fep->init_fep_state, expand->lmc_forced_nstart);
983 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0)
984 && (expand->elmcmove != elmcmoveNO));
985 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
986 CHECK((expand->lmc_forced_nstart < 0));
987 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)",
988 fep->init_fep_state);
989 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
991 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
992 CHECK((expand->init_wl_delta < 0));
993 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
994 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
995 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
996 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
998 /* if there is no temperature control, we need to specify an MC temperature */
999 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO)
1000 && (expand->mc_temp <= 0.0))
1003 "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set "
1004 "to a positive number");
1005 warning_error(wi, err_buf);
1007 if (expand->nstTij > 0)
1009 sprintf(err_buf, "nstlog must be non-zero");
1010 CHECK(ir->nstlog == 0);
1011 // Avoid modulus by zero in the case that already triggered an error exit.
1012 if (ir->nstlog != 0)
1015 "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
1016 expand->nstTij, ir->nstlog);
1017 CHECK((expand->nstTij % ir->nstlog) != 0);
1023 sprintf(err_buf, "walls only work with pbc=%s", c_pbcTypeNames[PbcType::XY].c_str());
1024 CHECK(ir->nwall && ir->pbcType != PbcType::XY);
1027 if (ir->pbcType != PbcType::Xyz && ir->nwall != 2)
1029 if (ir->pbcType == PbcType::No)
1031 if (ir->epc != epcNO)
1033 warning(wi, "Turning off pressure coupling for vacuum system");
1039 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
1040 c_pbcTypeNames[ir->pbcType].c_str());
1041 CHECK(ir->epc != epcNO);
1043 sprintf(err_buf, "Can not have Ewald with pbc=%s", c_pbcTypeNames[ir->pbcType].c_str());
1044 CHECK(EEL_FULL(ir->coulombtype));
1046 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
1047 c_pbcTypeNames[ir->pbcType].c_str());
1048 CHECK(ir->eDispCorr != edispcNO);
1051 if (ir->rlist == 0.0)
1054 "can only have neighborlist cut-off zero (=infinite)\n"
1055 "with coulombtype = %s or coulombtype = %s\n"
1056 "without periodic boundary conditions (pbc = %s) and\n"
1057 "rcoulomb and rvdw set to zero",
1058 eel_names[eelCUT], eel_names[eelUSER], c_pbcTypeNames[PbcType::No].c_str());
1059 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER))
1060 || (ir->pbcType != PbcType::No) || (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
1062 if (ir->nstlist > 0)
1065 "Simulating without cut-offs can be (slightly) faster with nstlist=0, "
1066 "nstype=simple and only one MPI rank");
1071 if (ir->nstcomm == 0)
1073 // TODO Change this behaviour. There should be exactly one way
1074 // to turn off an algorithm.
1075 ir->comm_mode = ecmNO;
1077 if (ir->comm_mode != ecmNO)
1079 if (ir->nstcomm < 0)
1081 // TODO Such input was once valid. Now that we've been
1082 // helpful for a few years, we should reject such input,
1083 // lest we have to support every historical decision
1086 "If you want to remove the rotation around the center of mass, you should set "
1087 "comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to "
1088 "its absolute value");
1089 ir->nstcomm = abs(ir->nstcomm);
1092 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
1095 "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting "
1096 "nstcomm to nstcalcenergy");
1097 ir->nstcomm = ir->nstcalcenergy;
1100 if (ir->comm_mode == ecmANGULAR)
1103 "Can not remove the rotation around the center of mass with periodic "
1105 CHECK(ir->bPeriodicMols);
1106 if (ir->pbcType != PbcType::No)
1109 "Removing the rotation around the center of mass in a periodic system, "
1110 "this can lead to artifacts. Only use this on a single (cluster of) "
1111 "molecules. This cluster should not cross periodic boundaries.");
1116 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->pbcType == PbcType::No && ir->comm_mode != ecmANGULAR)
1119 "Tumbling and flying ice-cubes: We are not removing rotation around center of mass "
1120 "in a non-periodic system. You should probably set comm_mode = ANGULAR or use "
1123 warning_note(wi, warn_buf);
1126 /* TEMPERATURE COUPLING */
1127 if (ir->etc == etcYES)
1129 ir->etc = etcBERENDSEN;
1131 "Old option for temperature coupling given: "
1132 "changing \"yes\" to \"Berendsen\"\n");
1135 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
1137 if (ir->opts.nhchainlength < 1)
1140 "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to "
1142 ir->opts.nhchainlength);
1143 ir->opts.nhchainlength = 1;
1144 warning(wi, warn_buf);
1147 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
1151 "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
1152 ir->opts.nhchainlength = 1;
1157 ir->opts.nhchainlength = 0;
1160 if (ir->eI == eiVVAK)
1163 "%s implemented primarily for validation, and requires nsttcouple = 1 and "
1166 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
1169 if (ETC_ANDERSEN(ir->etc))
1171 sprintf(err_buf, "%s temperature control not supported for integrator %s.",
1172 etcoupl_names[ir->etc], ei_names[ir->eI]);
1173 CHECK(!(EI_VV(ir->eI)));
1175 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
1178 "Center of mass removal not necessary for %s. All velocities of coupled "
1179 "groups are rerandomized periodically, so flying ice cube errors will not "
1181 etcoupl_names[ir->etc]);
1182 warning_note(wi, warn_buf);
1186 "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are "
1187 "randomized every time step",
1188 ir->nstcomm, etcoupl_names[ir->etc]);
1189 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
1192 if (ir->etc == etcBERENDSEN)
1195 "The %s thermostat does not generate the correct kinetic energy distribution. You "
1196 "might want to consider using the %s thermostat.",
1197 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
1198 warning_note(wi, warn_buf);
1201 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc)) && ir->epc == epcBERENDSEN)
1204 "Using Berendsen pressure coupling invalidates the "
1205 "true ensemble for the thermostat");
1206 warning(wi, warn_buf);
1209 /* PRESSURE COUPLING */
1210 if (ir->epc == epcISOTROPIC)
1212 ir->epc = epcBERENDSEN;
1214 "Old option for pressure coupling given: "
1215 "changing \"Isotropic\" to \"Berendsen\"\n");
1218 if (ir->epc != epcNO)
1220 dt_pcoupl = ir->nstpcouple * ir->delta_t;
1222 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1223 CHECK(ir->tau_p <= 0);
1225 if (ir->tau_p / dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10 * GMX_REAL_EPS)
1228 "For proper integration of the %s barostat, tau-p (%g) should be at least %d "
1229 "times larger than nstpcouple*dt (%g)",
1230 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1231 warning(wi, warn_buf);
1235 "compressibility must be > 0 when using pressure"
1237 EPCOUPLTYPE(ir->epc));
1238 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 || ir->compress[ZZ][ZZ] < 0
1239 || (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 && ir->compress[ZZ][XX] <= 0
1240 && ir->compress[ZZ][YY] <= 0));
1242 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1245 "You are generating velocities so I am assuming you "
1246 "are equilibrating a system. You are using "
1247 "%s pressure coupling, but this can be "
1248 "unstable for equilibration. If your system crashes, try "
1249 "equilibrating first with Berendsen pressure coupling. If "
1250 "you are not equilibrating the system, you can probably "
1251 "ignore this warning.",
1252 epcoupl_names[ir->epc]);
1253 warning(wi, warn_buf);
1259 if (ir->epc == epcMTTK)
1261 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1265 /* ELECTROSTATICS */
1266 /* More checks are in triple check (grompp.c) */
1268 if (ir->coulombtype == eelSWITCH)
1271 "coulombtype = %s is only for testing purposes and can lead to serious "
1272 "artifacts, advice: use coulombtype = %s",
1273 eel_names[ir->coulombtype], eel_names[eelRF_ZERO]);
1274 warning(wi, warn_buf);
1277 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1280 "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old "
1281 "format and exchanging epsilon-r and epsilon-rf",
1283 warning(wi, warn_buf);
1284 ir->epsilon_rf = ir->epsilon_r;
1285 ir->epsilon_r = 1.0;
1288 if (ir->epsilon_r == 0)
1291 "It is pointless to use long-range electrostatics with infinite relative "
1293 "Since you are effectively turning of electrostatics, a plain cutoff will be much "
1295 CHECK(EEL_FULL(ir->coulombtype));
1298 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1300 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1301 CHECK(ir->epsilon_r < 0);
1304 if (EEL_RF(ir->coulombtype))
1306 /* reaction field (at the cut-off) */
1308 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1311 "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1312 eel_names[ir->coulombtype]);
1313 warning(wi, warn_buf);
1314 ir->epsilon_rf = 0.0;
1317 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1318 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) || (ir->epsilon_r == 0));
1319 if (ir->epsilon_rf == ir->epsilon_r)
1321 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1322 eel_names[ir->coulombtype]);
1323 warning(wi, warn_buf);
1326 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1327 * means the interaction is zero outside rcoulomb, but it helps to
1328 * provide accurate energy conservation.
1330 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1332 if (ir_coulomb_switched(ir))
1335 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the "
1336 "potential modifier options!",
1337 eel_names[ir->coulombtype]);
1338 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1342 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1345 "Explicit switch/shift coulomb interactions cannot be used in combination with a "
1346 "secondary coulomb-modifier.");
1347 CHECK(ir->coulomb_modifier != eintmodNONE);
1349 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1352 "Explicit switch/shift vdw interactions cannot be used in combination with a "
1353 "secondary vdw-modifier.");
1354 CHECK(ir->vdw_modifier != eintmodNONE);
1357 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT || ir->vdwtype == evdwSWITCH
1358 || ir->vdwtype == evdwSHIFT)
1361 "The switch/shift interaction settings are just for compatibility; you will get "
1363 "performance from applying potential modifiers to your interactions!\n");
1364 warning_note(wi, warn_buf);
1367 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1369 if (ir->rcoulomb_switch / ir->rcoulomb < 0.9499)
1371 real percentage = 100 * (ir->rcoulomb - ir->rcoulomb_switch) / ir->rcoulomb;
1373 "The switching range should be 5%% or less (currently %.2f%% using a switching "
1374 "range of %4f-%4f) for accurate electrostatic energies, energy conservation "
1375 "will be good regardless, since ewald_rtol = %g.",
1376 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1377 warning(wi, warn_buf);
1381 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1383 if (ir->rvdw_switch == 0)
1386 "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones "
1387 "potential. This suggests it was not set in the mdp, which can lead to large "
1388 "energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw "
1389 "switching range.");
1390 warning(wi, warn_buf);
1394 if (EEL_FULL(ir->coulombtype))
1396 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER
1397 || ir->coulombtype == eelPMEUSERSWITCH)
1399 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1400 eel_names[ir->coulombtype]);
1401 CHECK(ir->rcoulomb > ir->rlist);
1405 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1407 // TODO: Move these checks into the ewald module with the options class
1409 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1411 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1413 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d",
1414 eel_names[ir->coulombtype], orderMin, orderMax);
1415 warning_error(wi, warn_buf);
1419 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1421 if (ir->ewald_geometry == eewg3D)
1423 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1424 c_pbcTypeNames[ir->pbcType].c_str(), eewg_names[eewg3DC]);
1425 warning(wi, warn_buf);
1427 /* This check avoids extra pbc coding for exclusion corrections */
1428 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1429 CHECK(ir->wall_ewald_zfac < 2);
1431 if ((ir->ewald_geometry == eewg3DC) && (ir->pbcType != PbcType::XY) && EEL_FULL(ir->coulombtype))
1433 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1434 eel_names[ir->coulombtype], eewg_names[eewg3DC], c_pbcTypeNames[PbcType::XY].c_str());
1435 warning(wi, warn_buf);
1437 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1439 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1440 CHECK(ir->bPeriodicMols);
1441 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1442 warning_note(wi, warn_buf);
1444 "With epsilon_surface > 0 you can only use domain decomposition "
1445 "when there are only small molecules with all bonds constrained (mdrun will check "
1447 warning_note(wi, warn_buf);
1450 if (ir_vdw_switched(ir))
1452 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1453 CHECK(ir->rvdw_switch >= ir->rvdw);
1455 if (ir->rvdw_switch < 0.5 * ir->rvdw)
1458 "You are applying a switch function to vdw forces or potentials from %g to %g "
1459 "nm, which is more than half the interaction range, whereas switch functions "
1460 "are intended to act only close to the cut-off.",
1461 ir->rvdw_switch, ir->rvdw);
1462 warning_note(wi, warn_buf);
1466 if (ir->vdwtype == evdwPME)
1468 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1470 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1471 evdw_names[ir->vdwtype], eintmod_names[eintmodPOTSHIFT], eintmod_names[eintmodNONE]);
1472 warning_error(wi, err_buf);
1476 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1479 "You have selected user tables with dispersion correction, the dispersion "
1480 "will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction "
1481 "between rvdw_switch and rvdw will not be double counted). Make sure that you "
1482 "really want dispersion correction to -C6/r^6.");
1485 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT) && ir->rvdw != 0)
1487 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1490 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1492 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1495 /* IMPLICIT SOLVENT */
1496 if (ir->coulombtype == eelGB_NOTUSED)
1498 sprintf(warn_buf, "Invalid option %s for coulombtype", eel_names[ir->coulombtype]);
1499 warning_error(wi, warn_buf);
1504 warning_error(wi, "The QMMM integration you are trying to use is no longer supported");
1509 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1512 // cosine acceleration is only supported in leap-frog
1513 if (ir->cos_accel != 0.0 && ir->eI != eiMD)
1515 warning_error(wi, "cos-acceleration is only supported by integrator = md");
1519 /* interpret a number of doubles from a string and put them in an array,
1520 after allocating space for them.
1521 str = the input string
1522 n = the (pre-allocated) number of doubles read
1523 r = the output array of doubles. */
1524 static void parse_n_real(char* str, int* n, real** r, warninp_t wi)
1526 auto values = gmx::splitString(str);
1530 for (int i = 0; i < *n; i++)
1534 (*r)[i] = gmx::fromString<real>(values[i]);
1536 catch (gmx::GromacsException&)
1538 warning_error(wi, "Invalid value " + values[i]
1539 + " in string in mdp file. Expected a real number.");
1545 static void do_fep_params(t_inputrec* ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1548 int i, j, max_n_lambda, nweights, nfep[efptNR];
1549 t_lambda* fep = ir->fepvals;
1550 t_expanded* expand = ir->expandedvals;
1551 real** count_fep_lambdas;
1552 bool bOneLambda = TRUE;
1554 snew(count_fep_lambdas, efptNR);
1556 /* FEP input processing */
1557 /* first, identify the number of lambda values for each type.
1558 All that are nonzero must have the same number */
1560 for (i = 0; i < efptNR; i++)
1562 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1565 /* now, determine the number of components. All must be either zero, or equal. */
1568 for (i = 0; i < efptNR; i++)
1570 if (nfep[i] > max_n_lambda)
1572 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1573 must have the same number if its not zero.*/
1578 for (i = 0; i < efptNR; i++)
1582 ir->fepvals->separate_dvdl[i] = FALSE;
1584 else if (nfep[i] == max_n_lambda)
1586 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute
1587 the derivative with respect to the temperature currently */
1589 ir->fepvals->separate_dvdl[i] = TRUE;
1595 "Number of lambdas (%d) for FEP type %s not equal to number of other types "
1597 nfep[i], efpt_names[i], max_n_lambda);
1600 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1601 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1603 /* the number of lambdas is the number we've read in, which is either zero
1604 or the same for all */
1605 fep->n_lambda = max_n_lambda;
1607 /* allocate space for the array of lambda values */
1608 snew(fep->all_lambda, efptNR);
1609 /* if init_lambda is defined, we need to set lambda */
1610 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1612 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1614 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1615 for (i = 0; i < efptNR; i++)
1617 snew(fep->all_lambda[i], fep->n_lambda);
1618 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1621 for (j = 0; j < fep->n_lambda; j++)
1623 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1625 sfree(count_fep_lambdas[i]);
1628 sfree(count_fep_lambdas);
1630 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for
1631 internal bookkeeping -- for now, init_lambda */
1633 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1635 for (i = 0; i < fep->n_lambda; i++)
1637 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1641 /* check to see if only a single component lambda is defined, and soft core is defined.
1642 In this case, turn on coulomb soft core */
1644 if (max_n_lambda == 0)
1650 for (i = 0; i < efptNR; i++)
1652 if ((nfep[i] != 0) && (i != efptFEP))
1658 if ((bOneLambda) && (fep->sc_alpha > 0))
1660 fep->bScCoul = TRUE;
1663 /* Fill in the others with the efptFEP if they are not explicitly
1664 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1665 they are all zero. */
1667 for (i = 0; i < efptNR; i++)
1669 if ((nfep[i] == 0) && (i != efptFEP))
1671 for (j = 0; j < fep->n_lambda; j++)
1673 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1679 /* now read in the weights */
1680 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1683 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1685 else if (nweights != fep->n_lambda)
1687 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1688 nweights, fep->n_lambda);
1690 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1692 expand->nstexpanded = fep->nstdhdl;
1693 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1698 static void do_simtemp_params(t_inputrec* ir)
1701 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1702 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1705 template<typename T>
1706 void convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, T* outputs)
1709 for (const auto& input : inputs)
1713 outputs[i] = gmx::fromStdString<T>(input);
1715 catch (gmx::GromacsException&)
1717 auto message = gmx::formatString(
1718 "Invalid value for mdp option %s. %s should only consist of integers separated "
1721 warning_error(wi, message);
1727 static void convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, real* outputs)
1730 for (const auto& input : inputs)
1734 outputs[i] = gmx::fromString<real>(input);
1736 catch (gmx::GromacsException&)
1738 auto message = gmx::formatString(
1739 "Invalid value for mdp option %s. %s should only consist of real numbers "
1740 "separated by spaces.",
1742 warning_error(wi, message);
1748 static void convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char* name, rvec* outputs)
1751 for (const auto& input : inputs)
1755 outputs[i][d] = gmx::fromString<real>(input);
1757 catch (gmx::GromacsException&)
1759 auto message = gmx::formatString(
1760 "Invalid value for mdp option %s. %s should only consist of real numbers "
1761 "separated by spaces.",
1763 warning_error(wi, message);
1774 static void do_wall_params(t_inputrec* ir, char* wall_atomtype, char* wall_density, t_gromppopts* opts, warninp_t wi)
1776 opts->wall_atomtype[0] = nullptr;
1777 opts->wall_atomtype[1] = nullptr;
1779 ir->wall_atomtype[0] = -1;
1780 ir->wall_atomtype[1] = -1;
1781 ir->wall_density[0] = 0;
1782 ir->wall_density[1] = 0;
1786 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1787 if (wallAtomTypes.size() != size_t(ir->nwall))
1789 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu", ir->nwall,
1790 wallAtomTypes.size());
1792 GMX_RELEASE_ASSERT(ir->nwall < 3, "Invalid number of walls");
1793 for (int i = 0; i < ir->nwall; i++)
1795 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1798 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1800 auto wallDensity = gmx::splitString(wall_density);
1801 if (wallDensity.size() != size_t(ir->nwall))
1803 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall,
1804 wallDensity.size());
1806 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1807 for (int i = 0; i < ir->nwall; i++)
1809 if (ir->wall_density[i] <= 0)
1811 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1818 static void add_wall_energrps(SimulationGroups* groups, int nwall, t_symtab* symtab)
1822 AtomGroupIndices* grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1823 for (int i = 0; i < nwall; i++)
1825 groups->groupNames.emplace_back(put_symtab(symtab, gmx::formatString("wall%d", i).c_str()));
1826 grps->emplace_back(groups->groupNames.size() - 1);
1831 static void read_expandedparams(std::vector<t_inpfile>* inp, t_expanded* expand, warninp_t wi)
1833 /* read expanded ensemble parameters */
1834 printStringNewline(inp, "expanded ensemble variables");
1835 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1836 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1837 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1838 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1839 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1840 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1841 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1842 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1843 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1844 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1845 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1846 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1847 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1848 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1849 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1850 expand->bSymmetrizedTMatrix =
1851 (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1852 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1853 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1854 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1855 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1856 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1857 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1858 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1861 /*! \brief Return whether an end state with the given coupling-lambda
1862 * value describes fully-interacting VDW.
1864 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1865 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1867 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1869 return (couple_lambda_value == ecouplamVDW || couple_lambda_value == ecouplamVDWQ);
1875 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1878 explicit MdpErrorHandler(warninp_t wi) : wi_(wi), mapping_(nullptr) {}
1880 void setBackMapping(const gmx::IKeyValueTreeBackMapping& mapping) { mapping_ = &mapping; }
1882 bool onError(gmx::UserInputError* ex, const gmx::KeyValueTreePath& context) override
1885 gmx::formatString("Error in mdp option \"%s\":", getOptionName(context).c_str()));
1886 std::string message = gmx::formatExceptionMessageToString(*ex);
1887 warning_error(wi_, message.c_str());
1892 std::string getOptionName(const gmx::KeyValueTreePath& context)
1894 if (mapping_ != nullptr)
1896 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1897 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1900 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1905 const gmx::IKeyValueTreeBackMapping* mapping_;
1910 void get_ir(const char* mdparin,
1911 const char* mdparout,
1912 gmx::MDModules* mdModules,
1915 WriteMdpHeader writeMdpHeader,
1919 double dumdub[2][6];
1921 char warn_buf[STRLEN];
1922 t_lambda* fep = ir->fepvals;
1923 t_expanded* expand = ir->expandedvals;
1925 const char* no_names[] = { "no", nullptr };
1927 init_inputrec_strings();
1928 gmx::TextInputFile stream(mdparin);
1929 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1931 snew(dumstr[0], STRLEN);
1932 snew(dumstr[1], STRLEN);
1934 /* ignore the following deprecated commands */
1935 replace_inp_entry(inp, "title", nullptr);
1936 replace_inp_entry(inp, "cpp", nullptr);
1937 replace_inp_entry(inp, "domain-decomposition", nullptr);
1938 replace_inp_entry(inp, "andersen-seed", nullptr);
1939 replace_inp_entry(inp, "dihre", nullptr);
1940 replace_inp_entry(inp, "dihre-fc", nullptr);
1941 replace_inp_entry(inp, "dihre-tau", nullptr);
1942 replace_inp_entry(inp, "nstdihreout", nullptr);
1943 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1944 replace_inp_entry(inp, "optimize-fft", nullptr);
1945 replace_inp_entry(inp, "adress_type", nullptr);
1946 replace_inp_entry(inp, "adress_const_wf", nullptr);
1947 replace_inp_entry(inp, "adress_ex_width", nullptr);
1948 replace_inp_entry(inp, "adress_hy_width", nullptr);
1949 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1950 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1951 replace_inp_entry(inp, "adress_site", nullptr);
1952 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1953 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1954 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1955 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1956 replace_inp_entry(inp, "rlistlong", nullptr);
1957 replace_inp_entry(inp, "nstcalclr", nullptr);
1958 replace_inp_entry(inp, "pull-print-com2", nullptr);
1959 replace_inp_entry(inp, "gb-algorithm", nullptr);
1960 replace_inp_entry(inp, "nstgbradii", nullptr);
1961 replace_inp_entry(inp, "rgbradii", nullptr);
1962 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1963 replace_inp_entry(inp, "gb-saltconc", nullptr);
1964 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1965 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1966 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1967 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1968 replace_inp_entry(inp, "sa-algorithm", nullptr);
1969 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1970 replace_inp_entry(inp, "ns-type", nullptr);
1972 /* replace the following commands with the clearer new versions*/
1973 replace_inp_entry(inp, "unconstrained-start", "continuation");
1974 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1975 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1976 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1977 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1978 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1979 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1981 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1982 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1983 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1984 setStringEntry(&inp, "include", opts->include, nullptr);
1985 printStringNoNewline(
1986 &inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1987 setStringEntry(&inp, "define", opts->define, nullptr);
1989 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1990 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1991 printStringNoNewline(&inp, "Start time and timestep in ps");
1992 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1993 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1994 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1995 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1996 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1997 printStringNoNewline(
1998 &inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1999 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
2000 printStringNoNewline(&inp, "Multiple time-stepping");
2001 ir->useMts = (get_eeenum(&inp, "mts", yesno_names, wi) != 0);
2004 opts->numMtsLevels = get_eint(&inp, "mts-levels", 2, wi);
2005 ir->mtsLevels.resize(2);
2006 gmx::MtsLevel& mtsLevel = ir->mtsLevels[1];
2007 opts->mtsLevel2Forces = setStringEntry(&inp, "mts-level2-forces", "longrange-nonbonded");
2008 mtsLevel.stepFactor = 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))
2014 ir->mtsLevels.clear();
2017 printStringNoNewline(&inp, "mode for center of mass motion removal");
2018 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
2019 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
2020 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
2021 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
2022 setStringEntry(&inp, "comm-grps", inputrecStrings->vcm, nullptr);
2024 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
2025 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
2026 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
2027 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
2030 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
2031 printStringNoNewline(&inp, "Force tolerance and initial step-size");
2032 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
2033 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
2034 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
2035 ir->niter = get_eint(&inp, "niter", 20, wi);
2036 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
2037 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
2038 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
2039 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
2040 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
2042 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
2043 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
2045 /* Output options */
2046 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
2047 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
2048 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
2049 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
2050 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
2051 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
2052 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
2053 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
2054 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
2055 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
2056 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
2057 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
2058 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
2059 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
2060 printStringNoNewline(&inp, "default, all atoms will be written.");
2061 setStringEntry(&inp, "compressed-x-grps", inputrecStrings->x_compressed_groups, nullptr);
2062 printStringNoNewline(&inp, "Selection of energy groups");
2063 setStringEntry(&inp, "energygrps", inputrecStrings->energy, nullptr);
2065 /* Neighbor searching */
2066 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
2067 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs)");
2068 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
2069 printStringNoNewline(&inp, "nblist update frequency");
2070 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
2071 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
2072 // TODO This conversion should be removed when proper std:string handling will be added to get_eeenum(...), etc.
2073 std::vector<const char*> pbcTypesNamesChar;
2074 for (const auto& pbcTypeName : c_pbcTypeNames)
2076 pbcTypesNamesChar.push_back(pbcTypeName.c_str());
2078 ir->pbcType = static_cast<PbcType>(get_eeenum(&inp, "pbc", pbcTypesNamesChar.data(), wi));
2079 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
2080 printStringNoNewline(&inp,
2081 "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
2082 printStringNoNewline(&inp, "a value of -1 means: use rlist");
2083 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
2084 printStringNoNewline(&inp, "nblist cut-off");
2085 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
2086 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
2088 /* Electrostatics */
2089 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
2090 printStringNoNewline(&inp, "Method for doing electrostatics");
2091 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
2092 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
2093 printStringNoNewline(&inp, "cut-off lengths");
2094 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
2095 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
2096 printStringNoNewline(&inp,
2097 "Relative dielectric constant for the medium and the reaction field");
2098 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
2099 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
2100 printStringNoNewline(&inp, "Method for doing Van der Waals");
2101 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
2102 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
2103 printStringNoNewline(&inp, "cut-off lengths");
2104 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
2105 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
2106 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
2107 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
2108 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
2109 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
2110 printStringNoNewline(&inp, "Separate tables between energy group pairs");
2111 setStringEntry(&inp, "energygrp-table", inputrecStrings->egptable, nullptr);
2112 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
2113 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
2114 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
2115 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
2116 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
2117 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
2118 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
2119 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
2120 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
2121 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
2122 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
2123 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
2124 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
2126 /* Implicit solvation is no longer supported, but we need grompp
2127 to be able to refuse old .mdp files that would have built a tpr
2128 to run it. Thus, only "no" is accepted. */
2129 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
2131 /* Coupling stuff */
2132 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
2133 printStringNoNewline(&inp, "Temperature coupling");
2134 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
2135 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
2136 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
2137 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
2138 printStringNoNewline(&inp, "Groups to couple separately");
2139 setStringEntry(&inp, "tc-grps", inputrecStrings->tcgrps, nullptr);
2140 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
2141 setStringEntry(&inp, "tau-t", inputrecStrings->tau_t, nullptr);
2142 setStringEntry(&inp, "ref-t", inputrecStrings->ref_t, nullptr);
2143 printStringNoNewline(&inp, "pressure coupling");
2144 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
2145 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
2146 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
2147 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
2148 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
2149 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
2150 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
2151 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
2152 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
2155 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
2156 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
2157 printStringNoNewline(&inp, "Groups treated with MiMiC");
2158 setStringEntry(&inp, "QMMM-grps", inputrecStrings->QMMM, nullptr);
2160 /* Simulated annealing */
2161 printStringNewline(&inp, "SIMULATED ANNEALING");
2162 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2163 setStringEntry(&inp, "annealing", inputrecStrings->anneal, nullptr);
2164 printStringNoNewline(&inp,
2165 "Number of time points to use for specifying annealing in each group");
2166 setStringEntry(&inp, "annealing-npoints", inputrecStrings->anneal_npoints, nullptr);
2167 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2168 setStringEntry(&inp, "annealing-time", inputrecStrings->anneal_time, nullptr);
2169 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2170 setStringEntry(&inp, "annealing-temp", inputrecStrings->anneal_temp, nullptr);
2173 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2174 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
2175 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2176 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2179 printStringNewline(&inp, "OPTIONS FOR BONDS");
2180 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2181 printStringNoNewline(&inp, "Type of constraint algorithm");
2182 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
2183 printStringNoNewline(&inp, "Do not constrain the start configuration");
2184 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
2185 printStringNoNewline(&inp,
2186 "Use successive overrelaxation to reduce the number of shake iterations");
2187 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
2188 printStringNoNewline(&inp, "Relative tolerance of shake");
2189 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2190 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2191 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2192 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2193 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2194 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2195 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2196 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2197 printStringNoNewline(&inp, "rotates over more degrees than");
2198 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2199 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2200 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
2202 /* Energy group exclusions */
2203 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2204 printStringNoNewline(
2205 &inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2206 setStringEntry(&inp, "energygrp-excl", inputrecStrings->egpexcl, nullptr);
2209 printStringNewline(&inp, "WALLS");
2210 printStringNoNewline(
2211 &inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2212 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2213 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
2214 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2215 setStringEntry(&inp, "wall-atomtype", inputrecStrings->wall_atomtype, nullptr);
2216 setStringEntry(&inp, "wall-density", inputrecStrings->wall_density, nullptr);
2217 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2220 printStringNewline(&inp, "COM PULLING");
2221 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2224 ir->pull = std::make_unique<pull_params_t>();
2225 inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
2229 for (int c = 0; c < ir->pull->ncoord; c++)
2231 if (ir->pull->coord[c].eType == epullCONSTRAINT)
2234 "Constraint COM pulling is not supported in combination with "
2235 "multiple time stepping");
2243 NOTE: needs COM pulling or free energy input */
2244 printStringNewline(&inp, "AWH biasing");
2245 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2248 ir->awhParams = gmx::readAwhParams(&inp, wi);
2251 /* Enforced rotation */
2252 printStringNewline(&inp, "ENFORCED ROTATION");
2253 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2254 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2258 inputrecStrings->rotateGroupNames = read_rotparams(&inp, ir->rot, wi);
2261 /* Interactive MD */
2263 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2264 setStringEntry(&inp, "IMD-group", inputrecStrings->imd_grp, nullptr);
2265 if (inputrecStrings->imd_grp[0] != '\0')
2272 printStringNewline(&inp, "NMR refinement stuff");
2273 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2274 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2275 printStringNoNewline(
2276 &inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2277 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2278 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2279 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2280 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2281 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2282 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2283 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2284 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2285 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2286 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2287 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2288 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2289 setStringEntry(&inp, "orire-fitgrp", inputrecStrings->orirefitgrp, nullptr);
2290 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2291 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2293 /* free energy variables */
2294 printStringNewline(&inp, "Free energy variables");
2295 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2296 setStringEntry(&inp, "couple-moltype", inputrecStrings->couple_moltype, nullptr);
2297 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2298 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2299 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2301 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2303 it was not entered */
2304 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2305 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2306 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2307 setStringEntry(&inp, "fep-lambdas", inputrecStrings->fep_lambda[efptFEP], nullptr);
2308 setStringEntry(&inp, "mass-lambdas", inputrecStrings->fep_lambda[efptMASS], nullptr);
2309 setStringEntry(&inp, "coul-lambdas", inputrecStrings->fep_lambda[efptCOUL], nullptr);
2310 setStringEntry(&inp, "vdw-lambdas", inputrecStrings->fep_lambda[efptVDW], nullptr);
2311 setStringEntry(&inp, "bonded-lambdas", inputrecStrings->fep_lambda[efptBONDED], nullptr);
2312 setStringEntry(&inp, "restraint-lambdas", inputrecStrings->fep_lambda[efptRESTRAINT], nullptr);
2313 setStringEntry(&inp, "temperature-lambdas", inputrecStrings->fep_lambda[efptTEMPERATURE], nullptr);
2314 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2315 setStringEntry(&inp, "init-lambda-weights", inputrecStrings->lambda_weights, nullptr);
2316 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2317 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2318 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2319 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2320 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2321 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2322 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2323 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2324 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2325 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2326 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2327 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2329 /* Non-equilibrium MD stuff */
2330 printStringNewline(&inp, "Non-equilibrium MD stuff");
2331 setStringEntry(&inp, "acc-grps", inputrecStrings->accgrps, nullptr);
2332 setStringEntry(&inp, "accelerate", inputrecStrings->acc, nullptr);
2333 setStringEntry(&inp, "freezegrps", inputrecStrings->freeze, nullptr);
2334 setStringEntry(&inp, "freezedim", inputrecStrings->frdim, nullptr);
2335 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2336 setStringEntry(&inp, "deform", inputrecStrings->deform, nullptr);
2338 /* simulated tempering variables */
2339 printStringNewline(&inp, "simulated tempering variables");
2340 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2341 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2342 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2343 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2345 /* expanded ensemble variables */
2346 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2348 read_expandedparams(&inp, expand, wi);
2351 /* Electric fields */
2353 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2354 gmx::KeyValueTreeTransformer transform;
2355 transform.rules()->addRule().keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2356 mdModules->initMdpTransform(transform.rules());
2357 for (const auto& path : transform.mappedPaths())
2359 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2360 mark_einp_set(inp, path[0].c_str());
2362 MdpErrorHandler errorHandler(wi);
2363 auto result = transform.transform(convertedValues, &errorHandler);
2364 ir->params = new gmx::KeyValueTreeObject(result.object());
2365 mdModules->adjustInputrecBasedOnModules(ir);
2366 errorHandler.setBackMapping(result.backMapping());
2367 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2370 /* Ion/water position swapping ("computational electrophysiology") */
2371 printStringNewline(&inp,
2372 "Ion/water position swapping for computational electrophysiology setups");
2373 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2374 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2375 if (ir->eSwapCoords != eswapNO)
2382 printStringNoNewline(&inp, "Swap attempt frequency");
2383 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2384 printStringNoNewline(&inp, "Number of ion types to be controlled");
2385 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2388 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2390 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2391 snew(ir->swap->grp, ir->swap->ngrp);
2392 for (i = 0; i < ir->swap->ngrp; i++)
2394 snew(ir->swap->grp[i].molname, STRLEN);
2396 printStringNoNewline(&inp,
2397 "Two index groups that contain the compartment-partitioning atoms");
2398 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2399 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2400 printStringNoNewline(&inp,
2401 "Use center of mass of split groups (yes/no), otherwise center of "
2402 "geometry is used");
2403 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2404 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2406 printStringNoNewline(&inp, "Name of solvent molecules");
2407 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2409 printStringNoNewline(&inp,
2410 "Split cylinder: radius, upper and lower extension (nm) (this will "
2411 "define the channels)");
2412 printStringNoNewline(&inp,
2413 "Note that the split cylinder settings do not have an influence on "
2414 "the swapping protocol,");
2415 printStringNoNewline(
2417 "however, if correctly defined, the permeation events are recorded per channel");
2418 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2419 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2420 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2421 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2422 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2423 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2425 printStringNoNewline(
2427 "Average the number of ions per compartment over these many swap attempt steps");
2428 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2430 printStringNoNewline(
2431 &inp, "Names of the ion types that can be exchanged with solvent molecules,");
2432 printStringNoNewline(
2433 &inp, "and the requested number of ions of this type in compartments A and B");
2434 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2435 for (i = 0; i < nIonTypes; i++)
2437 int ig = eSwapFixedGrpNR + i;
2439 sprintf(buf, "iontype%d-name", i);
2440 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2441 sprintf(buf, "iontype%d-in-A", i);
2442 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2443 sprintf(buf, "iontype%d-in-B", i);
2444 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2447 printStringNoNewline(
2449 "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2450 printStringNoNewline(
2452 "at maximum distance (= bulk concentration) to the split group layers. However,");
2453 printStringNoNewline(&inp,
2454 "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk "
2455 "layer from the middle at 0.0");
2456 printStringNoNewline(&inp,
2457 "towards one of the compartment-partitioning layers (at +/- 1.0).");
2458 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2459 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2460 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2461 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0))
2463 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2466 printStringNoNewline(
2467 &inp, "Start to swap ions if threshold difference to requested count is reached");
2468 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2471 /* AdResS is no longer supported, but we need grompp to be able to
2472 refuse to process old .mdp files that used it. */
2473 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2475 /* User defined thingies */
2476 printStringNewline(&inp, "User defined thingies");
2477 setStringEntry(&inp, "user1-grps", inputrecStrings->user1, nullptr);
2478 setStringEntry(&inp, "user2-grps", inputrecStrings->user2, nullptr);
2479 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2480 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2481 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2482 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2483 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2484 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2485 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2486 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2490 gmx::TextOutputFile stream(mdparout);
2491 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2493 // Transform module data into a flat key-value tree for output.
2494 gmx::KeyValueTreeBuilder builder;
2495 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2496 mdModules->buildMdpOutput(&builderObject);
2498 gmx::TextWriter writer(&stream);
2499 writeKeyValueTreeAsMdp(&writer, builder.build());
2504 /* Process options if necessary */
2505 for (m = 0; m < 2; m++)
2507 for (i = 0; i < 2 * DIM; i++)
2516 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2520 "Pressure coupling incorrect number of values (I need exactly 1)");
2522 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2524 case epctSEMIISOTROPIC:
2525 case epctSURFACETENSION:
2526 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2530 "Pressure coupling incorrect number of values (I need exactly 2)");
2532 dumdub[m][YY] = dumdub[m][XX];
2534 case epctANISOTROPIC:
2535 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf", &(dumdub[m][XX]), &(dumdub[m][YY]),
2536 &(dumdub[m][ZZ]), &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5]))
2541 "Pressure coupling incorrect number of values (I need exactly 6)");
2545 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2546 epcoupltype_names[ir->epct]);
2550 clear_mat(ir->ref_p);
2551 clear_mat(ir->compress);
2552 for (i = 0; i < DIM; i++)
2554 ir->ref_p[i][i] = dumdub[1][i];
2555 ir->compress[i][i] = dumdub[0][i];
2557 if (ir->epct == epctANISOTROPIC)
2559 ir->ref_p[XX][YY] = dumdub[1][3];
2560 ir->ref_p[XX][ZZ] = dumdub[1][4];
2561 ir->ref_p[YY][ZZ] = dumdub[1][5];
2562 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2565 "All off-diagonal reference pressures are non-zero. Are you sure you want to "
2566 "apply a threefold shear stress?\n");
2568 ir->compress[XX][YY] = dumdub[0][3];
2569 ir->compress[XX][ZZ] = dumdub[0][4];
2570 ir->compress[YY][ZZ] = dumdub[0][5];
2571 for (i = 0; i < DIM; i++)
2573 for (m = 0; m < i; m++)
2575 ir->ref_p[i][m] = ir->ref_p[m][i];
2576 ir->compress[i][m] = ir->compress[m][i];
2581 if (ir->comm_mode == ecmNO)
2586 opts->couple_moltype = nullptr;
2587 if (strlen(inputrecStrings->couple_moltype) > 0)
2589 if (ir->efep != efepNO)
2591 opts->couple_moltype = gmx_strdup(inputrecStrings->couple_moltype);
2592 if (opts->couple_lam0 == opts->couple_lam1)
2594 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2596 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam1 == ecouplamNONE))
2600 "For proper sampling of the (nearly) decoupled state, stochastic dynamics "
2607 "Free energy is turned off, so we will not decouple the molecule listed "
2611 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2612 if (ir->efep != efepNO)
2614 if (fep->delta_lambda != 0)
2616 ir->efep = efepSLOWGROWTH;
2620 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2622 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2624 "Old option for dhdl-print-energy given: "
2625 "changing \"yes\" to \"total\"\n");
2628 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2630 /* always print out the energy to dhdl if we are doing
2631 expanded ensemble, since we need the total energy for
2632 analysis if the temperature is changing. In some
2633 conditions one may only want the potential energy, so
2634 we will allow that if the appropriate mdp setting has
2635 been enabled. Otherwise, total it is:
2637 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2640 if ((ir->efep != efepNO) || ir->bSimTemp)
2642 ir->bExpanded = FALSE;
2643 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2645 ir->bExpanded = TRUE;
2647 do_fep_params(ir, inputrecStrings->fep_lambda, inputrecStrings->lambda_weights, wi);
2648 if (ir->bSimTemp) /* done after fep params */
2650 do_simtemp_params(ir);
2653 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2654 * setup and not on the old way of specifying the free-energy setup,
2655 * we should check for using soft-core when not needed, since that
2656 * can complicate the sampling significantly.
2657 * Note that we only check for the automated coupling setup.
2658 * If the (advanced) user does FEP through manual topology changes,
2659 * this check will not be triggered.
2661 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 && ir->fepvals->sc_alpha != 0
2662 && (couple_lambda_has_vdw_on(opts->couple_lam0) && couple_lambda_has_vdw_on(opts->couple_lam1)))
2665 "You are using soft-core interactions while the Van der Waals interactions are "
2666 "not decoupled (note that the sc-coul option is only active when using lambda "
2667 "states). Although this will not lead to errors, you will need much more "
2668 "sampling than without soft-core interactions. Consider using sc-alpha=0.");
2673 ir->fepvals->n_lambda = 0;
2676 /* WALL PARAMETERS */
2678 do_wall_params(ir, inputrecStrings->wall_atomtype, inputrecStrings->wall_density, opts, wi);
2680 /* ORIENTATION RESTRAINT PARAMETERS */
2682 if (opts->bOrire && gmx::splitString(inputrecStrings->orirefitgrp).size() != 1)
2684 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2687 /* DEFORMATION PARAMETERS */
2689 clear_mat(ir->deform);
2690 for (i = 0; i < 6; i++)
2695 double gmx_unused canary;
2696 int ndeform = sscanf(inputrecStrings->deform, "%lf %lf %lf %lf %lf %lf %lf", &(dumdub[0][0]),
2697 &(dumdub[0][1]), &(dumdub[0][2]), &(dumdub[0][3]), &(dumdub[0][4]),
2698 &(dumdub[0][5]), &canary);
2700 if (strlen(inputrecStrings->deform) > 0 && ndeform != 6)
2704 "Cannot parse exactly 6 box deformation velocities from string '%s'",
2705 inputrecStrings->deform)
2708 for (i = 0; i < 3; i++)
2710 ir->deform[i][i] = dumdub[0][i];
2712 ir->deform[YY][XX] = dumdub[0][3];
2713 ir->deform[ZZ][XX] = dumdub[0][4];
2714 ir->deform[ZZ][YY] = dumdub[0][5];
2715 if (ir->epc != epcNO)
2717 for (i = 0; i < 3; i++)
2719 for (j = 0; j <= i; j++)
2721 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2723 warning_error(wi, "A box element has deform set and compressibility > 0");
2727 for (i = 0; i < 3; i++)
2729 for (j = 0; j < i; j++)
2731 if (ir->deform[i][j] != 0)
2733 for (m = j; m < DIM; m++)
2735 if (ir->compress[m][j] != 0)
2738 "An off-diagonal box element has deform set while "
2739 "compressibility > 0 for the same component of another box "
2740 "vector, this might lead to spurious periodicity effects.");
2741 warning(wi, warn_buf);
2749 /* Ion/water position swapping checks */
2750 if (ir->eSwapCoords != eswapNO)
2752 if (ir->swap->nstswap < 1)
2754 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2756 if (ir->swap->nAverage < 1)
2758 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2760 if (ir->swap->threshold < 1.0)
2762 warning_error(wi, "Ion count threshold must be at least 1.\n");
2766 /* Set up MTS levels, this needs to happen before checking AWH parameters */
2769 setupMtsLevels(ir->mtsLevels, *ir, *opts, wi);
2774 gmx::checkAwhParams(ir->awhParams, ir, wi);
2781 /* We would like gn to be const as well, but C doesn't allow this */
2782 /* TODO this is utility functionality (search for the index of a
2783 string in a collection), so should be refactored and located more
2785 int search_string(const char* s, int ng, char* gn[])
2789 for (i = 0; (i < ng); i++)
2791 if (gmx_strcasecmp(s, gn[i]) == 0)
2798 "Group %s referenced in the .mdp file was not found in the index file.\n"
2799 "Group names must match either [moleculetype] names or custom index group\n"
2800 "names, in which case you must supply an index file to the '-n' option\n"
2805 static void atomGroupRangeValidation(int natoms, int groupIndex, const t_blocka& block)
2807 /* Now go over the atoms in the group */
2808 for (int j = block.index[groupIndex]; (j < block.index[groupIndex + 1]); j++)
2810 int aj = block.a[j];
2812 /* Range checking */
2813 if ((aj < 0) || (aj >= natoms))
2815 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2820 static void do_numbering(int natoms,
2821 SimulationGroups* groups,
2822 gmx::ArrayRef<std::string> groupsFromMdpFile,
2825 SimulationAtomGroupType gtype,
2831 unsigned short* cbuf;
2832 AtomGroupIndices* grps = &(groups->groups[gtype]);
2835 char warn_buf[STRLEN];
2837 title = shortName(gtype);
2840 /* Mark all id's as not set */
2841 for (int i = 0; (i < natoms); i++)
2846 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2848 /* Lookup the group name in the block structure */
2849 const int gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2850 if ((grptp != egrptpONE) || (i == 0))
2852 grps->emplace_back(gid);
2854 GMX_ASSERT(block, "Can't have a nullptr block");
2855 atomGroupRangeValidation(natoms, gid, *block);
2856 /* Now go over the atoms in the group */
2857 for (int j = block->index[gid]; (j < block->index[gid + 1]); j++)
2859 const int aj = block->a[j];
2860 /* Lookup up the old group number */
2861 const int ognr = cbuf[aj];
2864 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)", aj + 1, title,
2869 /* Store the group number in buffer */
2870 if (grptp == egrptpONE)
2883 /* Now check whether we have done all atoms */
2886 if (grptp == egrptpALL)
2888 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2890 else if (grptp == egrptpPART)
2892 sprintf(warn_buf, "%d atoms are not part of any of the %s groups", natoms - ntot, title);
2893 warning_note(wi, warn_buf);
2895 /* Assign all atoms currently unassigned to a rest group */
2896 for (int j = 0; (j < natoms); j++)
2898 if (cbuf[j] == NOGID)
2900 cbuf[j] = grps->size();
2903 if (grptp != egrptpPART)
2907 fprintf(stderr, "Making dummy/rest group for %s containing %d elements\n", title,
2910 /* Add group name "rest" */
2911 grps->emplace_back(restnm);
2913 /* Assign the rest name to all atoms not currently assigned to a group */
2914 for (int j = 0; (j < natoms); j++)
2916 if (cbuf[j] == NOGID)
2918 // group size was not updated before this here, so need to use -1.
2919 cbuf[j] = grps->size() - 1;
2925 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2927 /* All atoms are part of one (or no) group, no index required */
2928 groups->groupNumbers[gtype].clear();
2932 for (int j = 0; (j < natoms); j++)
2934 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2941 static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
2944 pull_params_t* pull;
2945 int natoms, imin, jmin;
2946 int * nrdf2, *na_vcm, na_tot;
2947 double * nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2952 * First calc 3xnr-atoms for each group
2953 * then subtract half a degree of freedom for each constraint
2955 * Only atoms and nuclei contribute to the degrees of freedom...
2960 const SimulationGroups& groups = mtop->groups;
2961 natoms = mtop->natoms;
2963 /* Allocate one more for a possible rest group */
2964 /* We need to sum degrees of freedom into doubles,
2965 * since floats give too low nrdf's above 3 million atoms.
2967 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size() + 1);
2968 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2969 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2970 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2971 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size() + 1);
2973 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2977 for (gmx::index i = 0;
2978 i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; i++)
2981 clear_ivec(dof_vcm[i]);
2983 nrdf_vcm_sub[i] = 0;
2985 snew(nrdf2, natoms);
2986 for (const AtomProxy atomP : AtomRange(*mtop))
2988 const t_atom& local = atomP.atom();
2989 int i = atomP.globalAtomNumber();
2991 if (local.ptype == eptAtom || local.ptype == eptNucleus)
2993 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
2994 for (int d = 0; d < DIM; d++)
2996 if (opts->nFreeze[g][d] == 0)
2998 /* Add one DOF for particle i (counted as 2*1) */
3000 /* VCM group i has dim d as a DOF */
3001 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] =
3005 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] +=
3007 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] +=
3013 for (const gmx_molblock_t& molb : mtop->molblock)
3015 const gmx_moltype_t& molt = mtop->moltype[molb.type];
3016 const t_atom* atom = molt.atoms.atom;
3017 for (int mol = 0; mol < molb.nmol; mol++)
3019 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
3021 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
3022 for (int i = 0; i < molt.ilist[ftype].size();)
3024 /* Subtract degrees of freedom for the constraints,
3025 * if the particles still have degrees of freedom left.
3026 * If one of the particles is a vsite or a shell, then all
3027 * constraint motion will go there, but since they do not
3028 * contribute to the constraints the degrees of freedom do not
3031 int ai = as + ia[i + 1];
3032 int aj = as + ia[i + 2];
3033 if (((atom[ia[i + 1]].ptype == eptNucleus) || (atom[ia[i + 1]].ptype == eptAtom))
3034 && ((atom[ia[i + 2]].ptype == eptNucleus) || (atom[ia[i + 2]].ptype == eptAtom)))
3052 imin = std::min(imin, nrdf2[ai]);
3053 jmin = std::min(jmin, nrdf2[aj]);
3056 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3058 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -=
3060 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3062 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -=
3065 i += interaction_function[ftype].nratoms + 1;
3068 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
3069 for (int i = 0; i < molt.ilist[F_SETTLE].size();)
3071 /* Subtract 1 dof from every atom in the SETTLE */
3072 for (int j = 0; j < 3; j++)
3074 int ai = as + ia[i + 1 + j];
3075 imin = std::min(2, nrdf2[ai]);
3077 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3079 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3084 as += molt.atoms.nr;
3090 /* Correct nrdf for the COM constraints.
3091 * We correct using the TC and VCM group of the first atom
3092 * in the reference and pull group. If atoms in one pull group
3093 * belong to different TC or VCM groups it is anyhow difficult
3094 * to determine the optimal nrdf assignment.
3096 pull = ir->pull.get();
3098 for (int i = 0; i < pull->ncoord; i++)
3100 if (pull->coord[i].eType != epullCONSTRAINT)
3107 for (int j = 0; j < 2; j++)
3109 const t_pull_group* pgrp;
3111 pgrp = &pull->group[pull->coord[i].group[j]];
3113 if (!pgrp->ind.empty())
3115 /* Subtract 1/2 dof from each group */
3116 int ai = pgrp->ind[0];
3117 nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -=
3119 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -=
3121 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
3124 "Center of mass pulling constraints caused the number of degrees "
3125 "of freedom for temperature coupling group %s to be negative",
3126 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(
3127 groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
3132 /* We need to subtract the whole DOF from group j=1 */
3139 if (ir->nstcomm != 0)
3141 GMX_RELEASE_ASSERT(!groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty(),
3142 "Expect at least one group when removing COM motion");
3144 /* We remove COM motion up to dim ndof_com() */
3145 const int ndim_rm_vcm = ndof_com(ir);
3147 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3148 * the number of degrees of freedom in each vcm group when COM
3149 * translation is removed and 6 when rotation is removed as well.
3150 * Note that we do not and should not include the rest group here.
3152 for (gmx::index j = 0;
3153 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]); j++)
3155 switch (ir->comm_mode)
3158 case ecmLINEAR_ACCELERATION_CORRECTION:
3159 nrdf_vcm_sub[j] = 0;
3160 for (int d = 0; d < ndim_rm_vcm; d++)
3168 case ecmANGULAR: nrdf_vcm_sub[j] = 6; break;
3169 default: gmx_incons("Checking comm_mode");
3173 for (gmx::index i = 0;
3174 i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
3176 /* Count the number of atoms of TC group i for every VCM group */
3177 for (gmx::index j = 0;
3178 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3183 for (int ai = 0; ai < natoms; ai++)
3185 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
3187 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
3191 /* Correct for VCM removal according to the fraction of each VCM
3192 * group present in this TC group.
3194 nrdf_uc = nrdf_tc[i];
3196 for (gmx::index j = 0;
3197 j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval]) + 1; j++)
3199 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3201 nrdf_tc[i] += nrdf_uc * (static_cast<double>(na_vcm[j]) / static_cast<double>(na_tot))
3202 * (nrdf_vcm[j] - nrdf_vcm_sub[j]) / nrdf_vcm[j];
3207 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
3209 opts->nrdf[i] = nrdf_tc[i];
3210 if (opts->nrdf[i] < 0)
3214 fprintf(stderr, "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3215 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
3223 sfree(nrdf_vcm_sub);
3226 static bool do_egp_flag(t_inputrec* ir, SimulationGroups* groups, const char* option, const char* val, int flag)
3228 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3229 * But since this is much larger than STRLEN, such a line can not be parsed.
3230 * The real maximum is the number of names that fit in a string: STRLEN/2.
3232 #define EGP_MAX (STRLEN / 2)
3236 auto names = gmx::splitString(val);
3237 if (names.size() % 2 != 0)
3239 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3241 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3243 for (size_t i = 0; i < names.size() / 2; i++)
3245 // TODO this needs to be replaced by a solution using std::find_if
3249 names[2 * i].c_str(),
3250 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
3256 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i].c_str(), option);
3261 names[2 * i + 1].c_str(),
3262 *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
3268 gmx_fatal(FARGS, "%s in %s is not an energy group\n", names[2 * i + 1].c_str(), option);
3270 if ((j < nr) && (k < nr))
3272 ir->opts.egp_flags[nr * j + k] |= flag;
3273 ir->opts.egp_flags[nr * k + j] |= flag;
3282 static void make_swap_groups(t_swapcoords* swap, t_blocka* grps, char** gnames)
3284 int ig = -1, i = 0, gind;
3288 /* Just a quick check here, more thorough checks are in mdrun */
3289 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3291 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3294 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3295 for (ig = 0; ig < swap->ngrp; ig++)
3297 swapg = &swap->grp[ig];
3298 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3299 swapg->nat = grps->index[gind + 1] - grps->index[gind];
3303 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3304 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap", swap->grp[ig].molname, swapg->nat);
3305 snew(swapg->ind, swapg->nat);
3306 for (i = 0; i < swapg->nat; i++)
3308 swapg->ind[i] = grps->a[grps->index[gind] + i];
3313 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3319 static void make_IMD_group(t_IMD* IMDgroup, char* IMDgname, t_blocka* grps, char** gnames)
3324 ig = search_string(IMDgname, grps->nr, gnames);
3325 IMDgroup->nat = grps->index[ig + 1] - grps->index[ig];
3327 if (IMDgroup->nat > 0)
3330 "Group '%s' with %d atoms can be activated for interactive molecular dynamics "
3332 IMDgname, IMDgroup->nat);
3333 snew(IMDgroup->ind, IMDgroup->nat);
3334 for (i = 0; i < IMDgroup->nat; i++)
3336 IMDgroup->ind[i] = grps->a[grps->index[ig] + i];
3341 /* Checks whether atoms are both part of a COM removal group and frozen.
3342 * If a fully frozen atom is part of a COM removal group, it is removed
3343 * from the COM removal group. A note is issued if such atoms are present.
3344 * A warning is issued for atom with one or two dimensions frozen that
3345 * are part of a COM removal group (mdrun would need to compute COM mass
3346 * per dimension to handle this correctly).
3347 * Also issues a warning when non-frozen atoms are not part of a COM
3348 * removal group while COM removal is active.
3350 static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
3352 const t_grpopts& opts,
3355 const int vcmRestGroup =
3356 std::max(int(groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()), 1);
3358 int numFullyFrozenVcmAtoms = 0;
3359 int numPartiallyFrozenVcmAtoms = 0;
3360 int numNonVcmAtoms = 0;
3361 for (int a = 0; a < numAtoms; a++)
3363 const int freezeGroup = getGroupType(*groups, SimulationAtomGroupType::Freeze, a);
3364 int numFrozenDims = 0;
3365 for (int d = 0; d < DIM; d++)
3367 numFrozenDims += opts.nFreeze[freezeGroup][d];
3370 const int vcmGroup = getGroupType(*groups, SimulationAtomGroupType::MassCenterVelocityRemoval, a);
3371 if (vcmGroup < vcmRestGroup)
3373 if (numFrozenDims == DIM)
3375 /* Do not remove COM motion for this fully frozen atom */
3376 if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
3378 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
3381 groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
3382 numFullyFrozenVcmAtoms++;
3384 else if (numFrozenDims > 0)
3386 numPartiallyFrozenVcmAtoms++;
3389 else if (numFrozenDims < DIM)
3395 if (numFullyFrozenVcmAtoms > 0)
3397 std::string warningText = gmx::formatString(
3398 "There are %d atoms that are fully frozen and part of COMM removal group(s), "
3399 "removing these atoms from the COMM removal group(s)",
3400 numFullyFrozenVcmAtoms);
3401 warning_note(wi, warningText.c_str());
3403 if (numPartiallyFrozenVcmAtoms > 0 && numPartiallyFrozenVcmAtoms < numAtoms)
3405 std::string warningText = gmx::formatString(
3406 "There are %d atoms that are frozen along less then %d dimensions and part of COMM "
3407 "removal group(s), due to limitations in the code these still contribute to the "
3408 "mass of the COM along frozen dimensions and therefore the COMM correction will be "
3410 numPartiallyFrozenVcmAtoms, DIM);
3411 warning(wi, warningText.c_str());
3413 if (numNonVcmAtoms > 0)
3415 std::string warningText = gmx::formatString(
3416 "%d atoms are not part of any center of mass motion removal group.\n"
3417 "This may lead to artifacts.\n"
3418 "In most cases one should use one group for the whole system.",
3420 warning(wi, warningText.c_str());
3424 void do_index(const char* mdparin,
3428 const gmx::MdModulesNotifier& notifier,
3432 t_blocka* defaultIndexGroups;
3440 int i, j, k, restnm;
3441 bool bExcl, bTable, bAnneal;
3442 char warn_buf[STRLEN];
3446 fprintf(stderr, "processing index file...\n");
3450 snew(defaultIndexGroups, 1);
3451 snew(defaultIndexGroups->index, 1);
3453 atoms_all = gmx_mtop_global_atoms(mtop);
3454 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3455 done_atom(&atoms_all);
3459 defaultIndexGroups = init_index(ndx, &gnames);
3462 SimulationGroups* groups = &mtop->groups;
3463 natoms = mtop->natoms;
3464 symtab = &mtop->symtab;
3466 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3468 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3470 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3471 restnm = groups->groupNames.size() - 1;
3472 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3473 srenew(gnames, defaultIndexGroups->nr + 1);
3474 gnames[restnm] = *(groups->groupNames.back());
3476 set_warning_line(wi, mdparin, -1);
3478 auto temperatureCouplingTauValues = gmx::splitString(inputrecStrings->tau_t);
3479 auto temperatureCouplingReferenceValues = gmx::splitString(inputrecStrings->ref_t);
3480 auto temperatureCouplingGroupNames = gmx::splitString(inputrecStrings->tcgrps);
3481 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size()
3482 || temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3485 "Invalid T coupling input: %zu groups, %zu ref-t values and "
3487 temperatureCouplingGroupNames.size(), temperatureCouplingReferenceValues.size(),
3488 temperatureCouplingTauValues.size());
3491 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3492 do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
3493 SimulationAtomGroupType::TemperatureCoupling, restnm,
3494 useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3495 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3497 snew(ir->opts.nrdf, nr);
3498 snew(ir->opts.tau_t, nr);
3499 snew(ir->opts.ref_t, nr);
3500 if (ir->eI == eiBD && ir->bd_fric == 0)
3502 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3505 if (useReferenceTemperature)
3507 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3509 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3513 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3514 for (i = 0; (i < nr); i++)
3516 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3518 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3519 warning_error(wi, warn_buf);
3522 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3526 "tau-t = -1 is the value to signal that a group should not have "
3527 "temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3530 if (ir->opts.tau_t[i] >= 0)
3532 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3535 if (ir->etc != etcNO && ir->nsttcouple == -1)
3537 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3542 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3545 "Cannot do Nose-Hoover temperature with Berendsen pressure control with "
3546 "md-vv; use either vrescale temperature with berendsen pressure or "
3547 "Nose-Hoover temperature with MTTK pressure");
3549 if (ir->epc == epcMTTK)
3551 if (ir->etc != etcNOSEHOOVER)
3554 "Cannot do MTTK pressure coupling without Nose-Hoover temperature "
3559 if (ir->nstpcouple != ir->nsttcouple)
3561 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3562 ir->nstpcouple = ir->nsttcouple = mincouple;
3564 "for current Trotter decomposition methods with vv, nsttcouple and "
3565 "nstpcouple must be equal. Both have been reset to "
3566 "min(nsttcouple,nstpcouple) = %d",
3568 warning_note(wi, warn_buf);
3573 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3574 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3576 if (ETC_ANDERSEN(ir->etc))
3578 if (ir->nsttcouple != 1)
3582 "Andersen temperature control methods assume nsttcouple = 1; there is no "
3583 "need for larger nsttcouple > 1, since no global parameters are computed. "
3584 "nsttcouple has been reset to 1");
3585 warning_note(wi, warn_buf);
3588 nstcmin = tcouple_min_integration_steps(ir->etc);
3591 if (tau_min / (ir->delta_t * ir->nsttcouple) < nstcmin - 10 * GMX_REAL_EPS)
3594 "For proper integration of the %s thermostat, tau-t (%g) should be at "
3595 "least %d times larger than nsttcouple*dt (%g)",
3596 ETCOUPLTYPE(ir->etc), tau_min, nstcmin, ir->nsttcouple * ir->delta_t);
3597 warning(wi, warn_buf);
3600 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3601 for (i = 0; (i < nr); i++)
3603 if (ir->opts.ref_t[i] < 0)
3605 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3608 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3609 if we are in this conditional) if mc_temp is negative */
3610 if (ir->expandedvals->mc_temp < 0)
3612 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3616 /* Simulated annealing for each group. There are nr groups */
3617 auto simulatedAnnealingGroupNames = gmx::splitString(inputrecStrings->anneal);
3618 if (simulatedAnnealingGroupNames.size() == 1
3619 && gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3621 simulatedAnnealingGroupNames.resize(0);
3623 if (!simulatedAnnealingGroupNames.empty() && gmx::ssize(simulatedAnnealingGroupNames) != nr)
3625 gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
3626 simulatedAnnealingGroupNames.size(), nr);
3630 snew(ir->opts.annealing, nr);
3631 snew(ir->opts.anneal_npoints, nr);
3632 snew(ir->opts.anneal_time, nr);
3633 snew(ir->opts.anneal_temp, nr);
3634 for (i = 0; i < nr; i++)
3636 ir->opts.annealing[i] = eannNO;
3637 ir->opts.anneal_npoints[i] = 0;
3638 ir->opts.anneal_time[i] = nullptr;
3639 ir->opts.anneal_temp[i] = nullptr;
3641 if (!simulatedAnnealingGroupNames.empty())
3644 for (i = 0; i < nr; i++)
3646 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3648 ir->opts.annealing[i] = eannNO;
3650 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3652 ir->opts.annealing[i] = eannSINGLE;
3655 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3657 ir->opts.annealing[i] = eannPERIODIC;
3663 /* Read the other fields too */
3664 auto simulatedAnnealingPoints = gmx::splitString(inputrecStrings->anneal_npoints);
3665 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3667 gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
3668 simulatedAnnealingPoints.size(), simulatedAnnealingGroupNames.size());
3670 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3671 size_t numSimulatedAnnealingFields = 0;
3672 for (i = 0; i < nr; i++)
3674 if (ir->opts.anneal_npoints[i] == 1)
3678 "Please specify at least a start and an end point for annealing\n");
3680 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3681 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3682 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3685 auto simulatedAnnealingTimes = gmx::splitString(inputrecStrings->anneal_time);
3687 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3689 gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
3690 simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
3692 auto simulatedAnnealingTemperatures = gmx::splitString(inputrecStrings->anneal_temp);
3693 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3695 gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
3696 simulatedAnnealingTemperatures.size(), numSimulatedAnnealingFields);
3699 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3700 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3701 convertReals(wi, simulatedAnnealingTimes, "anneal-time",
3702 allSimulatedAnnealingTimes.data());
3703 convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp",
3704 allSimulatedAnnealingTemperatures.data());
3705 for (i = 0, k = 0; i < nr; i++)
3707 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3709 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3710 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3713 if (ir->opts.anneal_time[i][0] > (ir->init_t + GMX_REAL_EPS))
3715 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3721 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j - 1])
3724 "Annealing timepoints out of order: t=%f comes after "
3726 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j - 1]);
3729 if (ir->opts.anneal_temp[i][j] < 0)
3731 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n",
3732 ir->opts.anneal_temp[i][j]);
3737 /* Print out some summary information, to make sure we got it right */
3738 for (i = 0; i < nr; i++)
3740 if (ir->opts.annealing[i] != eannNO)
3742 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3743 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3744 *(groups->groupNames[j]), eann_names[ir->opts.annealing[i]],
3745 ir->opts.anneal_npoints[i]);
3746 fprintf(stderr, "Time (ps) Temperature (K)\n");
3747 /* All terms except the last one */
3748 for (j = 0; j < (ir->opts.anneal_npoints[i] - 1); j++)
3750 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3751 ir->opts.anneal_temp[i][j]);
3754 /* Finally the last one */
3755 j = ir->opts.anneal_npoints[i] - 1;
3756 if (ir->opts.annealing[i] == eannSINGLE)
3758 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j],
3759 ir->opts.anneal_temp[i][j]);
3763 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j],
3764 ir->opts.anneal_temp[i][j]);
3765 if (std::fabs(ir->opts.anneal_temp[i][j] - ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3768 "There is a temperature jump when your annealing "
3780 for (int i = 1; i < ir->pull->ngroup; i++)
3782 const int gid = search_string(inputrecStrings->pullGroupNames[i].c_str(),
3783 defaultIndexGroups->nr, gnames);
3784 GMX_ASSERT(defaultIndexGroups, "Must have initialized default index groups");
3785 atomGroupRangeValidation(natoms, gid, *defaultIndexGroups);
3788 process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
3790 checkPullCoords(ir->pull->group, ir->pull->coord);
3795 make_rotation_groups(ir->rot, inputrecStrings->rotateGroupNames, defaultIndexGroups, gnames);
3798 if (ir->eSwapCoords != eswapNO)
3800 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3803 /* Make indices for IMD session */
3806 make_IMD_group(ir->imd, inputrecStrings->imd_grp, defaultIndexGroups, gnames);
3809 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3810 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3811 notifier.preProcessingNotifications_.notify(defaultIndexGroupsAndNames);
3813 auto accelerations = gmx::splitString(inputrecStrings->acc);
3814 auto accelerationGroupNames = gmx::splitString(inputrecStrings->accgrps);
3815 if (accelerationGroupNames.size() * DIM != accelerations.size())
3817 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3818 accelerationGroupNames.size(), accelerations.size());
3820 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3821 SimulationAtomGroupType::Acceleration, restnm, egrptpALL_GENREST, bVerbose, wi);
3822 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3823 snew(ir->opts.acc, nr);
3824 ir->opts.ngacc = nr;
3826 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3828 auto freezeDims = gmx::splitString(inputrecStrings->frdim);
3829 auto freezeGroupNames = gmx::splitString(inputrecStrings->freeze);
3830 if (freezeDims.size() != DIM * freezeGroupNames.size())
3832 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3833 freezeGroupNames.size(), freezeDims.size());
3835 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3836 SimulationAtomGroupType::Freeze, restnm, egrptpALL_GENREST, bVerbose, wi);
3837 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3838 ir->opts.ngfrz = nr;
3839 snew(ir->opts.nFreeze, nr);
3840 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3842 for (j = 0; (j < DIM); j++, k++)
3844 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3845 if (!ir->opts.nFreeze[i][j])
3847 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3850 "Please use Y(ES) or N(O) for freezedim only "
3852 freezeDims[k].c_str());
3853 warning(wi, warn_buf);
3858 for (; (i < nr); i++)
3860 for (j = 0; (j < DIM); j++)
3862 ir->opts.nFreeze[i][j] = 0;
3866 auto energyGroupNames = gmx::splitString(inputrecStrings->energy);
3867 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3868 SimulationAtomGroupType::EnergyOutput, restnm, egrptpALL_GENREST, bVerbose, wi);
3869 add_wall_energrps(groups, ir->nwall, symtab);
3870 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3871 auto vcmGroupNames = gmx::splitString(inputrecStrings->vcm);
3872 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3873 SimulationAtomGroupType::MassCenterVelocityRemoval, restnm,
3874 vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3876 if (ir->comm_mode != ecmNO)
3878 checkAndUpdateVcmFreezeGroupConsistency(groups, natoms, ir->opts, wi);
3881 /* Now we have filled the freeze struct, so we can calculate NRDF */
3882 calc_nrdf(mtop, ir, gnames);
3884 auto user1GroupNames = gmx::splitString(inputrecStrings->user1);
3885 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3886 SimulationAtomGroupType::User1, restnm, egrptpALL_GENREST, bVerbose, wi);
3887 auto user2GroupNames = gmx::splitString(inputrecStrings->user2);
3888 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3889 SimulationAtomGroupType::User2, restnm, egrptpALL_GENREST, bVerbose, wi);
3890 auto compressedXGroupNames = gmx::splitString(inputrecStrings->x_compressed_groups);
3891 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3892 SimulationAtomGroupType::CompressedPositionOutput, restnm, egrptpONE, bVerbose, wi);
3893 auto orirefFitGroupNames = gmx::splitString(inputrecStrings->orirefitgrp);
3894 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3895 SimulationAtomGroupType::OrientationRestraintsFit, restnm, egrptpALL_GENREST,
3898 /* MiMiC QMMM input processing */
3899 auto qmGroupNames = gmx::splitString(inputrecStrings->QMMM);
3900 if (qmGroupNames.size() > 1)
3902 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3904 /* group rest, if any, is always MM! */
3905 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3906 SimulationAtomGroupType::QuantumMechanics, restnm, egrptpALL_GENREST, bVerbose, wi);
3907 ir->opts.ngQM = qmGroupNames.size();
3909 /* end of MiMiC QMMM input */
3913 for (auto group : gmx::keysOf(groups->groups))
3915 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3916 for (const auto& entry : groups->groups[group])
3918 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3920 fprintf(stderr, "\n");
3924 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3925 snew(ir->opts.egp_flags, nr * nr);
3927 bExcl = do_egp_flag(ir, groups, "energygrp-excl", inputrecStrings->egpexcl, EGP_EXCL);
3928 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3930 warning_error(wi, "Energy group exclusions are currently not supported");
3932 if (bExcl && EEL_FULL(ir->coulombtype))
3934 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3937 bTable = do_egp_flag(ir, groups, "energygrp-table", inputrecStrings->egptable, EGP_TABLE);
3938 if (bTable && !(ir->vdwtype == evdwUSER) && !(ir->coulombtype == eelUSER)
3939 && !(ir->coulombtype == eelPMEUSER) && !(ir->coulombtype == eelPMEUSERSWITCH))
3942 "Can only have energy group pair tables in combination with user tables for VdW "
3946 /* final check before going out of scope if simulated tempering variables
3947 * need to be set to default values.
3949 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3951 ir->expandedvals->nstexpanded = 2 * static_cast<int>(ir->opts.tau_t[0] / ir->delta_t);
3952 warning(wi, gmx::formatString(
3953 "the value for nstexpanded was not specified for "
3954 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3955 "by default, but it is recommended to set it to an explicit value!",
3956 ir->expandedvals->nstexpanded));
3958 for (i = 0; (i < defaultIndexGroups->nr); i++)
3963 done_blocka(defaultIndexGroups);
3964 sfree(defaultIndexGroups);
3968 static void check_disre(const gmx_mtop_t* mtop)
3970 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3972 const gmx_ffparams_t& ffparams = mtop->ffparams;
3975 for (int i = 0; i < ffparams.numTypes(); i++)
3977 int ftype = ffparams.functype[i];
3978 if (ftype == F_DISRES)
3980 int label = ffparams.iparams[i].disres.label;
3981 if (label == old_label)
3983 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3992 "Found %d double distance restraint indices,\n"
3993 "probably the parameters for multiple pairs in one restraint "
3994 "are not identical\n",
4000 //! Returns whether dimensions have an absolute reference due to walls, pbc or freezing
4001 static BasicVector<bool> haveAbsoluteReference(const t_inputrec& ir)
4003 BasicVector<bool> absRef = { false, false, false };
4005 /* Check the degrees of freedom of the COM (not taking COMM removal into account) */
4006 for (int d = 0; d < DIM; d++)
4008 absRef[d] = (d >= ndof_com(&ir));
4010 /* Check for freeze groups */
4011 for (int g = 0; g < ir.opts.ngfrz; g++)
4013 for (int d = 0; d < DIM; d++)
4015 if (ir.opts.nFreeze[g][d] != 0)
4025 //! Returns whether position restraints are used for dimensions
4026 static BasicVector<bool> havePositionRestraints(const gmx_mtop_t& sys)
4028 BasicVector<bool> havePosres = { false, false, false };
4030 gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(&sys);
4032 while (const InteractionLists* ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
4034 if (nmol > 0 && (!havePosres[XX] || !havePosres[YY] || !havePosres[ZZ]))
4036 for (int i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
4038 const t_iparams& pr = sys.ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
4039 for (int d = 0; d < DIM; d++)
4041 if (pr.posres.fcA[d] != 0)
4043 havePosres[d] = true;
4047 for (int i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
4049 /* Check for flat-bottom posres */
4050 const t_iparams& pr = sys.ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
4051 if (pr.fbposres.k != 0)
4053 switch (pr.fbposres.geom)
4055 case efbposresSPHERE: havePosres = { true, true, true }; break;
4056 case efbposresCYLINDERX: havePosres[YY] = havePosres[ZZ] = true; break;
4057 case efbposresCYLINDERY: havePosres[XX] = havePosres[ZZ] = true; break;
4058 case efbposresCYLINDER:
4059 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
4060 case efbposresCYLINDERZ: havePosres[XX] = havePosres[YY] = true; break;
4061 case efbposresX: /* d=XX */
4062 case efbposresY: /* d=YY */
4063 case efbposresZ: /* d=ZZ */
4064 havePosres[pr.fbposres.geom - efbposresX] = true;
4068 "Invalid geometry for flat-bottom position restraint.\n"
4069 "Expected nr between 1 and %d. Found %d\n",
4070 efbposresNR - 1, pr.fbposres.geom);
4080 static void check_combination_rule_differences(const gmx_mtop_t* mtop,
4082 bool* bC6ParametersWorkWithGeometricRules,
4083 bool* bC6ParametersWorkWithLBRules,
4084 bool* bLBRulesPossible)
4086 int ntypes, tpi, tpj;
4089 double c6i, c6j, c12i, c12j;
4090 double c6, c6_geometric, c6_LB;
4091 double sigmai, sigmaj, epsi, epsj;
4092 bool bCanDoLBRules, bCanDoGeometricRules;
4095 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
4096 * force-field floating point parameters.
4099 ptr = getenv("GMX_LJCOMB_TOL");
4103 double gmx_unused canary;
4105 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
4108 "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4113 *bC6ParametersWorkWithLBRules = TRUE;
4114 *bC6ParametersWorkWithGeometricRules = TRUE;
4115 bCanDoLBRules = TRUE;
4116 ntypes = mtop->ffparams.atnr;
4117 snew(typecount, ntypes);
4118 gmx_mtop_count_atomtypes(mtop, state, typecount);
4119 *bLBRulesPossible = TRUE;
4120 for (tpi = 0; tpi < ntypes; ++tpi)
4122 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4123 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4124 for (tpj = tpi; tpj < ntypes; ++tpj)
4126 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4127 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4128 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4129 c6_geometric = std::sqrt(c6i * c6j);
4130 if (!gmx_numzero(c6_geometric))
4132 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4134 sigmai = gmx::sixthroot(c12i / c6i);
4135 sigmaj = gmx::sixthroot(c12j / c6j);
4136 epsi = c6i * c6i / (4.0 * c12i);
4137 epsj = c6j * c6j / (4.0 * c12j);
4138 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4142 *bLBRulesPossible = FALSE;
4143 c6_LB = c6_geometric;
4145 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4150 *bC6ParametersWorkWithLBRules = FALSE;
4153 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4155 if (!bCanDoGeometricRules)
4157 *bC6ParametersWorkWithGeometricRules = FALSE;
4164 static void check_combination_rules(const t_inputrec* ir, const gmx_mtop_t* mtop, warninp_t wi)
4166 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4168 check_combination_rule_differences(mtop, 0, &bC6ParametersWorkWithGeometricRules,
4169 &bC6ParametersWorkWithLBRules, &bLBRulesPossible);
4170 if (ir->ljpme_combination_rule == eljpmeLB)
4172 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
4175 "You are using arithmetic-geometric combination rules "
4176 "in LJ-PME, but your non-bonded C6 parameters do not "
4177 "follow these rules.");
4182 if (!bC6ParametersWorkWithGeometricRules)
4184 if (ir->eDispCorr != edispcNO)
4187 "You are using geometric combination rules in "
4188 "LJ-PME, but your non-bonded C6 parameters do "
4189 "not follow these rules. "
4190 "This will introduce very small errors in the forces and energies in "
4191 "your simulations. Dispersion correction will correct total energy "
4192 "and/or pressure for isotropic systems, but not forces or surface "
4198 "You are using geometric combination rules in "
4199 "LJ-PME, but your non-bonded C6 parameters do "
4200 "not follow these rules. "
4201 "This will introduce very small errors in the forces and energies in "
4202 "your simulations. If your system is homogeneous, consider using "
4203 "dispersion correction "
4204 "for the total energy and pressure.");
4210 static bool allTrue(const BasicVector<bool>& boolVector)
4212 return boolVector[0] && boolVector[1] && boolVector[2];
4215 void triple_check(const char* mdparin, t_inputrec* ir, gmx_mtop_t* sys, warninp_t wi)
4217 // Not meeting MTS requirements should have resulted in a fatal error, so we can assert here
4218 gmx::assertMtsRequirements(*ir);
4220 char err_buf[STRLEN];
4225 gmx_mtop_atomloop_block_t aloopb;
4226 char warn_buf[STRLEN];
4228 set_warning_line(wi, mdparin, -1);
4230 if (allTrue(haveAbsoluteReference(*ir)) && allTrue(havePositionRestraints(*sys)))
4233 "Removing center of mass motion in the presence of position restraints might "
4234 "cause artifacts. When you are using position restraints to equilibrate a "
4235 "macro-molecule, the artifacts are usually negligible.");
4238 if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0 && ir->nstlist > 1
4239 && ((EI_MD(ir->eI) || EI_SD(ir->eI)) && (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4241 /* Check if a too small Verlet buffer might potentially
4242 * cause more drift than the thermostat can couple off.
4244 /* Temperature error fraction for warning and suggestion */
4245 const real T_error_warn = 0.002;
4246 const real T_error_suggest = 0.001;
4247 /* For safety: 2 DOF per atom (typical with constraints) */
4248 const real nrdf_at = 2;
4249 real T, tau, max_T_error;
4254 for (i = 0; i < ir->opts.ngtc; i++)
4256 T = std::max(T, ir->opts.ref_t[i]);
4257 tau = std::max(tau, ir->opts.tau_t[i]);
4261 /* This is a worst case estimate of the temperature error,
4262 * assuming perfect buffer estimation and no cancelation
4263 * of errors. The factor 0.5 is because energy distributes
4264 * equally over Ekin and Epot.
4266 max_T_error = 0.5 * tau * ir->verletbuf_tol / (nrdf_at * BOLTZ * T);
4267 if (max_T_error > T_error_warn)
4270 "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature "
4271 "of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. "
4272 "To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to "
4273 "%.0e or decrease tau_t.",
4274 ir->verletbuf_tol, T, tau, 100 * max_T_error, 100 * T_error_suggest,
4275 ir->verletbuf_tol * T_error_suggest / max_T_error);
4276 warning(wi, warn_buf);
4281 if (ETC_ANDERSEN(ir->etc))
4285 for (i = 0; i < ir->opts.ngtc; i++)
4288 "all tau_t must currently be equal using Andersen temperature control, "
4289 "violated for group %d",
4291 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4293 "all tau_t must be positive using Andersen temperature control, "
4295 i, ir->opts.tau_t[i]);
4296 CHECK(ir->opts.tau_t[i] < 0);
4299 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
4301 for (i = 0; i < ir->opts.ngtc; i++)
4303 int nsteps = gmx::roundToInt(ir->opts.tau_t[i] / ir->delta_t);
4305 "tau_t/delta_t for group %d for temperature control method %s must be a "
4306 "multiple of nstcomm (%d), as velocities of atoms in coupled groups are "
4307 "randomized every time step. The input tau_t (%8.3f) leads to %d steps per "
4309 i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
4310 CHECK(nsteps % ir->nstcomm != 0);
4315 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD && ir->comm_mode == ecmNO
4316 && !(allTrue(haveAbsoluteReference(*ir)) || allTrue(havePositionRestraints(*sys)) || ir->nsteps <= 10)
4317 && !ETC_ANDERSEN(ir->etc))
4320 "You are not using center of mass motion removal (mdp option comm-mode), numerical "
4321 "rounding errors can lead to build up of kinetic energy of the center of mass");
4324 if (ir->epc == epcPARRINELLORAHMAN && ir->etc == etcNOSEHOOVER)
4327 for (int g = 0; g < ir->opts.ngtc; g++)
4329 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
4331 if (ir->tau_p < 1.9 * tau_t_max)
4333 std::string message = gmx::formatString(
4334 "With %s T-coupling and %s p-coupling, "
4335 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
4336 etcoupl_names[ir->etc], epcoupl_names[ir->epc], "tau-p", ir->tau_p, "tau-t",
4338 warning(wi, message.c_str());
4342 /* Check for pressure coupling with absolute position restraints */
4343 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4345 const BasicVector<bool> havePosres = havePositionRestraints(*sys);
4347 for (m = 0; m < DIM; m++)
4349 if (havePosres[m] && norm2(ir->compress[m]) > 0)
4352 "You are using pressure coupling with absolute position restraints, "
4353 "this will give artifacts. Use the refcoord_scaling option.");
4361 aloopb = gmx_mtop_atomloop_block_init(sys);
4363 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4365 if (atom->q != 0 || atom->qB != 0)
4373 if (EEL_FULL(ir->coulombtype))
4376 "You are using full electrostatics treatment %s for a system without charges.\n"
4377 "This costs a lot of performance for just processing zeros, consider using %s "
4379 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4380 warning(wi, err_buf);
4385 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4388 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4389 "You might want to consider using %s electrostatics.\n",
4391 warning_note(wi, err_buf);
4395 /* Check if combination rules used in LJ-PME are the same as in the force field */
4396 if (EVDW_PME(ir->vdwtype))
4398 check_combination_rules(ir, sys, wi);
4401 /* Generalized reaction field */
4402 if (ir->coulombtype == eelGRF_NOTUSED)
4405 "Generalized reaction-field electrostatics is no longer supported. "
4406 "You can use normal reaction-field instead and compute the reaction-field "
4407 "constant by hand.");
4411 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4413 for (m = 0; (m < DIM); m++)
4415 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4424 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4425 for (const AtomProxy atomP : AtomRange(*sys))
4427 const t_atom& local = atomP.atom();
4428 int i = atomP.globalAtomNumber();
4429 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4432 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4434 for (m = 0; (m < DIM); m++)
4436 acc[m] += ir->opts.acc[i][m] * mgrp[i];
4440 for (m = 0; (m < DIM); m++)
4442 if (fabs(acc[m]) > 1e-6)
4444 const char* dim[DIM] = { "X", "Y", "Z" };
4445 fprintf(stderr, "Net Acceleration in %s direction, will %s be corrected\n", dim[m],
4446 ir->nstcomm != 0 ? "" : "not");
4447 if (ir->nstcomm != 0 && m < ndof_com(ir))
4451 (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4453 ir->opts.acc[i][m] -= acc[m];
4461 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0
4462 && !gmx_within_tol(sys->ffparams.reppow, 12.0, 10 * GMX_DOUBLE_EPS))
4464 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4472 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4474 if (ir->pull->coord[i].group[0] == 0 || ir->pull->coord[i].group[1] == 0)
4476 const auto absRef = haveAbsoluteReference(*ir);
4477 const auto havePosres = havePositionRestraints(*sys);
4478 for (m = 0; m < DIM; m++)
4480 if (ir->pull->coord[i].dim[m] && !(absRef[m] || havePosres[m]))
4483 "You are using an absolute reference for pulling, but the rest of "
4484 "the system does not have an absolute reference. This will lead to "
4493 for (i = 0; i < 3; i++)
4495 for (m = 0; m <= i; m++)
4497 if ((ir->epc != epcNO && ir->compress[i][m] != 0) || ir->deform[i][m] != 0)
4499 for (c = 0; c < ir->pull->ncoord; c++)
4501 if (ir->pull->coord[c].eGeom == epullgDIRPBC && ir->pull->coord[c].vec[m] != 0)
4504 "Can not have dynamic box while using pull geometry '%s' "
4506 EPULLGEOM(ir->pull->coord[c].eGeom), 'x' + m);
4517 void double_check(t_inputrec* ir, matrix box, bool bHasNormalConstraints, bool bHasAnyConstraints, warninp_t wi)
4519 char warn_buf[STRLEN];
4522 ptr = check_box(ir->pbcType, box);
4525 warning_error(wi, ptr);
4528 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4530 if (ir->shake_tol <= 0.0)
4532 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n", ir->shake_tol);
4533 warning_error(wi, warn_buf);
4537 if ((ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4539 /* If we have Lincs constraints: */
4540 if (ir->eI == eiMD && ir->etc == etcNO && ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4543 "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4544 warning_note(wi, warn_buf);
4547 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4550 "For accurate %s with LINCS constraints, lincs-order should be 8 or more.",
4552 warning_note(wi, warn_buf);
4554 if (ir->epc == epcMTTK)
4556 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4560 if (bHasAnyConstraints && ir->epc == epcMTTK)
4562 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4565 if (ir->LincsWarnAngle > 90.0)
4567 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4568 warning(wi, warn_buf);
4569 ir->LincsWarnAngle = 90.0;
4572 if (ir->pbcType != PbcType::No)
4574 if (ir->nstlist == 0)
4577 "With nstlist=0 atoms are only put into the box at step 0, therefore drifting "
4578 "atoms might cause the simulation to crash.");
4580 if (gmx::square(ir->rlist) >= max_cutoff2(ir->pbcType, box))
4583 "ERROR: The cut-off length is longer than half the shortest box vector or "
4584 "longer than the smallest box diagonal element. Increase the box size or "
4585 "decrease rlist.\n");
4586 warning_error(wi, warn_buf);