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,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/awh/read_params.h"
50 #include "gromacs/fileio/readinp.h"
51 #include "gromacs/fileio/warninp.h"
52 #include "gromacs/gmxlib/chargegroup.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/gmxpreprocess/toputil.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/calc_verletbuf.h"
59 #include "gromacs/mdrun/mdmodules.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/mdtypes/pull_params.h"
63 #include "gromacs/options/options.h"
64 #include "gromacs/options/treesupport.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/selection/indexutil.h"
67 #include "gromacs/topology/block.h"
68 #include "gromacs/topology/ifunc.h"
69 #include "gromacs/topology/index.h"
70 #include "gromacs/topology/mtop_util.h"
71 #include "gromacs/topology/symtab.h"
72 #include "gromacs/topology/topology.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/exceptions.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/filestream.h"
77 #include "gromacs/utility/gmxassert.h"
78 #include "gromacs/utility/ikeyvaluetreeerror.h"
79 #include "gromacs/utility/keyvaluetree.h"
80 #include "gromacs/utility/keyvaluetreebuilder.h"
81 #include "gromacs/utility/keyvaluetreemdpwriter.h"
82 #include "gromacs/utility/keyvaluetreetransform.h"
83 #include "gromacs/utility/smalloc.h"
84 #include "gromacs/utility/strconvert.h"
85 #include "gromacs/utility/stringcompare.h"
86 #include "gromacs/utility/stringutil.h"
87 #include "gromacs/utility/textwriter.h"
92 /* Resource parameters
93 * Do not change any of these until you read the instruction
94 * in readinp.h. Some cpp's do not take spaces after the backslash
95 * (like the c-shell), which will give you a very weird compiler
99 typedef struct t_inputrec_strings
101 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
102 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
103 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
104 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
105 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
107 char fep_lambda[efptNR][STRLEN];
108 char lambda_weights[STRLEN];
111 char anneal[STRLEN], anneal_npoints[STRLEN],
112 anneal_time[STRLEN], anneal_temp[STRLEN];
113 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
114 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
115 SAoff[STRLEN], SAsteps[STRLEN];
117 } gmx_inputrec_strings;
119 static gmx_inputrec_strings *is = nullptr;
121 void init_inputrec_strings()
125 gmx_incons("Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time.");
130 void done_inputrec_strings()
138 egrptpALL, /* All particles have to be a member of a group. */
139 egrptpALL_GENREST, /* A rest group with name is generated for particles *
140 * that are not part of any group. */
141 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
142 * for the rest group. */
143 egrptpONE /* Merge all selected groups into one group, *
144 * make a rest group for the remaining particles. */
147 static const char *constraints[eshNR+1] = {
148 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
151 static const char *couple_lam[ecouplamNR+1] = {
152 "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] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
167 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
169 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
171 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
173 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
178 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
179 gmx_fatal(FARGS, "%s", errorstr);
186 static void _low_check(bool b, const char *s, warninp_t wi)
190 warning_error(wi, s);
194 static void check_nst(const char *desc_nst, int nst,
195 const char *desc_p, int *p,
200 if (*p > 0 && *p % nst != 0)
202 /* Round up to the next multiple of nst */
203 *p = ((*p)/nst + 1)*nst;
204 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
205 desc_p, desc_nst, desc_p, *p);
210 static bool ir_NVE(const t_inputrec *ir)
212 return (EI_MD(ir->eI) && ir->etc == etcNO);
215 static int lcd(int n1, int n2)
220 for (i = 2; (i <= n1 && i <= n2); i++)
222 if (n1 % i == 0 && n2 % i == 0)
231 //! Convert legacy mdp entries to modern ones.
232 static void process_interaction_modifier(int *eintmod)
234 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
236 *eintmod = eintmodPOTSHIFT;
240 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
242 /* Check internal consistency.
243 * NOTE: index groups are not set here yet, don't check things
244 * like temperature coupling group options here, but in triple_check
247 /* Strange macro: first one fills the err_buf, and then one can check
248 * the condition, which will print the message and increase the error
251 #define CHECK(b) _low_check(b, err_buf, wi)
252 char err_buf[256], warn_buf[STRLEN];
255 t_lambda *fep = ir->fepvals;
256 t_expanded *expand = ir->expandedvals;
258 set_warning_line(wi, mdparin, -1);
260 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
262 sprintf(warn_buf, "%s electrostatics is no longer supported",
263 eel_names[eelRF_NEC_UNSUPPORTED]);
264 warning_error(wi, warn_buf);
267 /* BASIC CUT-OFF STUFF */
268 if (ir->rcoulomb < 0)
270 warning_error(wi, "rcoulomb should be >= 0");
274 warning_error(wi, "rvdw should be >= 0");
277 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
279 warning_error(wi, "rlist should be >= 0");
281 sprintf(err_buf, "nstlist can not be smaller than 0. (If you were trying to use the heuristic neighbour-list update scheme for efficient buffering for improved energy conservation, please use the Verlet cut-off scheme instead.)");
282 CHECK(ir->nstlist < 0);
284 process_interaction_modifier(&ir->coulomb_modifier);
285 process_interaction_modifier(&ir->vdw_modifier);
287 if (ir->cutoff_scheme == ecutsGROUP)
290 "The group cutoff scheme has been removed since GROMACS 2020. "
291 "Please use the Verlet cutoff scheme.");
293 if (ir->cutoff_scheme == ecutsVERLET)
297 /* Normal Verlet type neighbor-list, currently only limited feature support */
298 if (inputrec2nboundeddim(ir) < 3)
300 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
303 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
304 if (ir->rcoulomb != ir->rvdw)
306 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
307 // for PME load balancing, we can support this exception.
308 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
309 ir->vdwtype == evdwCUT &&
310 ir->rcoulomb > ir->rvdw);
311 if (!bUsesPmeTwinRangeKernel)
313 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
317 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
319 if (ir->vdw_modifier == eintmodNONE ||
320 ir->vdw_modifier == eintmodPOTSHIFT)
322 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
324 sprintf(warn_buf, "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
325 warning_note(wi, warn_buf);
327 ir->vdwtype = evdwCUT;
331 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
332 warning_error(wi, warn_buf);
336 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
338 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
340 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
341 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
343 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
345 if (!(ir->coulomb_modifier == eintmodNONE ||
346 ir->coulomb_modifier == eintmodPOTSHIFT))
348 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
349 warning_error(wi, warn_buf);
352 if (EEL_USER(ir->coulombtype))
354 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme", eel_names[ir->coulombtype]);
355 warning_error(wi, warn_buf);
358 if (ir->nstlist <= 0)
360 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
363 if (ir->nstlist < 10)
365 warning_note(wi, "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
368 rc_max = std::max(ir->rvdw, ir->rcoulomb);
370 if (ir->verletbuf_tol <= 0)
372 if (ir->verletbuf_tol == 0)
374 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
377 if (ir->rlist < rc_max)
379 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
382 if (ir->rlist == rc_max && ir->nstlist > 1)
384 warning_note(wi, "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
389 if (ir->rlist > rc_max)
391 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
394 if (ir->nstlist == 1)
396 /* No buffer required */
401 if (EI_DYNAMICS(ir->eI))
403 if (inputrec2nboundeddim(ir) < 3)
405 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
407 /* Set rlist temporarily so we can continue processing */
412 /* Set the buffer to 5% of the cut-off */
413 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
419 /* GENERAL INTEGRATOR STUFF */
422 if (ir->etc != etcNO)
424 if (EI_RANDOM(ir->eI))
426 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling implicitly. See the documentation for more information on which parameters affect temperature for %s.", etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
430 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
432 warning_note(wi, warn_buf);
436 if (ir->eI == eiVVAK)
438 sprintf(warn_buf, "Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s", ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
439 warning_note(wi, warn_buf);
441 if (!EI_DYNAMICS(ir->eI))
443 if (ir->epc != epcNO)
445 sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
446 warning_note(wi, warn_buf);
450 if (EI_DYNAMICS(ir->eI))
452 if (ir->nstcalcenergy < 0)
454 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
455 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
457 /* nstcalcenergy larger than nstener does not make sense.
458 * We ideally want nstcalcenergy=nstener.
462 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
466 ir->nstcalcenergy = ir->nstenergy;
470 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
471 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
472 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
475 const char *nsten = "nstenergy";
476 const char *nstdh = "nstdhdl";
477 const char *min_name = nsten;
478 int min_nst = ir->nstenergy;
480 /* find the smallest of ( nstenergy, nstdhdl ) */
481 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
482 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
484 min_nst = ir->fepvals->nstdhdl;
487 /* If the user sets nstenergy small, we should respect that */
489 "Setting nstcalcenergy (%d) equal to %s (%d)",
490 ir->nstcalcenergy, min_name, min_nst);
491 warning_note(wi, warn_buf);
492 ir->nstcalcenergy = min_nst;
495 if (ir->epc != epcNO)
497 if (ir->nstpcouple < 0)
499 ir->nstpcouple = ir_optimal_nstpcouple(ir);
503 if (ir->nstcalcenergy > 0)
505 if (ir->efep != efepNO)
507 /* nstdhdl should be a multiple of nstcalcenergy */
508 check_nst("nstcalcenergy", ir->nstcalcenergy,
509 "nstdhdl", &ir->fepvals->nstdhdl, wi);
513 /* nstexpanded should be a multiple of nstcalcenergy */
514 check_nst("nstcalcenergy", ir->nstcalcenergy,
515 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
517 /* for storing exact averages nstenergy should be
518 * a multiple of nstcalcenergy
520 check_nst("nstcalcenergy", ir->nstcalcenergy,
521 "nstenergy", &ir->nstenergy, wi);
525 if (ir->nsteps == 0 && !ir->bContinuation)
527 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
531 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
532 ir->bContinuation && ir->ld_seed != -1)
534 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
540 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
541 CHECK(ir->ePBC != epbcXYZ);
542 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
543 CHECK(ir->ns_type != ensGRID);
544 sprintf(err_buf, "with TPI nstlist should be larger than zero");
545 CHECK(ir->nstlist <= 0);
546 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
547 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
548 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
549 CHECK(ir->cutoff_scheme == ecutsVERLET);
553 if ( (opts->nshake > 0) && (opts->bMorse) )
556 "Using morse bond-potentials while constraining bonds is useless");
557 warning(wi, warn_buf);
560 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
561 ir->bContinuation && ir->ld_seed != -1)
563 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
565 /* verify simulated tempering options */
569 bool bAllTempZero = TRUE;
570 for (i = 0; i < fep->n_lambda; i++)
572 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
573 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
574 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
576 bAllTempZero = FALSE;
579 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
580 CHECK(bAllTempZero == TRUE);
582 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
583 CHECK(ir->eI != eiVV);
585 /* check compatability of the temperature coupling with simulated tempering */
587 if (ir->etc == etcNOSEHOOVER)
589 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
590 warning_note(wi, warn_buf);
593 /* check that the temperatures make sense */
595 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= than the simulated tempering lower temperature (%g)", ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
596 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
598 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
599 CHECK(ir->simtempvals->simtemp_high <= 0);
601 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
602 CHECK(ir->simtempvals->simtemp_low <= 0);
605 /* verify free energy options */
607 if (ir->efep != efepNO)
610 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
612 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
614 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
615 static_cast<int>(fep->sc_r_power));
616 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
618 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
619 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
621 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
622 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
624 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
625 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
627 sprintf(err_buf, "Free-energy not implemented for Ewald");
628 CHECK(ir->coulombtype == eelEWALD);
630 /* check validty of lambda inputs */
631 if (fep->n_lambda == 0)
633 /* Clear output in case of no states:*/
634 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
635 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
639 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
640 CHECK((fep->init_fep_state >= fep->n_lambda));
643 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
644 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
646 sprintf(err_buf, "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with init-lambda-state or with init-lambda, but not both",
647 fep->init_lambda, fep->init_fep_state);
648 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
652 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
656 for (i = 0; i < efptNR; i++)
658 if (fep->separate_dvdl[i])
663 if (n_lambda_terms > 1)
665 sprintf(warn_buf, "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't use init-lambda to set lambda state (except for slow growth). Use init-lambda-state instead.");
666 warning(wi, warn_buf);
669 if (n_lambda_terms < 2 && fep->n_lambda > 0)
672 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
676 for (j = 0; j < efptNR; j++)
678 for (i = 0; i < fep->n_lambda; i++)
680 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[j], fep->all_lambda[j][i]);
681 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
685 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
687 for (i = 0; i < fep->n_lambda; i++)
689 sprintf(err_buf, "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to crashes, and is not supported.", i, fep->all_lambda[efptVDW][i],
690 fep->all_lambda[efptCOUL][i]);
691 CHECK((fep->sc_alpha > 0) &&
692 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
693 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
694 ((fep->all_lambda[efptVDW][i] > 0.0) &&
695 (fep->all_lambda[efptVDW][i] < 1.0))));
699 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
701 real sigma, lambda, r_sc;
704 /* Maximum estimate for A and B charges equal with lambda power 1 */
706 r_sc = std::pow(lambda*fep->sc_alpha*std::pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
707 sprintf(warn_buf, "With PME there is a minor soft core effect present at the cut-off, proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on energy conservation, but usually other effects dominate. With a common sigma value of %g nm the fraction of the particle-particle potential at the cut-off at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
709 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
710 warning_note(wi, warn_buf);
713 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
714 be treated differently, but that's the next step */
716 for (i = 0; i < efptNR; i++)
718 for (j = 0; j < fep->n_lambda; j++)
720 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
721 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
726 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
730 /* checking equilibration of weights inputs for validity */
732 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
733 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
734 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
736 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
737 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
738 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
740 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
741 expand->equil_steps, elmceq_names[elmceqSTEPS]);
742 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
744 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
745 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
746 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
748 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
749 expand->equil_ratio, elmceq_names[elmceqRATIO]);
750 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
752 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
753 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
754 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
756 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
757 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
758 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
760 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
761 expand->equil_steps, elmceq_names[elmceqSTEPS]);
762 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
764 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
765 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
766 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
768 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
769 expand->equil_ratio, elmceq_names[elmceqRATIO]);
770 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
772 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
773 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
774 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
776 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
777 CHECK((expand->lmc_repeats <= 0));
778 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
779 CHECK((expand->minvarmin <= 0));
780 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
781 CHECK((expand->c_range < 0));
782 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
783 fep->init_fep_state, expand->lmc_forced_nstart);
784 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
785 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
786 CHECK((expand->lmc_forced_nstart < 0));
787 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
788 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
790 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
791 CHECK((expand->init_wl_delta < 0));
792 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
793 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
794 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
795 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
797 /* if there is no temperature control, we need to specify an MC temperature */
798 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO) && (expand->mc_temp <= 0.0))
800 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set to a positive number");
801 warning_error(wi, err_buf);
803 if (expand->nstTij > 0)
805 sprintf(err_buf, "nstlog must be non-zero");
806 CHECK(ir->nstlog == 0);
807 // Avoid modulus by zero in the case that already triggered an error exit.
810 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
811 expand->nstTij, ir->nstlog);
812 CHECK((expand->nstTij % ir->nstlog) != 0);
818 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
819 CHECK(ir->nwall && ir->ePBC != epbcXY);
822 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
824 if (ir->ePBC == epbcNONE)
826 if (ir->epc != epcNO)
828 warning(wi, "Turning off pressure coupling for vacuum system");
834 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
835 epbc_names[ir->ePBC]);
836 CHECK(ir->epc != epcNO);
838 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
839 CHECK(EEL_FULL(ir->coulombtype));
841 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
842 epbc_names[ir->ePBC]);
843 CHECK(ir->eDispCorr != edispcNO);
846 if (ir->rlist == 0.0)
848 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
849 "with coulombtype = %s or coulombtype = %s\n"
850 "without periodic boundary conditions (pbc = %s) and\n"
851 "rcoulomb and rvdw set to zero",
852 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
853 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
854 (ir->ePBC != epbcNONE) ||
855 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
859 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
864 if (ir->nstcomm == 0)
866 // TODO Change this behaviour. There should be exactly one way
867 // to turn off an algorithm.
868 ir->comm_mode = ecmNO;
870 if (ir->comm_mode != ecmNO)
874 // TODO Such input was once valid. Now that we've been
875 // helpful for a few years, we should reject such input,
876 // lest we have to support every historical decision
878 warning(wi, "If you want to remove the rotation around the center of mass, you should set comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to its absolute value");
879 ir->nstcomm = abs(ir->nstcomm);
882 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
884 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
885 ir->nstcomm = ir->nstcalcenergy;
888 if (ir->comm_mode == ecmANGULAR)
890 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
891 CHECK(ir->bPeriodicMols);
892 if (ir->ePBC != epbcNONE)
894 warning(wi, "Removing the rotation around the center of mass in a periodic system, this can lead to artifacts. Only use this on a single (cluster of) molecules. This cluster should not cross periodic boundaries.");
899 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
901 sprintf(warn_buf, "Tumbling and flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR or use integrator = %s.", ei_names[eiSD1]);
902 warning_note(wi, warn_buf);
905 /* TEMPERATURE COUPLING */
906 if (ir->etc == etcYES)
908 ir->etc = etcBERENDSEN;
909 warning_note(wi, "Old option for temperature coupling given: "
910 "changing \"yes\" to \"Berendsen\"\n");
913 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
915 if (ir->opts.nhchainlength < 1)
917 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
918 ir->opts.nhchainlength = 1;
919 warning(wi, warn_buf);
922 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
924 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
925 ir->opts.nhchainlength = 1;
930 ir->opts.nhchainlength = 0;
933 if (ir->eI == eiVVAK)
935 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
937 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
940 if (ETC_ANDERSEN(ir->etc))
942 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
943 CHECK(!(EI_VV(ir->eI)));
945 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
947 sprintf(warn_buf, "Center of mass removal not necessary for %s. All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.", etcoupl_names[ir->etc]);
948 warning_note(wi, warn_buf);
951 sprintf(err_buf, "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step", ir->nstcomm, etcoupl_names[ir->etc]);
952 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
955 if (ir->etc == etcBERENDSEN)
957 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
958 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
959 warning_note(wi, warn_buf);
962 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
963 && ir->epc == epcBERENDSEN)
965 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
966 "true ensemble for the thermostat");
967 warning(wi, warn_buf);
970 /* PRESSURE COUPLING */
971 if (ir->epc == epcISOTROPIC)
973 ir->epc = epcBERENDSEN;
974 warning_note(wi, "Old option for pressure coupling given: "
975 "changing \"Isotropic\" to \"Berendsen\"\n");
978 if (ir->epc != epcNO)
980 dt_pcoupl = ir->nstpcouple*ir->delta_t;
982 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
983 CHECK(ir->tau_p <= 0);
985 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
987 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
988 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
989 warning(wi, warn_buf);
992 sprintf(err_buf, "compressibility must be > 0 when using pressure"
993 " coupling %s\n", EPCOUPLTYPE(ir->epc));
994 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
995 ir->compress[ZZ][ZZ] < 0 ||
996 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
997 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
999 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1002 "You are generating velocities so I am assuming you "
1003 "are equilibrating a system. You are using "
1004 "%s pressure coupling, but this can be "
1005 "unstable for equilibration. If your system crashes, try "
1006 "equilibrating first with Berendsen pressure coupling. If "
1007 "you are not equilibrating the system, you can probably "
1008 "ignore this warning.",
1009 epcoupl_names[ir->epc]);
1010 warning(wi, warn_buf);
1016 if (ir->epc > epcNO)
1018 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1020 warning_error(wi, "for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
1026 if (ir->epc == epcMTTK)
1028 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1032 /* ELECTROSTATICS */
1033 /* More checks are in triple check (grompp.c) */
1035 if (ir->coulombtype == eelSWITCH)
1037 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1038 "artifacts, advice: use coulombtype = %s",
1039 eel_names[ir->coulombtype],
1040 eel_names[eelRF_ZERO]);
1041 warning(wi, warn_buf);
1044 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1046 sprintf(warn_buf, "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old format and exchanging epsilon-r and epsilon-rf", ir->epsilon_r);
1047 warning(wi, warn_buf);
1048 ir->epsilon_rf = ir->epsilon_r;
1049 ir->epsilon_r = 1.0;
1052 if (ir->epsilon_r == 0)
1055 "It is pointless to use long-range electrostatics with infinite relative permittivity."
1056 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1057 CHECK(EEL_FULL(ir->coulombtype));
1060 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1062 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1063 CHECK(ir->epsilon_r < 0);
1066 if (EEL_RF(ir->coulombtype))
1068 /* reaction field (at the cut-off) */
1070 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1072 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1073 eel_names[ir->coulombtype]);
1074 warning(wi, warn_buf);
1075 ir->epsilon_rf = 0.0;
1078 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1079 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1080 (ir->epsilon_r == 0));
1081 if (ir->epsilon_rf == ir->epsilon_r)
1083 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1084 eel_names[ir->coulombtype]);
1085 warning(wi, warn_buf);
1088 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1089 * means the interaction is zero outside rcoulomb, but it helps to
1090 * provide accurate energy conservation.
1092 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1094 if (ir_coulomb_switched(ir))
1097 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1098 eel_names[ir->coulombtype]);
1099 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1103 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1106 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1107 CHECK( ir->coulomb_modifier != eintmodNONE);
1109 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1112 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1113 CHECK( ir->vdw_modifier != eintmodNONE);
1116 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1117 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1120 "The switch/shift interaction settings are just for compatibility; you will get better "
1121 "performance from applying potential modifiers to your interactions!\n");
1122 warning_note(wi, warn_buf);
1125 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1127 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1129 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1130 sprintf(warn_buf, "The switching range should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
1131 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1132 warning(wi, warn_buf);
1136 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1138 if (ir->rvdw_switch == 0)
1140 sprintf(warn_buf, "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones potential. This suggests it was not set in the mdp, which can lead to large energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw switching range.");
1141 warning(wi, warn_buf);
1145 if (EEL_FULL(ir->coulombtype))
1147 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1148 ir->coulombtype == eelPMEUSERSWITCH)
1150 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1151 eel_names[ir->coulombtype]);
1152 CHECK(ir->rcoulomb > ir->rlist);
1156 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1158 // TODO: Move these checks into the ewald module with the options class
1160 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1162 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1164 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
1165 warning_error(wi, warn_buf);
1169 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1171 if (ir->ewald_geometry == eewg3D)
1173 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1174 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1175 warning(wi, warn_buf);
1177 /* This check avoids extra pbc coding for exclusion corrections */
1178 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1179 CHECK(ir->wall_ewald_zfac < 2);
1181 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1182 EEL_FULL(ir->coulombtype))
1184 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1185 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1186 warning(wi, warn_buf);
1188 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1190 if (ir->cutoff_scheme == ecutsVERLET)
1192 sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
1193 eel_names[ir->coulombtype]);
1194 warning(wi, warn_buf);
1198 sprintf(warn_buf, "Dipole corrections to %s electrostatics only work if all charge groups that can cross PBC boundaries are dipoles. If this is not the case set epsilon_surface to 0",
1199 eel_names[ir->coulombtype]);
1200 warning_note(wi, warn_buf);
1204 if (ir_vdw_switched(ir))
1206 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1207 CHECK(ir->rvdw_switch >= ir->rvdw);
1209 if (ir->rvdw_switch < 0.5*ir->rvdw)
1211 sprintf(warn_buf, "You are applying a switch function to vdw forces or potentials from %g to %g nm, which is more than half the interaction range, whereas switch functions are intended to act only close to the cut-off.",
1212 ir->rvdw_switch, ir->rvdw);
1213 warning_note(wi, warn_buf);
1217 if (ir->vdwtype == evdwPME)
1219 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1221 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1222 evdw_names[ir->vdwtype],
1223 eintmod_names[eintmodPOTSHIFT],
1224 eintmod_names[eintmodNONE]);
1225 warning_error(wi, err_buf);
1229 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1231 warning_note(wi, "You have selected user tables with dispersion correction, the dispersion will be corrected to -C6/r^6 beyond rvdw_switch (the tabulated interaction between rvdw_switch and rvdw will not be double counted). Make sure that you really want dispersion correction to -C6/r^6.");
1234 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1237 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1240 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1242 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1245 /* IMPLICIT SOLVENT */
1246 if (ir->coulombtype == eelGB_NOTUSED)
1248 sprintf(warn_buf, "Invalid option %s for coulombtype",
1249 eel_names[ir->coulombtype]);
1250 warning_error(wi, warn_buf);
1255 if (ir->cutoff_scheme != ecutsGROUP)
1257 warning_error(wi, "QMMM is currently only supported with cutoff-scheme=group");
1259 if (!EI_DYNAMICS(ir->eI))
1262 sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
1263 warning_error(wi, buf);
1269 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1273 /* interpret a number of doubles from a string and put them in an array,
1274 after allocating space for them.
1275 str = the input string
1276 n = the (pre-allocated) number of doubles read
1277 r = the output array of doubles. */
1278 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1280 auto values = gmx::splitString(str);
1284 for (int i = 0; i < *n; i++)
1288 (*r)[i] = gmx::fromString<real>(values[i]);
1290 catch (gmx::GromacsException &)
1292 warning_error(wi, "Invalid value " + values[i] + " in string in mdp file. Expected a real number.");
1298 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1301 int i, j, max_n_lambda, nweights, nfep[efptNR];
1302 t_lambda *fep = ir->fepvals;
1303 t_expanded *expand = ir->expandedvals;
1304 real **count_fep_lambdas;
1305 bool bOneLambda = TRUE;
1307 snew(count_fep_lambdas, efptNR);
1309 /* FEP input processing */
1310 /* first, identify the number of lambda values for each type.
1311 All that are nonzero must have the same number */
1313 for (i = 0; i < efptNR; i++)
1315 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1318 /* now, determine the number of components. All must be either zero, or equal. */
1321 for (i = 0; i < efptNR; i++)
1323 if (nfep[i] > max_n_lambda)
1325 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1326 must have the same number if its not zero.*/
1331 for (i = 0; i < efptNR; i++)
1335 ir->fepvals->separate_dvdl[i] = FALSE;
1337 else if (nfep[i] == max_n_lambda)
1339 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1340 respect to the temperature currently */
1342 ir->fepvals->separate_dvdl[i] = TRUE;
1347 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1348 nfep[i], efpt_names[i], max_n_lambda);
1351 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1352 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1354 /* the number of lambdas is the number we've read in, which is either zero
1355 or the same for all */
1356 fep->n_lambda = max_n_lambda;
1358 /* allocate space for the array of lambda values */
1359 snew(fep->all_lambda, efptNR);
1360 /* if init_lambda is defined, we need to set lambda */
1361 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1363 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1365 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1366 for (i = 0; i < efptNR; i++)
1368 snew(fep->all_lambda[i], fep->n_lambda);
1369 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1372 for (j = 0; j < fep->n_lambda; j++)
1374 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1376 sfree(count_fep_lambdas[i]);
1379 sfree(count_fep_lambdas);
1381 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1382 bookkeeping -- for now, init_lambda */
1384 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1386 for (i = 0; i < fep->n_lambda; i++)
1388 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1392 /* check to see if only a single component lambda is defined, and soft core is defined.
1393 In this case, turn on coulomb soft core */
1395 if (max_n_lambda == 0)
1401 for (i = 0; i < efptNR; i++)
1403 if ((nfep[i] != 0) && (i != efptFEP))
1409 if ((bOneLambda) && (fep->sc_alpha > 0))
1411 fep->bScCoul = TRUE;
1414 /* Fill in the others with the efptFEP if they are not explicitly
1415 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1416 they are all zero. */
1418 for (i = 0; i < efptNR; i++)
1420 if ((nfep[i] == 0) && (i != efptFEP))
1422 for (j = 0; j < fep->n_lambda; j++)
1424 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1430 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1431 if (fep->sc_r_power == 48)
1433 if (fep->sc_alpha > 0.1)
1435 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1439 /* now read in the weights */
1440 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1443 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1445 else if (nweights != fep->n_lambda)
1447 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1448 nweights, fep->n_lambda);
1450 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1452 expand->nstexpanded = fep->nstdhdl;
1453 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1458 static void do_simtemp_params(t_inputrec *ir)
1461 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1462 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1466 convertYesNos(warninp_t /*wi*/, gmx::ArrayRef<const std::string> inputs, const char * /*name*/, gmx_bool *outputs)
1469 for (const auto &input : inputs)
1471 outputs[i] = gmx::equalCaseInsensitive(input, "Y", 1);
1476 template <typename T> void
1477 convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, T *outputs)
1480 for (const auto &input : inputs)
1484 outputs[i] = gmx::fromStdString<T>(input);
1486 catch (gmx::GromacsException &)
1488 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of integers separated by spaces.",
1490 warning_error(wi, message);
1497 convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, real *outputs)
1500 for (const auto &input : inputs)
1504 outputs[i] = gmx::fromString<real>(input);
1506 catch (gmx::GromacsException &)
1508 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of real numbers separated by spaces.",
1510 warning_error(wi, message);
1517 convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, rvec *outputs)
1520 for (const auto &input : inputs)
1524 outputs[i][d] = gmx::fromString<real>(input);
1526 catch (gmx::GromacsException &)
1528 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of real numbers separated by spaces.",
1530 warning_error(wi, message);
1541 static void do_wall_params(t_inputrec *ir,
1542 char *wall_atomtype, char *wall_density,
1546 opts->wall_atomtype[0] = nullptr;
1547 opts->wall_atomtype[1] = nullptr;
1549 ir->wall_atomtype[0] = -1;
1550 ir->wall_atomtype[1] = -1;
1551 ir->wall_density[0] = 0;
1552 ir->wall_density[1] = 0;
1556 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1557 if (wallAtomTypes.size() != size_t(ir->nwall))
1559 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu",
1560 ir->nwall, wallAtomTypes.size());
1562 for (int i = 0; i < ir->nwall; i++)
1564 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1567 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1569 auto wallDensity = gmx::splitString(wall_density);
1570 if (wallDensity.size() != size_t(ir->nwall))
1572 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall, wallDensity.size());
1574 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1575 for (int i = 0; i < ir->nwall; i++)
1577 if (ir->wall_density[i] <= 0)
1579 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1586 static void add_wall_energrps(SimulationGroups *groups, int nwall, t_symtab *symtab)
1590 AtomGroupIndices *grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1591 for (int i = 0; i < nwall; i++)
1593 groups->groupNames.emplace_back(
1596 gmx::formatString("wall%d", i).c_str()));
1597 grps->emplace_back(groups->groupNames.size() - 1);
1602 static void read_expandedparams(std::vector<t_inpfile> *inp,
1603 t_expanded *expand, warninp_t wi)
1605 /* read expanded ensemble parameters */
1606 printStringNewline(inp, "expanded ensemble variables");
1607 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1608 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1609 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1610 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1611 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1612 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1613 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1614 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1615 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1616 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1617 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1618 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1619 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1620 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1621 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1622 expand->bSymmetrizedTMatrix = (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1623 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1624 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1625 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1626 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1627 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1628 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1629 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1632 /*! \brief Return whether an end state with the given coupling-lambda
1633 * value describes fully-interacting VDW.
1635 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1636 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1638 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1640 return (couple_lambda_value == ecouplamVDW ||
1641 couple_lambda_value == ecouplamVDWQ);
1647 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1650 explicit MdpErrorHandler(warninp_t wi)
1651 : wi_(wi), mapping_(nullptr)
1655 void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
1657 mapping_ = &mapping;
1660 bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context) override
1662 ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
1663 getOptionName(context).c_str()));
1664 std::string message = gmx::formatExceptionMessageToString(*ex);
1665 warning_error(wi_, message.c_str());
1670 std::string getOptionName(const gmx::KeyValueTreePath &context)
1672 if (mapping_ != nullptr)
1674 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1675 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1678 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1683 const gmx::IKeyValueTreeBackMapping *mapping_;
1688 void get_ir(const char *mdparin, const char *mdparout,
1689 gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
1690 WriteMdpHeader writeMdpHeader, warninp_t wi)
1693 double dumdub[2][6];
1695 char warn_buf[STRLEN];
1696 t_lambda *fep = ir->fepvals;
1697 t_expanded *expand = ir->expandedvals;
1699 const char *no_names[] = { "no", nullptr };
1701 init_inputrec_strings();
1702 gmx::TextInputFile stream(mdparin);
1703 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1705 snew(dumstr[0], STRLEN);
1706 snew(dumstr[1], STRLEN);
1708 if (-1 == search_einp(inp, "cutoff-scheme"))
1711 "%s did not specify a value for the .mdp option "
1712 "\"cutoff-scheme\". As of GROMACS 2020, the Verlet scheme "
1713 "is the only cutoff scheme supported. This may affect your "
1714 "simulation if you are using an old mdp file that assumes use "
1715 "of the (removed) group cutoff scheme.", mdparin);
1716 warning_note(wi, warn_buf);
1719 /* ignore the following deprecated commands */
1720 replace_inp_entry(inp, "title", nullptr);
1721 replace_inp_entry(inp, "cpp", nullptr);
1722 replace_inp_entry(inp, "domain-decomposition", nullptr);
1723 replace_inp_entry(inp, "andersen-seed", nullptr);
1724 replace_inp_entry(inp, "dihre", nullptr);
1725 replace_inp_entry(inp, "dihre-fc", nullptr);
1726 replace_inp_entry(inp, "dihre-tau", nullptr);
1727 replace_inp_entry(inp, "nstdihreout", nullptr);
1728 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1729 replace_inp_entry(inp, "optimize-fft", nullptr);
1730 replace_inp_entry(inp, "adress_type", nullptr);
1731 replace_inp_entry(inp, "adress_const_wf", nullptr);
1732 replace_inp_entry(inp, "adress_ex_width", nullptr);
1733 replace_inp_entry(inp, "adress_hy_width", nullptr);
1734 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1735 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1736 replace_inp_entry(inp, "adress_site", nullptr);
1737 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1738 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1739 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1740 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1741 replace_inp_entry(inp, "rlistlong", nullptr);
1742 replace_inp_entry(inp, "nstcalclr", nullptr);
1743 replace_inp_entry(inp, "pull-print-com2", nullptr);
1744 replace_inp_entry(inp, "gb-algorithm", nullptr);
1745 replace_inp_entry(inp, "nstgbradii", nullptr);
1746 replace_inp_entry(inp, "rgbradii", nullptr);
1747 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1748 replace_inp_entry(inp, "gb-saltconc", nullptr);
1749 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1750 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1751 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1752 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1753 replace_inp_entry(inp, "sa-algorithm", nullptr);
1754 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1756 /* replace the following commands with the clearer new versions*/
1757 replace_inp_entry(inp, "unconstrained-start", "continuation");
1758 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1759 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1760 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1761 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1762 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1763 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1765 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1766 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1767 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1768 setStringEntry(&inp, "include", opts->include, nullptr);
1769 printStringNoNewline(&inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1770 setStringEntry(&inp, "define", opts->define, nullptr);
1772 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1773 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1774 printStringNoNewline(&inp, "Start time and timestep in ps");
1775 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1776 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1777 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1778 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1779 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1780 printStringNoNewline(&inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1781 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1782 printStringNoNewline(&inp, "mode for center of mass motion removal");
1783 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
1784 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
1785 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
1786 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
1787 setStringEntry(&inp, "comm-grps", is->vcm, nullptr);
1789 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
1790 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
1791 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
1792 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
1795 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
1796 printStringNoNewline(&inp, "Force tolerance and initial step-size");
1797 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
1798 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
1799 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
1800 ir->niter = get_eint(&inp, "niter", 20, wi);
1801 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
1802 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
1803 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
1804 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
1805 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
1807 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
1808 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
1810 /* Output options */
1811 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
1812 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
1813 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
1814 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
1815 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
1816 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
1817 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
1818 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
1819 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
1820 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
1821 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
1822 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
1823 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
1824 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
1825 printStringNoNewline(&inp, "default, all atoms will be written.");
1826 setStringEntry(&inp, "compressed-x-grps", is->x_compressed_groups, nullptr);
1827 printStringNoNewline(&inp, "Selection of energy groups");
1828 setStringEntry(&inp, "energygrps", is->energy, nullptr);
1830 /* Neighbor searching */
1831 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
1832 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1833 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
1834 printStringNoNewline(&inp, "nblist update frequency");
1835 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
1836 printStringNoNewline(&inp, "ns algorithm (simple or grid)");
1837 ir->ns_type = get_eeenum(&inp, "ns-type", ens_names, wi);
1838 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
1839 ir->ePBC = get_eeenum(&inp, "pbc", epbc_names, wi);
1840 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
1841 printStringNoNewline(&inp, "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1842 printStringNoNewline(&inp, "a value of -1 means: use rlist");
1843 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
1844 printStringNoNewline(&inp, "nblist cut-off");
1845 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
1846 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
1848 /* Electrostatics */
1849 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
1850 printStringNoNewline(&inp, "Method for doing electrostatics");
1851 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
1852 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
1853 printStringNoNewline(&inp, "cut-off lengths");
1854 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
1855 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
1856 printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
1857 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
1858 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
1859 printStringNoNewline(&inp, "Method for doing Van der Waals");
1860 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
1861 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
1862 printStringNoNewline(&inp, "cut-off lengths");
1863 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
1864 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
1865 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
1866 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
1867 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
1868 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
1869 printStringNoNewline(&inp, "Separate tables between energy group pairs");
1870 setStringEntry(&inp, "energygrp-table", is->egptable, nullptr);
1871 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
1872 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
1873 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
1874 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
1875 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
1876 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
1877 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
1878 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
1879 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
1880 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
1881 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
1882 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
1883 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
1885 /* Implicit solvation is no longer supported, but we need grompp
1886 to be able to refuse old .mdp files that would have built a tpr
1887 to run it. Thus, only "no" is accepted. */
1888 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
1890 /* Coupling stuff */
1891 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
1892 printStringNoNewline(&inp, "Temperature coupling");
1893 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
1894 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
1895 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
1896 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
1897 printStringNoNewline(&inp, "Groups to couple separately");
1898 setStringEntry(&inp, "tc-grps", is->tcgrps, nullptr);
1899 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
1900 setStringEntry(&inp, "tau-t", is->tau_t, nullptr);
1901 setStringEntry(&inp, "ref-t", is->ref_t, nullptr);
1902 printStringNoNewline(&inp, "pressure coupling");
1903 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
1904 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
1905 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
1906 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
1907 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
1908 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
1909 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
1910 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
1911 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
1914 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
1915 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
1916 printStringNoNewline(&inp, "Groups treated Quantum Mechanically");
1917 setStringEntry(&inp, "QMMM-grps", is->QMMM, nullptr);
1918 printStringNoNewline(&inp, "QM method");
1919 setStringEntry(&inp, "QMmethod", is->QMmethod, nullptr);
1920 printStringNoNewline(&inp, "QMMM scheme");
1921 ir->QMMMscheme = get_eeenum(&inp, "QMMMscheme", eQMMMscheme_names, wi);
1922 printStringNoNewline(&inp, "QM basisset");
1923 setStringEntry(&inp, "QMbasis", is->QMbasis, nullptr);
1924 printStringNoNewline(&inp, "QM charge");
1925 setStringEntry(&inp, "QMcharge", is->QMcharge, nullptr);
1926 printStringNoNewline(&inp, "QM multiplicity");
1927 setStringEntry(&inp, "QMmult", is->QMmult, nullptr);
1928 printStringNoNewline(&inp, "Surface Hopping");
1929 setStringEntry(&inp, "SH", is->bSH, nullptr);
1930 printStringNoNewline(&inp, "CAS space options");
1931 setStringEntry(&inp, "CASorbitals", is->CASorbitals, nullptr);
1932 setStringEntry(&inp, "CASelectrons", is->CASelectrons, nullptr);
1933 setStringEntry(&inp, "SAon", is->SAon, nullptr);
1934 setStringEntry(&inp, "SAoff", is->SAoff, nullptr);
1935 setStringEntry(&inp, "SAsteps", is->SAsteps, nullptr);
1936 printStringNoNewline(&inp, "Scale factor for MM charges");
1937 ir->scalefactor = get_ereal(&inp, "MMChargeScaleFactor", 1.0, wi);
1939 /* Simulated annealing */
1940 printStringNewline(&inp, "SIMULATED ANNEALING");
1941 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
1942 setStringEntry(&inp, "annealing", is->anneal, nullptr);
1943 printStringNoNewline(&inp, "Number of time points to use for specifying annealing in each group");
1944 setStringEntry(&inp, "annealing-npoints", is->anneal_npoints, nullptr);
1945 printStringNoNewline(&inp, "List of times at the annealing points for each group");
1946 setStringEntry(&inp, "annealing-time", is->anneal_time, nullptr);
1947 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
1948 setStringEntry(&inp, "annealing-temp", is->anneal_temp, nullptr);
1951 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
1952 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
1953 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
1954 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
1957 printStringNewline(&inp, "OPTIONS FOR BONDS");
1958 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
1959 printStringNoNewline(&inp, "Type of constraint algorithm");
1960 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
1961 printStringNoNewline(&inp, "Do not constrain the start configuration");
1962 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
1963 printStringNoNewline(&inp, "Use successive overrelaxation to reduce the number of shake iterations");
1964 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
1965 printStringNoNewline(&inp, "Relative tolerance of shake");
1966 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
1967 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
1968 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
1969 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
1970 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
1971 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
1972 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
1973 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
1974 printStringNoNewline(&inp, "rotates over more degrees than");
1975 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
1976 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
1977 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
1979 /* Energy group exclusions */
1980 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
1981 printStringNoNewline(&inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
1982 setStringEntry(&inp, "energygrp-excl", is->egpexcl, nullptr);
1985 printStringNewline(&inp, "WALLS");
1986 printStringNoNewline(&inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1987 ir->nwall = get_eint(&inp, "nwall", 0, wi);
1988 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
1989 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
1990 setStringEntry(&inp, "wall-atomtype", is->wall_atomtype, nullptr);
1991 setStringEntry(&inp, "wall-density", is->wall_density, nullptr);
1992 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
1995 printStringNewline(&inp, "COM PULLING");
1996 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2000 is->pull_grp = read_pullparams(&inp, ir->pull, wi);
2004 NOTE: needs COM pulling input */
2005 printStringNewline(&inp, "AWH biasing");
2006 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2011 ir->awhParams = gmx::readAndCheckAwhParams(&inp, ir, wi);
2015 gmx_fatal(FARGS, "AWH biasing is only compatible with COM pulling turned on");
2019 /* Enforced rotation */
2020 printStringNewline(&inp, "ENFORCED ROTATION");
2021 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2022 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2026 is->rot_grp = read_rotparams(&inp, ir->rot, wi);
2029 /* Interactive MD */
2031 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2032 setStringEntry(&inp, "IMD-group", is->imd_grp, nullptr);
2033 if (is->imd_grp[0] != '\0')
2040 printStringNewline(&inp, "NMR refinement stuff");
2041 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2042 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2043 printStringNoNewline(&inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2044 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2045 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2046 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2047 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2048 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2049 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2050 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2051 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2052 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2053 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2054 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2055 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2056 setStringEntry(&inp, "orire-fitgrp", is->orirefitgrp, nullptr);
2057 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2058 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2060 /* free energy variables */
2061 printStringNewline(&inp, "Free energy variables");
2062 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2063 setStringEntry(&inp, "couple-moltype", is->couple_moltype, nullptr);
2064 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2065 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2066 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2068 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2070 it was not entered */
2071 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2072 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2073 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2074 setStringEntry(&inp, "fep-lambdas", is->fep_lambda[efptFEP], nullptr);
2075 setStringEntry(&inp, "mass-lambdas", is->fep_lambda[efptMASS], nullptr);
2076 setStringEntry(&inp, "coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
2077 setStringEntry(&inp, "vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
2078 setStringEntry(&inp, "bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
2079 setStringEntry(&inp, "restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
2080 setStringEntry(&inp, "temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
2081 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2082 setStringEntry(&inp, "init-lambda-weights", is->lambda_weights, nullptr);
2083 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2084 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2085 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2086 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2087 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2088 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2089 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2090 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2091 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2092 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2093 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2094 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2096 /* Non-equilibrium MD stuff */
2097 printStringNewline(&inp, "Non-equilibrium MD stuff");
2098 setStringEntry(&inp, "acc-grps", is->accgrps, nullptr);
2099 setStringEntry(&inp, "accelerate", is->acc, nullptr);
2100 setStringEntry(&inp, "freezegrps", is->freeze, nullptr);
2101 setStringEntry(&inp, "freezedim", is->frdim, nullptr);
2102 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2103 setStringEntry(&inp, "deform", is->deform, nullptr);
2105 /* simulated tempering variables */
2106 printStringNewline(&inp, "simulated tempering variables");
2107 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2108 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2109 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2110 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2112 /* expanded ensemble variables */
2113 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2115 read_expandedparams(&inp, expand, wi);
2118 /* Electric fields */
2120 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2121 gmx::KeyValueTreeTransformer transform;
2122 transform.rules()->addRule()
2123 .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2124 mdModules->initMdpTransform(transform.rules());
2125 for (const auto &path : transform.mappedPaths())
2127 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2128 mark_einp_set(inp, path[0].c_str());
2130 MdpErrorHandler errorHandler(wi);
2132 = transform.transform(convertedValues, &errorHandler);
2133 ir->params = new gmx::KeyValueTreeObject(result.object());
2134 mdModules->adjustInputrecBasedOnModules(ir);
2135 errorHandler.setBackMapping(result.backMapping());
2136 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2139 /* Ion/water position swapping ("computational electrophysiology") */
2140 printStringNewline(&inp, "Ion/water position swapping for computational electrophysiology setups");
2141 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2142 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2143 if (ir->eSwapCoords != eswapNO)
2150 printStringNoNewline(&inp, "Swap attempt frequency");
2151 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2152 printStringNoNewline(&inp, "Number of ion types to be controlled");
2153 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2156 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2158 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2159 snew(ir->swap->grp, ir->swap->ngrp);
2160 for (i = 0; i < ir->swap->ngrp; i++)
2162 snew(ir->swap->grp[i].molname, STRLEN);
2164 printStringNoNewline(&inp, "Two index groups that contain the compartment-partitioning atoms");
2165 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2166 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2167 printStringNoNewline(&inp, "Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2168 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2169 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2171 printStringNoNewline(&inp, "Name of solvent molecules");
2172 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2174 printStringNoNewline(&inp, "Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2175 printStringNoNewline(&inp, "Note that the split cylinder settings do not have an influence on the swapping protocol,");
2176 printStringNoNewline(&inp, "however, if correctly defined, the permeation events are recorded per channel");
2177 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2178 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2179 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2180 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2181 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2182 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2184 printStringNoNewline(&inp, "Average the number of ions per compartment over these many swap attempt steps");
2185 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2187 printStringNoNewline(&inp, "Names of the ion types that can be exchanged with solvent molecules,");
2188 printStringNoNewline(&inp, "and the requested number of ions of this type in compartments A and B");
2189 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2190 for (i = 0; i < nIonTypes; i++)
2192 int ig = eSwapFixedGrpNR + i;
2194 sprintf(buf, "iontype%d-name", i);
2195 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2196 sprintf(buf, "iontype%d-in-A", i);
2197 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2198 sprintf(buf, "iontype%d-in-B", i);
2199 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2202 printStringNoNewline(&inp, "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2203 printStringNoNewline(&inp, "at maximum distance (= bulk concentration) to the split group layers. However,");
2204 printStringNoNewline(&inp, "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2205 printStringNoNewline(&inp, "towards one of the compartment-partitioning layers (at +/- 1.0).");
2206 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2207 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2208 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2209 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2211 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2214 printStringNoNewline(&inp, "Start to swap ions if threshold difference to requested count is reached");
2215 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2218 /* AdResS is no longer supported, but we need grompp to be able to
2219 refuse to process old .mdp files that used it. */
2220 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2222 /* User defined thingies */
2223 printStringNewline(&inp, "User defined thingies");
2224 setStringEntry(&inp, "user1-grps", is->user1, nullptr);
2225 setStringEntry(&inp, "user2-grps", is->user2, nullptr);
2226 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2227 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2228 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2229 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2230 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2231 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2232 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2233 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2237 gmx::TextOutputFile stream(mdparout);
2238 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2240 // Transform module data into a flat key-value tree for output.
2241 gmx::KeyValueTreeBuilder builder;
2242 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2243 mdModules->buildMdpOutput(&builderObject);
2245 gmx::TextWriter writer(&stream);
2246 writeKeyValueTreeAsMdp(&writer, builder.build());
2251 /* Process options if necessary */
2252 for (m = 0; m < 2; m++)
2254 for (i = 0; i < 2*DIM; i++)
2263 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2265 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
2267 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2269 case epctSEMIISOTROPIC:
2270 case epctSURFACETENSION:
2271 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2273 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
2275 dumdub[m][YY] = dumdub[m][XX];
2277 case epctANISOTROPIC:
2278 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2279 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2280 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2282 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
2286 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2287 epcoupltype_names[ir->epct]);
2291 clear_mat(ir->ref_p);
2292 clear_mat(ir->compress);
2293 for (i = 0; i < DIM; i++)
2295 ir->ref_p[i][i] = dumdub[1][i];
2296 ir->compress[i][i] = dumdub[0][i];
2298 if (ir->epct == epctANISOTROPIC)
2300 ir->ref_p[XX][YY] = dumdub[1][3];
2301 ir->ref_p[XX][ZZ] = dumdub[1][4];
2302 ir->ref_p[YY][ZZ] = dumdub[1][5];
2303 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2305 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2307 ir->compress[XX][YY] = dumdub[0][3];
2308 ir->compress[XX][ZZ] = dumdub[0][4];
2309 ir->compress[YY][ZZ] = dumdub[0][5];
2310 for (i = 0; i < DIM; i++)
2312 for (m = 0; m < i; m++)
2314 ir->ref_p[i][m] = ir->ref_p[m][i];
2315 ir->compress[i][m] = ir->compress[m][i];
2320 if (ir->comm_mode == ecmNO)
2325 opts->couple_moltype = nullptr;
2326 if (strlen(is->couple_moltype) > 0)
2328 if (ir->efep != efepNO)
2330 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2331 if (opts->couple_lam0 == opts->couple_lam1)
2333 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2335 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2336 opts->couple_lam1 == ecouplamNONE))
2338 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2343 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2346 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2347 if (ir->efep != efepNO)
2349 if (fep->delta_lambda > 0)
2351 ir->efep = efepSLOWGROWTH;
2355 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2357 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2358 warning_note(wi, "Old option for dhdl-print-energy given: "
2359 "changing \"yes\" to \"total\"\n");
2362 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2364 /* always print out the energy to dhdl if we are doing
2365 expanded ensemble, since we need the total energy for
2366 analysis if the temperature is changing. In some
2367 conditions one may only want the potential energy, so
2368 we will allow that if the appropriate mdp setting has
2369 been enabled. Otherwise, total it is:
2371 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2374 if ((ir->efep != efepNO) || ir->bSimTemp)
2376 ir->bExpanded = FALSE;
2377 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2379 ir->bExpanded = TRUE;
2381 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2382 if (ir->bSimTemp) /* done after fep params */
2384 do_simtemp_params(ir);
2387 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2388 * setup and not on the old way of specifying the free-energy setup,
2389 * we should check for using soft-core when not needed, since that
2390 * can complicate the sampling significantly.
2391 * Note that we only check for the automated coupling setup.
2392 * If the (advanced) user does FEP through manual topology changes,
2393 * this check will not be triggered.
2395 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2396 ir->fepvals->sc_alpha != 0 &&
2397 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2398 couple_lambda_has_vdw_on(opts->couple_lam1)))
2400 warning(wi, "You are using soft-core interactions while the Van der Waals interactions are not decoupled (note that the sc-coul option is only active when using lambda states). Although this will not lead to errors, you will need much more sampling than without soft-core interactions. Consider using sc-alpha=0.");
2405 ir->fepvals->n_lambda = 0;
2408 /* WALL PARAMETERS */
2410 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts, wi);
2412 /* ORIENTATION RESTRAINT PARAMETERS */
2414 if (opts->bOrire && gmx::splitString(is->orirefitgrp).size() != 1)
2416 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2419 /* DEFORMATION PARAMETERS */
2421 clear_mat(ir->deform);
2422 for (i = 0; i < 6; i++)
2427 double gmx_unused canary;
2428 int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
2429 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2430 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
2432 if (strlen(is->deform) > 0 && ndeform != 6)
2434 warning_error(wi, gmx::formatString("Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform).c_str());
2436 for (i = 0; i < 3; i++)
2438 ir->deform[i][i] = dumdub[0][i];
2440 ir->deform[YY][XX] = dumdub[0][3];
2441 ir->deform[ZZ][XX] = dumdub[0][4];
2442 ir->deform[ZZ][YY] = dumdub[0][5];
2443 if (ir->epc != epcNO)
2445 for (i = 0; i < 3; i++)
2447 for (j = 0; j <= i; j++)
2449 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2451 warning_error(wi, "A box element has deform set and compressibility > 0");
2455 for (i = 0; i < 3; i++)
2457 for (j = 0; j < i; j++)
2459 if (ir->deform[i][j] != 0)
2461 for (m = j; m < DIM; m++)
2463 if (ir->compress[m][j] != 0)
2465 sprintf(warn_buf, "An off-diagonal box element has deform set while compressibility > 0 for the same component of another box vector, this might lead to spurious periodicity effects.");
2466 warning(wi, warn_buf);
2474 /* Ion/water position swapping checks */
2475 if (ir->eSwapCoords != eswapNO)
2477 if (ir->swap->nstswap < 1)
2479 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2481 if (ir->swap->nAverage < 1)
2483 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2485 if (ir->swap->threshold < 1.0)
2487 warning_error(wi, "Ion count threshold must be at least 1.\n");
2495 static int search_QMstring(const char *s, int ng, const char *gn[])
2497 /* same as normal search_string, but this one searches QM strings */
2500 for (i = 0; (i < ng); i++)
2502 if (gmx_strcasecmp(s, gn[i]) == 0)
2508 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2509 } /* search_QMstring */
2511 /* We would like gn to be const as well, but C doesn't allow this */
2512 /* TODO this is utility functionality (search for the index of a
2513 string in a collection), so should be refactored and located more
2515 int search_string(const char *s, int ng, char *gn[])
2519 for (i = 0; (i < ng); i++)
2521 if (gmx_strcasecmp(s, gn[i]) == 0)
2528 "Group %s referenced in the .mdp file was not found in the index file.\n"
2529 "Group names must match either [moleculetype] names or custom index group\n"
2530 "names, in which case you must supply an index file to the '-n' option\n"
2535 static bool do_numbering(int natoms, SimulationGroups *groups,
2536 gmx::ArrayRef<std::string> groupsFromMdpFile,
2537 t_blocka *block, char *gnames[],
2538 SimulationAtomGroupType gtype, int restnm,
2539 int grptp, bool bVerbose,
2542 unsigned short *cbuf;
2543 AtomGroupIndices *grps = &(groups->groups[gtype]);
2544 int j, gid, aj, ognr, ntot = 0;
2547 char warn_buf[STRLEN];
2549 title = shortName(gtype);
2552 /* Mark all id's as not set */
2553 for (int i = 0; (i < natoms); i++)
2558 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2560 /* Lookup the group name in the block structure */
2561 gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2562 if ((grptp != egrptpONE) || (i == 0))
2564 grps->emplace_back(gid);
2567 /* Now go over the atoms in the group */
2568 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2573 /* Range checking */
2574 if ((aj < 0) || (aj >= natoms))
2576 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2578 /* Lookup up the old group number */
2582 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2583 aj+1, title, ognr+1, i+1);
2587 /* Store the group number in buffer */
2588 if (grptp == egrptpONE)
2601 /* Now check whether we have done all atoms */
2605 if (grptp == egrptpALL)
2607 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2608 natoms-ntot, title);
2610 else if (grptp == egrptpPART)
2612 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2613 natoms-ntot, title);
2614 warning_note(wi, warn_buf);
2616 /* Assign all atoms currently unassigned to a rest group */
2617 for (j = 0; (j < natoms); j++)
2619 if (cbuf[j] == NOGID)
2621 cbuf[j] = grps->size();
2625 if (grptp != egrptpPART)
2630 "Making dummy/rest group for %s containing %d elements\n",
2631 title, natoms-ntot);
2633 /* Add group name "rest" */
2634 grps->emplace_back(restnm);
2636 /* Assign the rest name to all atoms not currently assigned to a group */
2637 for (j = 0; (j < natoms); j++)
2639 if (cbuf[j] == NOGID)
2641 // group size was not updated before this here, so need to use -1.
2642 cbuf[j] = grps->size() - 1;
2648 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2650 /* All atoms are part of one (or no) group, no index required */
2651 groups->groupNumbers[gtype].clear();
2655 for (int j = 0; (j < natoms); j++)
2657 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2663 return (bRest && grptp == egrptpPART);
2666 static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2669 pull_params_t *pull;
2670 int natoms, imin, jmin;
2671 int *nrdf2, *na_vcm, na_tot;
2672 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2677 * First calc 3xnr-atoms for each group
2678 * then subtract half a degree of freedom for each constraint
2680 * Only atoms and nuclei contribute to the degrees of freedom...
2685 const SimulationGroups &groups = mtop->groups;
2686 natoms = mtop->natoms;
2688 /* Allocate one more for a possible rest group */
2689 /* We need to sum degrees of freedom into doubles,
2690 * since floats give too low nrdf's above 3 million atoms.
2692 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size()+1);
2693 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2694 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2695 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2696 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2698 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2702 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; i++)
2705 clear_ivec(dof_vcm[i]);
2707 nrdf_vcm_sub[i] = 0;
2709 snew(nrdf2, natoms);
2710 for (const AtomProxy atomP : AtomRange(*mtop))
2712 const t_atom &local = atomP.atom();
2713 int i = atomP.globalAtomNumber();
2715 if (local.ptype == eptAtom || local.ptype == eptNucleus)
2717 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
2718 for (int d = 0; d < DIM; d++)
2720 if (opts->nFreeze[g][d] == 0)
2722 /* Add one DOF for particle i (counted as 2*1) */
2724 /* VCM group i has dim d as a DOF */
2725 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] = 1;
2728 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] += 0.5*nrdf2[i];
2729 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] += 0.5*nrdf2[i];
2734 for (const gmx_molblock_t &molb : mtop->molblock)
2736 const gmx_moltype_t &molt = mtop->moltype[molb.type];
2737 const t_atom *atom = molt.atoms.atom;
2738 for (int mol = 0; mol < molb.nmol; mol++)
2740 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2742 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
2743 for (int i = 0; i < molt.ilist[ftype].size(); )
2745 /* Subtract degrees of freedom for the constraints,
2746 * if the particles still have degrees of freedom left.
2747 * If one of the particles is a vsite or a shell, then all
2748 * constraint motion will go there, but since they do not
2749 * contribute to the constraints the degrees of freedom do not
2752 int ai = as + ia[i + 1];
2753 int aj = as + ia[i + 2];
2754 if (((atom[ia[i + 1]].ptype == eptNucleus) ||
2755 (atom[ia[i + 1]].ptype == eptAtom)) &&
2756 ((atom[ia[i + 2]].ptype == eptNucleus) ||
2757 (atom[ia[i + 2]].ptype == eptAtom)))
2775 imin = std::min(imin, nrdf2[ai]);
2776 jmin = std::min(jmin, nrdf2[aj]);
2779 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2780 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -= 0.5*jmin;
2781 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2782 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -= 0.5*jmin;
2784 i += interaction_function[ftype].nratoms+1;
2787 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
2788 for (int i = 0; i < molt.ilist[F_SETTLE].size(); )
2790 /* Subtract 1 dof from every atom in the SETTLE */
2791 for (int j = 0; j < 3; j++)
2793 int ai = as + ia[i + 1 + j];
2794 imin = std::min(2, nrdf2[ai]);
2796 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2797 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2801 as += molt.atoms.nr;
2807 /* Correct nrdf for the COM constraints.
2808 * We correct using the TC and VCM group of the first atom
2809 * in the reference and pull group. If atoms in one pull group
2810 * belong to different TC or VCM groups it is anyhow difficult
2811 * to determine the optimal nrdf assignment.
2815 for (int i = 0; i < pull->ncoord; i++)
2817 if (pull->coord[i].eType != epullCONSTRAINT)
2824 for (int j = 0; j < 2; j++)
2826 const t_pull_group *pgrp;
2828 pgrp = &pull->group[pull->coord[i].group[j]];
2832 /* Subtract 1/2 dof from each group */
2833 int ai = pgrp->ind[0];
2834 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2835 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2836 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
2838 gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)]]);
2843 /* We need to subtract the whole DOF from group j=1 */
2850 if (ir->nstcomm != 0)
2854 /* We remove COM motion up to dim ndof_com() */
2855 ndim_rm_vcm = ndof_com(ir);
2857 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
2858 * the number of degrees of freedom in each vcm group when COM
2859 * translation is removed and 6 when rotation is removed as well.
2861 for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2863 switch (ir->comm_mode)
2866 case ecmLINEAR_ACCELERATION_CORRECTION:
2867 nrdf_vcm_sub[j] = 0;
2868 for (int d = 0; d < ndim_rm_vcm; d++)
2877 nrdf_vcm_sub[j] = 6;
2880 gmx_incons("Checking comm_mode");
2884 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2886 /* Count the number of atoms of TC group i for every VCM group */
2887 for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2892 for (int ai = 0; ai < natoms; ai++)
2894 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
2896 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
2900 /* Correct for VCM removal according to the fraction of each VCM
2901 * group present in this TC group.
2903 nrdf_uc = nrdf_tc[i];
2905 for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2907 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
2909 nrdf_tc[i] += nrdf_uc*(static_cast<double>(na_vcm[j])/static_cast<double>(na_tot))*
2910 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
2915 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
2917 opts->nrdf[i] = nrdf_tc[i];
2918 if (opts->nrdf[i] < 0)
2923 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2924 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
2932 sfree(nrdf_vcm_sub);
2935 static bool do_egp_flag(t_inputrec *ir, SimulationGroups *groups,
2936 const char *option, const char *val, int flag)
2938 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2939 * But since this is much larger than STRLEN, such a line can not be parsed.
2940 * The real maximum is the number of names that fit in a string: STRLEN/2.
2942 #define EGP_MAX (STRLEN/2)
2946 auto names = gmx::splitString(val);
2947 if (names.size() % 2 != 0)
2949 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2951 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
2953 for (size_t i = 0; i < names.size() / 2; i++)
2955 // TODO this needs to be replaced by a solution using std::find_if
2958 gmx_strcasecmp(names[2*i].c_str(), *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
2964 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2965 names[2*i].c_str(), option);
2969 gmx_strcasecmp(names[2*i+1].c_str(), *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
2975 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2976 names[2*i+1].c_str(), option);
2978 if ((j < nr) && (k < nr))
2980 ir->opts.egp_flags[nr*j+k] |= flag;
2981 ir->opts.egp_flags[nr*k+j] |= flag;
2990 static void make_swap_groups(
2995 int ig = -1, i = 0, gind;
2999 /* Just a quick check here, more thorough checks are in mdrun */
3000 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3002 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3005 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3006 for (ig = 0; ig < swap->ngrp; ig++)
3008 swapg = &swap->grp[ig];
3009 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3010 swapg->nat = grps->index[gind+1] - grps->index[gind];
3014 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3015 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3016 swap->grp[ig].molname, swapg->nat);
3017 snew(swapg->ind, swapg->nat);
3018 for (i = 0; i < swapg->nat; i++)
3020 swapg->ind[i] = grps->a[grps->index[gind]+i];
3025 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3031 static void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3036 ig = search_string(IMDgname, grps->nr, gnames);
3037 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3039 if (IMDgroup->nat > 0)
3041 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3042 IMDgname, IMDgroup->nat);
3043 snew(IMDgroup->ind, IMDgroup->nat);
3044 for (i = 0; i < IMDgroup->nat; i++)
3046 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3051 void do_index(const char* mdparin, const char *ndx,
3054 const gmx::MDModules::notifier_type ¬ifier,
3058 t_blocka *defaultIndexGroups;
3062 char warnbuf[STRLEN], **gnames;
3066 int i, j, k, restnm;
3067 bool bExcl, bTable, bAnneal, bRest;
3068 char warn_buf[STRLEN];
3072 fprintf(stderr, "processing index file...\n");
3076 snew(defaultIndexGroups, 1);
3077 snew(defaultIndexGroups->index, 1);
3079 atoms_all = gmx_mtop_global_atoms(mtop);
3080 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3081 done_atom(&atoms_all);
3085 defaultIndexGroups = init_index(ndx, &gnames);
3088 SimulationGroups *groups = &mtop->groups;
3089 natoms = mtop->natoms;
3090 symtab = &mtop->symtab;
3092 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3094 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3096 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3097 restnm = groups->groupNames.size() - 1;
3098 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3099 srenew(gnames, defaultIndexGroups->nr+1);
3100 gnames[restnm] = *(groups->groupNames.back());
3102 set_warning_line(wi, mdparin, -1);
3104 auto temperatureCouplingTauValues = gmx::splitString(is->tau_t);
3105 auto temperatureCouplingReferenceValues = gmx::splitString(is->ref_t);
3106 auto temperatureCouplingGroupNames = gmx::splitString(is->tcgrps);
3107 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size() ||
3108 temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3110 gmx_fatal(FARGS, "Invalid T coupling input: %zu groups, %zu ref-t values and "
3112 temperatureCouplingGroupNames.size(),
3113 temperatureCouplingReferenceValues.size(),
3114 temperatureCouplingTauValues.size());
3117 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3118 do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
3119 SimulationAtomGroupType::TemperatureCoupling,
3120 restnm, useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3121 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3123 snew(ir->opts.nrdf, nr);
3124 snew(ir->opts.tau_t, nr);
3125 snew(ir->opts.ref_t, nr);
3126 if (ir->eI == eiBD && ir->bd_fric == 0)
3128 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3131 if (useReferenceTemperature)
3133 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3135 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3139 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3140 for (i = 0; (i < nr); i++)
3142 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3144 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3145 warning_error(wi, warn_buf);
3148 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3150 warning_note(wi, "tau-t = -1 is the value to signal that a group should not have temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3153 if (ir->opts.tau_t[i] >= 0)
3155 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3158 if (ir->etc != etcNO && ir->nsttcouple == -1)
3160 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3165 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3167 gmx_fatal(FARGS, "Cannot do Nose-Hoover temperature with Berendsen pressure control with md-vv; use either vrescale temperature with berendsen pressure or Nose-Hoover temperature with MTTK pressure");
3169 if (ir->epc == epcMTTK)
3171 if (ir->etc != etcNOSEHOOVER)
3173 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3177 if (ir->nstpcouple != ir->nsttcouple)
3179 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3180 ir->nstpcouple = ir->nsttcouple = mincouple;
3181 sprintf(warn_buf, "for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple) = %d", mincouple);
3182 warning_note(wi, warn_buf);
3187 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3188 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3190 if (ETC_ANDERSEN(ir->etc))
3192 if (ir->nsttcouple != 1)
3195 sprintf(warn_buf, "Andersen temperature control methods assume nsttcouple = 1; there is no need for larger nsttcouple > 1, since no global parameters are computed. nsttcouple has been reset to 1");
3196 warning_note(wi, warn_buf);
3199 nstcmin = tcouple_min_integration_steps(ir->etc);
3202 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3204 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3205 ETCOUPLTYPE(ir->etc),
3207 ir->nsttcouple*ir->delta_t);
3208 warning(wi, warn_buf);
3211 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3212 for (i = 0; (i < nr); i++)
3214 if (ir->opts.ref_t[i] < 0)
3216 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3219 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3220 if we are in this conditional) if mc_temp is negative */
3221 if (ir->expandedvals->mc_temp < 0)
3223 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3227 /* Simulated annealing for each group. There are nr groups */
3228 auto simulatedAnnealingGroupNames = gmx::splitString(is->anneal);
3229 if (simulatedAnnealingGroupNames.size() == 1 &&
3230 gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3232 simulatedAnnealingGroupNames.resize(0);
3234 if (!simulatedAnnealingGroupNames.empty() &&
3235 gmx::ssize(simulatedAnnealingGroupNames) != nr)
3237 gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
3238 simulatedAnnealingGroupNames.size(), nr);
3242 snew(ir->opts.annealing, nr);
3243 snew(ir->opts.anneal_npoints, nr);
3244 snew(ir->opts.anneal_time, nr);
3245 snew(ir->opts.anneal_temp, nr);
3246 for (i = 0; i < nr; i++)
3248 ir->opts.annealing[i] = eannNO;
3249 ir->opts.anneal_npoints[i] = 0;
3250 ir->opts.anneal_time[i] = nullptr;
3251 ir->opts.anneal_temp[i] = nullptr;
3253 if (!simulatedAnnealingGroupNames.empty())
3256 for (i = 0; i < nr; i++)
3258 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3260 ir->opts.annealing[i] = eannNO;
3262 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3264 ir->opts.annealing[i] = eannSINGLE;
3267 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3269 ir->opts.annealing[i] = eannPERIODIC;
3275 /* Read the other fields too */
3276 auto simulatedAnnealingPoints = gmx::splitString(is->anneal_npoints);
3277 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3279 gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
3280 simulatedAnnealingPoints.size(), simulatedAnnealingGroupNames.size());
3282 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3283 size_t numSimulatedAnnealingFields = 0;
3284 for (i = 0; i < nr; i++)
3286 if (ir->opts.anneal_npoints[i] == 1)
3288 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3290 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3291 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3292 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3295 auto simulatedAnnealingTimes = gmx::splitString(is->anneal_time);
3297 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3299 gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
3300 simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
3302 auto simulatedAnnealingTemperatures = gmx::splitString(is->anneal_temp);
3303 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3305 gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
3306 simulatedAnnealingTemperatures.size(), numSimulatedAnnealingFields);
3309 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3310 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3311 convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
3312 convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp", allSimulatedAnnealingTemperatures.data());
3313 for (i = 0, k = 0; i < nr; i++)
3315 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3317 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3318 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3321 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3323 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3329 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3331 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3332 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3335 if (ir->opts.anneal_temp[i][j] < 0)
3337 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3342 /* Print out some summary information, to make sure we got it right */
3343 for (i = 0; i < nr; i++)
3345 if (ir->opts.annealing[i] != eannNO)
3347 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3348 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3349 *(groups->groupNames[j]), eann_names[ir->opts.annealing[i]],
3350 ir->opts.anneal_npoints[i]);
3351 fprintf(stderr, "Time (ps) Temperature (K)\n");
3352 /* All terms except the last one */
3353 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3355 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3358 /* Finally the last one */
3359 j = ir->opts.anneal_npoints[i]-1;
3360 if (ir->opts.annealing[i] == eannSINGLE)
3362 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3366 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3367 if (std::fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3369 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3380 make_pull_groups(ir->pull, is->pull_grp, defaultIndexGroups, gnames);
3382 make_pull_coords(ir->pull);
3387 make_rotation_groups(ir->rot, is->rot_grp, defaultIndexGroups, gnames);
3390 if (ir->eSwapCoords != eswapNO)
3392 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3395 /* Make indices for IMD session */
3398 make_IMD_group(ir->imd, is->imd_grp, defaultIndexGroups, gnames);
3401 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3402 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3403 notifier.notify(defaultIndexGroupsAndNames);
3405 auto accelerations = gmx::splitString(is->acc);
3406 auto accelerationGroupNames = gmx::splitString(is->accgrps);
3407 if (accelerationGroupNames.size() * DIM != accelerations.size())
3409 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3410 accelerationGroupNames.size(), accelerations.size());
3412 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3413 SimulationAtomGroupType::Acceleration,
3414 restnm, egrptpALL_GENREST, bVerbose, wi);
3415 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3416 snew(ir->opts.acc, nr);
3417 ir->opts.ngacc = nr;
3419 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3421 auto freezeDims = gmx::splitString(is->frdim);
3422 auto freezeGroupNames = gmx::splitString(is->freeze);
3423 if (freezeDims.size() != DIM * freezeGroupNames.size())
3425 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3426 freezeGroupNames.size(), freezeDims.size());
3428 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3429 SimulationAtomGroupType::Freeze,
3430 restnm, egrptpALL_GENREST, bVerbose, wi);
3431 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3432 ir->opts.ngfrz = nr;
3433 snew(ir->opts.nFreeze, nr);
3434 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3436 for (j = 0; (j < DIM); j++, k++)
3438 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3439 if (!ir->opts.nFreeze[i][j])
3441 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3443 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3444 "(not %s)", freezeDims[k].c_str());
3445 warning(wi, warn_buf);
3450 for (; (i < nr); i++)
3452 for (j = 0; (j < DIM); j++)
3454 ir->opts.nFreeze[i][j] = 0;
3458 auto energyGroupNames = gmx::splitString(is->energy);
3459 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3460 SimulationAtomGroupType::EnergyOutput,
3461 restnm, egrptpALL_GENREST, bVerbose, wi);
3462 add_wall_energrps(groups, ir->nwall, symtab);
3463 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3464 auto vcmGroupNames = gmx::splitString(is->vcm);
3466 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3467 SimulationAtomGroupType::MassCenterVelocityRemoval,
3468 restnm, vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3471 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3472 "This may lead to artifacts.\n"
3473 "In most cases one should use one group for the whole system.");
3476 /* Now we have filled the freeze struct, so we can calculate NRDF */
3477 calc_nrdf(mtop, ir, gnames);
3479 auto user1GroupNames = gmx::splitString(is->user1);
3480 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3481 SimulationAtomGroupType::User1,
3482 restnm, egrptpALL_GENREST, bVerbose, wi);
3483 auto user2GroupNames = gmx::splitString(is->user2);
3484 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3485 SimulationAtomGroupType::User2,
3486 restnm, egrptpALL_GENREST, bVerbose, wi);
3487 auto compressedXGroupNames = gmx::splitString(is->x_compressed_groups);
3488 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3489 SimulationAtomGroupType::CompressedPositionOutput,
3490 restnm, egrptpONE, bVerbose, wi);
3491 auto orirefFitGroupNames = gmx::splitString(is->orirefitgrp);
3492 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3493 SimulationAtomGroupType::OrientationRestraintsFit,
3494 restnm, egrptpALL_GENREST, bVerbose, wi);
3496 /* QMMM input processing */
3497 auto qmGroupNames = gmx::splitString(is->QMMM);
3498 auto qmMethods = gmx::splitString(is->QMmethod);
3499 auto qmBasisSets = gmx::splitString(is->QMbasis);
3500 if (ir->eI != eiMimic)
3502 if (qmMethods.size() != qmGroupNames.size() ||
3503 qmBasisSets.size() != qmGroupNames.size())
3505 gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets"
3506 " and %zu methods\n", qmGroupNames.size(),
3507 qmBasisSets.size(), qmMethods.size());
3509 /* group rest, if any, is always MM! */
3510 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3511 SimulationAtomGroupType::QuantumMechanics,
3512 restnm, egrptpALL_GENREST, bVerbose, wi);
3513 nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
3514 ir->opts.ngQM = qmGroupNames.size();
3515 snew(ir->opts.QMmethod, nr);
3516 snew(ir->opts.QMbasis, nr);
3517 for (i = 0; i < nr; i++)
3519 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3520 * converted to the corresponding enum in names.c
3522 ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(),
3525 ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(),
3530 auto qmMultiplicities = gmx::splitString(is->QMmult);
3531 auto qmCharges = gmx::splitString(is->QMcharge);
3532 auto qmbSH = gmx::splitString(is->bSH);
3533 snew(ir->opts.QMmult, nr);
3534 snew(ir->opts.QMcharge, nr);
3535 snew(ir->opts.bSH, nr);
3536 convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult);
3537 convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge);
3538 convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH);
3540 auto CASelectrons = gmx::splitString(is->CASelectrons);
3541 auto CASorbitals = gmx::splitString(is->CASorbitals);
3542 snew(ir->opts.CASelectrons, nr);
3543 snew(ir->opts.CASorbitals, nr);
3544 convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons);
3545 convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals);
3547 auto SAon = gmx::splitString(is->SAon);
3548 auto SAoff = gmx::splitString(is->SAoff);
3549 auto SAsteps = gmx::splitString(is->SAsteps);
3550 snew(ir->opts.SAon, nr);
3551 snew(ir->opts.SAoff, nr);
3552 snew(ir->opts.SAsteps, nr);
3553 convertInts(wi, SAon, "SAon", ir->opts.SAon);
3554 convertInts(wi, SAoff, "SAoff", ir->opts.SAoff);
3555 convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps);
3560 if (qmGroupNames.size() > 1)
3562 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3564 /* group rest, if any, is always MM! */
3565 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3566 SimulationAtomGroupType::QuantumMechanics,
3567 restnm, egrptpALL_GENREST, bVerbose, wi);
3569 ir->opts.ngQM = qmGroupNames.size();
3572 /* end of QMMM input */
3576 for (auto group : gmx::keysOf(groups->groups))
3578 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3579 for (const auto &entry : groups->groups[group])
3581 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3583 fprintf(stderr, "\n");
3587 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3588 snew(ir->opts.egp_flags, nr*nr);
3590 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3591 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3593 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3595 if (bExcl && EEL_FULL(ir->coulombtype))
3597 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3600 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3601 if (bTable && !(ir->vdwtype == evdwUSER) &&
3602 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3603 !(ir->coulombtype == eelPMEUSERSWITCH))
3605 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3608 /* final check before going out of scope if simulated tempering variables
3609 * need to be set to default values.
3611 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3613 ir->expandedvals->nstexpanded = 2*static_cast<int>(ir->opts.tau_t[0]/ir->delta_t);
3614 warning(wi, gmx::formatString("the value for nstexpanded was not specified for "
3615 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3616 "by default, but it is recommended to set it to an explicit value!",
3617 ir->expandedvals->nstexpanded));
3619 for (i = 0; (i < defaultIndexGroups->nr); i++)
3624 done_blocka(defaultIndexGroups);
3625 sfree(defaultIndexGroups);
3631 static void check_disre(const gmx_mtop_t *mtop)
3633 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3635 const gmx_ffparams_t &ffparams = mtop->ffparams;
3638 for (int i = 0; i < ffparams.numTypes(); i++)
3640 int ftype = ffparams.functype[i];
3641 if (ftype == F_DISRES)
3643 int label = ffparams.iparams[i].disres.label;
3644 if (label == old_label)
3646 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3654 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3655 "probably the parameters for multiple pairs in one restraint "
3656 "are not identical\n", ndouble);
3661 static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3666 gmx_mtop_ilistloop_t iloop;
3675 for (d = 0; d < DIM; d++)
3677 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3679 /* Check for freeze groups */
3680 for (g = 0; g < ir->opts.ngfrz; g++)
3682 for (d = 0; d < DIM; d++)
3684 if (ir->opts.nFreeze[g][d] != 0)
3692 /* Check for position restraints */
3693 iloop = gmx_mtop_ilistloop_init(sys);
3694 while (const InteractionLists *ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
3697 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3699 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
3701 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
3702 for (d = 0; d < DIM; d++)
3704 if (pr->posres.fcA[d] != 0)
3710 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
3712 /* Check for flat-bottom posres */
3713 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
3714 if (pr->fbposres.k != 0)
3716 switch (pr->fbposres.geom)
3718 case efbposresSPHERE:
3719 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3721 case efbposresCYLINDERX:
3722 AbsRef[YY] = AbsRef[ZZ] = 1;
3724 case efbposresCYLINDERY:
3725 AbsRef[XX] = AbsRef[ZZ] = 1;
3727 case efbposresCYLINDER:
3728 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3729 case efbposresCYLINDERZ:
3730 AbsRef[XX] = AbsRef[YY] = 1;
3732 case efbposresX: /* d=XX */
3733 case efbposresY: /* d=YY */
3734 case efbposresZ: /* d=ZZ */
3735 d = pr->fbposres.geom - efbposresX;
3739 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3740 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3748 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3752 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3753 bool *bC6ParametersWorkWithGeometricRules,
3754 bool *bC6ParametersWorkWithLBRules,
3755 bool *bLBRulesPossible)
3757 int ntypes, tpi, tpj;
3760 double c6i, c6j, c12i, c12j;
3761 double c6, c6_geometric, c6_LB;
3762 double sigmai, sigmaj, epsi, epsj;
3763 bool bCanDoLBRules, bCanDoGeometricRules;
3766 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3767 * force-field floating point parameters.
3770 ptr = getenv("GMX_LJCOMB_TOL");
3774 double gmx_unused canary;
3776 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3778 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
3783 *bC6ParametersWorkWithLBRules = TRUE;
3784 *bC6ParametersWorkWithGeometricRules = TRUE;
3785 bCanDoLBRules = TRUE;
3786 ntypes = mtop->ffparams.atnr;
3787 snew(typecount, ntypes);
3788 gmx_mtop_count_atomtypes(mtop, state, typecount);
3789 *bLBRulesPossible = TRUE;
3790 for (tpi = 0; tpi < ntypes; ++tpi)
3792 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3793 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3794 for (tpj = tpi; tpj < ntypes; ++tpj)
3796 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3797 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3798 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3799 c6_geometric = std::sqrt(c6i * c6j);
3800 if (!gmx_numzero(c6_geometric))
3802 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3804 sigmai = gmx::sixthroot(c12i / c6i);
3805 sigmaj = gmx::sixthroot(c12j / c6j);
3806 epsi = c6i * c6i /(4.0 * c12i);
3807 epsj = c6j * c6j /(4.0 * c12j);
3808 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
3812 *bLBRulesPossible = FALSE;
3813 c6_LB = c6_geometric;
3815 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3820 *bC6ParametersWorkWithLBRules = FALSE;
3823 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3825 if (!bCanDoGeometricRules)
3827 *bC6ParametersWorkWithGeometricRules = FALSE;
3835 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3838 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3840 check_combination_rule_differences(mtop, 0,
3841 &bC6ParametersWorkWithGeometricRules,
3842 &bC6ParametersWorkWithLBRules,
3844 if (ir->ljpme_combination_rule == eljpmeLB)
3846 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
3848 warning(wi, "You are using arithmetic-geometric combination rules "
3849 "in LJ-PME, but your non-bonded C6 parameters do not "
3850 "follow these rules.");
3855 if (!bC6ParametersWorkWithGeometricRules)
3857 if (ir->eDispCorr != edispcNO)
3859 warning_note(wi, "You are using geometric combination rules in "
3860 "LJ-PME, but your non-bonded C6 parameters do "
3861 "not follow these rules. "
3862 "This will introduce very small errors in the forces and energies in "
3863 "your simulations. Dispersion correction will correct total energy "
3864 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3868 warning_note(wi, "You are using geometric combination rules in "
3869 "LJ-PME, but your non-bonded C6 parameters do "
3870 "not follow these rules. "
3871 "This will introduce very small errors in the forces and energies in "
3872 "your simulations. If your system is homogeneous, consider using dispersion correction "
3873 "for the total energy and pressure.");
3879 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3882 char err_buf[STRLEN];
3887 gmx_mtop_atomloop_block_t aloopb;
3889 char warn_buf[STRLEN];
3891 set_warning_line(wi, mdparin, -1);
3893 if (ir->cutoff_scheme == ecutsVERLET &&
3894 ir->verletbuf_tol > 0 &&
3896 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
3897 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
3899 /* Check if a too small Verlet buffer might potentially
3900 * cause more drift than the thermostat can couple off.
3902 /* Temperature error fraction for warning and suggestion */
3903 const real T_error_warn = 0.002;
3904 const real T_error_suggest = 0.001;
3905 /* For safety: 2 DOF per atom (typical with constraints) */
3906 const real nrdf_at = 2;
3907 real T, tau, max_T_error;
3912 for (i = 0; i < ir->opts.ngtc; i++)
3914 T = std::max(T, ir->opts.ref_t[i]);
3915 tau = std::max(tau, ir->opts.tau_t[i]);
3919 /* This is a worst case estimate of the temperature error,
3920 * assuming perfect buffer estimation and no cancelation
3921 * of errors. The factor 0.5 is because energy distributes
3922 * equally over Ekin and Epot.
3924 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
3925 if (max_T_error > T_error_warn)
3927 sprintf(warn_buf, "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to %.0e or decrease tau_t.",
3928 ir->verletbuf_tol, T, tau,
3930 100*T_error_suggest,
3931 ir->verletbuf_tol*T_error_suggest/max_T_error);
3932 warning(wi, warn_buf);
3937 if (ETC_ANDERSEN(ir->etc))
3941 for (i = 0; i < ir->opts.ngtc; i++)
3943 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
3944 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
3945 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
3946 i, ir->opts.tau_t[i]);
3947 CHECK(ir->opts.tau_t[i] < 0);
3950 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
3952 for (i = 0; i < ir->opts.ngtc; i++)
3954 int nsteps = gmx::roundToInt(ir->opts.tau_t[i]/ir->delta_t);
3955 sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
3956 CHECK(nsteps % ir->nstcomm != 0);
3961 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3962 ir->comm_mode == ecmNO &&
3963 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
3964 !ETC_ANDERSEN(ir->etc))
3966 warning(wi, "You are not using center of mass motion removal (mdp option comm-mode), numerical rounding errors can lead to build up of kinetic energy of the center of mass");
3969 if (ir->epc == epcPARRINELLORAHMAN &&
3970 ir->etc == etcNOSEHOOVER)
3973 for (int g = 0; g < ir->opts.ngtc; g++)
3975 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
3977 if (ir->tau_p < 1.9*tau_t_max)
3979 std::string message =
3980 gmx::formatString("With %s T-coupling and %s p-coupling, "
3981 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
3982 etcoupl_names[ir->etc],
3983 epcoupl_names[ir->epc],
3985 "tau-t", tau_t_max);
3986 warning(wi, message.c_str());
3990 /* Check for pressure coupling with absolute position restraints */
3991 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3993 absolute_reference(ir, sys, TRUE, AbsRef);
3995 for (m = 0; m < DIM; m++)
3997 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3999 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4007 aloopb = gmx_mtop_atomloop_block_init(sys);
4009 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4011 if (atom->q != 0 || atom->qB != 0)
4019 if (EEL_FULL(ir->coulombtype))
4022 "You are using full electrostatics treatment %s for a system without charges.\n"
4023 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4024 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4025 warning(wi, err_buf);
4030 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4033 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4034 "You might want to consider using %s electrostatics.\n",
4036 warning_note(wi, err_buf);
4040 /* Check if combination rules used in LJ-PME are the same as in the force field */
4041 if (EVDW_PME(ir->vdwtype))
4043 check_combination_rules(ir, sys, wi);
4046 /* Generalized reaction field */
4047 if (ir->coulombtype == eelGRF_NOTUSED)
4049 warning_error(wi, "Generalized reaction-field electrostatics is no longer supported. "
4050 "You can use normal reaction-field instead and compute the reaction-field constant by hand.");
4054 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4056 for (m = 0; (m < DIM); m++)
4058 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4067 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4068 for (const AtomProxy atomP : AtomRange(*sys))
4070 const t_atom &local = atomP.atom();
4071 int i = atomP.globalAtomNumber();
4072 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4075 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4077 for (m = 0; (m < DIM); m++)
4079 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4083 for (m = 0; (m < DIM); m++)
4085 if (fabs(acc[m]) > 1e-6)
4087 const char *dim[DIM] = { "X", "Y", "Z" };
4089 "Net Acceleration in %s direction, will %s be corrected\n",
4090 dim[m], ir->nstcomm != 0 ? "" : "not");
4091 if (ir->nstcomm != 0 && m < ndof_com(ir))
4094 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4096 ir->opts.acc[i][m] -= acc[m];
4104 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4105 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4107 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4115 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4117 if (ir->pull->coord[i].group[0] == 0 ||
4118 ir->pull->coord[i].group[1] == 0)
4120 absolute_reference(ir, sys, FALSE, AbsRef);
4121 for (m = 0; m < DIM; m++)
4123 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4125 warning(wi, "You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
4133 for (i = 0; i < 3; i++)
4135 for (m = 0; m <= i; m++)
4137 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4138 ir->deform[i][m] != 0)
4140 for (c = 0; c < ir->pull->ncoord; c++)
4142 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4143 ir->pull->coord[c].vec[m] != 0)
4145 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4156 void double_check(t_inputrec *ir, matrix box,
4157 bool bHasNormalConstraints,
4158 bool bHasAnyConstraints,
4162 char warn_buf[STRLEN];
4165 ptr = check_box(ir->ePBC, box);
4168 warning_error(wi, ptr);
4171 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4173 if (ir->shake_tol <= 0.0)
4175 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4177 warning_error(wi, warn_buf);
4181 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4183 /* If we have Lincs constraints: */
4184 if (ir->eI == eiMD && ir->etc == etcNO &&
4185 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4187 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4188 warning_note(wi, warn_buf);
4191 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4193 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4194 warning_note(wi, warn_buf);
4196 if (ir->epc == epcMTTK)
4198 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4202 if (bHasAnyConstraints && ir->epc == epcMTTK)
4204 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4207 if (ir->LincsWarnAngle > 90.0)
4209 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4210 warning(wi, warn_buf);
4211 ir->LincsWarnAngle = 90.0;
4214 if (ir->ePBC != epbcNONE)
4216 if (ir->nstlist == 0)
4218 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4220 if (ir->ns_type == ensGRID)
4222 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4224 sprintf(warn_buf, "ERROR: The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease rlist.\n");
4225 warning_error(wi, warn_buf);
4230 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4231 if (2*ir->rlist >= min_size)
4233 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4234 warning_error(wi, warn_buf);
4237 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4244 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4248 real rvdw1, rvdw2, rcoul1, rcoul2;
4249 char warn_buf[STRLEN];
4251 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4255 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4260 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4266 if (rvdw1 + rvdw2 > ir->rlist ||
4267 rcoul1 + rcoul2 > ir->rlist)
4270 "The sum of the two largest charge group radii (%f) "
4271 "is larger than rlist (%f)\n",
4272 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4273 warning(wi, warn_buf);
4277 /* Here we do not use the zero at cut-off macro,
4278 * since user defined interactions might purposely
4279 * not be zero at the cut-off.
4281 if (ir_vdw_is_zero_at_cutoff(ir) &&
4282 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4284 sprintf(warn_buf, "The sum of the two largest charge group "
4285 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4286 "With exact cut-offs, better performance can be "
4287 "obtained with cutoff-scheme = %s, because it "
4288 "does not use charge groups at all.",
4290 ir->rlist, ir->rvdw,
4291 ecutscheme_names[ecutsVERLET]);
4294 warning(wi, warn_buf);
4298 warning_note(wi, warn_buf);
4301 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4302 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4304 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4305 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4307 ir->rlist, ir->rcoulomb,
4308 ecutscheme_names[ecutsVERLET]);
4311 warning(wi, warn_buf);
4315 warning_note(wi, warn_buf);