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/topology/block.h"
67 #include "gromacs/topology/ifunc.h"
68 #include "gromacs/topology/index.h"
69 #include "gromacs/topology/mtop_util.h"
70 #include "gromacs/topology/symtab.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/utility/cstringutil.h"
73 #include "gromacs/utility/exceptions.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/filestream.h"
76 #include "gromacs/utility/gmxassert.h"
77 #include "gromacs/utility/ikeyvaluetreeerror.h"
78 #include "gromacs/utility/keyvaluetree.h"
79 #include "gromacs/utility/keyvaluetreebuilder.h"
80 #include "gromacs/utility/keyvaluetreemdpwriter.h"
81 #include "gromacs/utility/keyvaluetreetransform.h"
82 #include "gromacs/utility/smalloc.h"
83 #include "gromacs/utility/strconvert.h"
84 #include "gromacs/utility/stringcompare.h"
85 #include "gromacs/utility/stringutil.h"
86 #include "gromacs/utility/textwriter.h"
91 /* Resource parameters
92 * Do not change any of these until you read the instruction
93 * in readinp.h. Some cpp's do not take spaces after the backslash
94 * (like the c-shell), which will give you a very weird compiler
98 typedef struct t_inputrec_strings
100 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
101 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
102 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
103 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
104 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
106 char fep_lambda[efptNR][STRLEN];
107 char lambda_weights[STRLEN];
110 char anneal[STRLEN], anneal_npoints[STRLEN],
111 anneal_time[STRLEN], anneal_temp[STRLEN];
112 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
113 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
114 SAoff[STRLEN], SAsteps[STRLEN];
116 } gmx_inputrec_strings;
118 static gmx_inputrec_strings *is = nullptr;
120 void init_inputrec_strings()
124 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.");
129 void done_inputrec_strings()
137 egrptpALL, /* All particles have to be a member of a group. */
138 egrptpALL_GENREST, /* A rest group with name is generated for particles *
139 * that are not part of any group. */
140 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
141 * for the rest group. */
142 egrptpONE /* Merge all selected groups into one group, *
143 * make a rest group for the remaining particles. */
146 static const char *constraints[eshNR+1] = {
147 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
150 static const char *couple_lam[ecouplamNR+1] = {
151 "vdw-q", "vdw", "q", "none", nullptr
154 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
159 for (i = 0; i < ntemps; i++)
161 /* simple linear scaling -- allows more control */
162 if (simtemp->eSimTempScale == esimtempLINEAR)
164 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
166 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
168 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
170 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
172 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
177 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
178 gmx_fatal(FARGS, "%s", errorstr);
185 static void _low_check(bool b, const char *s, warninp_t wi)
189 warning_error(wi, s);
193 static void check_nst(const char *desc_nst, int nst,
194 const char *desc_p, int *p,
199 if (*p > 0 && *p % nst != 0)
201 /* Round up to the next multiple of nst */
202 *p = ((*p)/nst + 1)*nst;
203 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
204 desc_p, desc_nst, desc_p, *p);
209 static bool ir_NVE(const t_inputrec *ir)
211 return (EI_MD(ir->eI) && ir->etc == etcNO);
214 static int lcd(int n1, int n2)
219 for (i = 2; (i <= n1 && i <= n2); i++)
221 if (n1 % i == 0 && n2 % i == 0)
230 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
232 if (*eintmod == eintmodPOTSHIFT_VERLET)
234 if (ir->cutoff_scheme == ecutsVERLET)
236 *eintmod = eintmodPOTSHIFT;
240 *eintmod = eintmodNONE;
245 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
247 /* Check internal consistency.
248 * NOTE: index groups are not set here yet, don't check things
249 * like temperature coupling group options here, but in triple_check
252 /* Strange macro: first one fills the err_buf, and then one can check
253 * the condition, which will print the message and increase the error
256 #define CHECK(b) _low_check(b, err_buf, wi)
257 char err_buf[256], warn_buf[STRLEN];
260 t_lambda *fep = ir->fepvals;
261 t_expanded *expand = ir->expandedvals;
263 set_warning_line(wi, mdparin, -1);
265 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
267 sprintf(warn_buf, "%s electrostatics is no longer supported",
268 eel_names[eelRF_NEC_UNSUPPORTED]);
269 warning_error(wi, warn_buf);
272 /* BASIC CUT-OFF STUFF */
273 if (ir->rcoulomb < 0)
275 warning_error(wi, "rcoulomb should be >= 0");
279 warning_error(wi, "rvdw should be >= 0");
282 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
284 warning_error(wi, "rlist should be >= 0");
286 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.)");
287 CHECK(ir->nstlist < 0);
289 process_interaction_modifier(ir, &ir->coulomb_modifier);
290 process_interaction_modifier(ir, &ir->vdw_modifier);
292 if (ir->cutoff_scheme == ecutsGROUP)
295 "The group cutoff scheme has been removed since GROMACS 2020. "
296 "Please use the Verlet cutoff scheme.");
298 if (ir->cutoff_scheme == ecutsVERLET)
302 /* Normal Verlet type neighbor-list, currently only limited feature support */
303 if (inputrec2nboundeddim(ir) < 3)
305 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
308 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
309 if (ir->rcoulomb != ir->rvdw)
311 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
312 // for PME load balancing, we can support this exception.
313 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
314 ir->vdwtype == evdwCUT &&
315 ir->rcoulomb > ir->rvdw);
316 if (!bUsesPmeTwinRangeKernel)
318 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
322 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
324 if (ir->vdw_modifier == eintmodNONE ||
325 ir->vdw_modifier == eintmodPOTSHIFT)
327 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
329 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]);
330 warning_note(wi, warn_buf);
332 ir->vdwtype = evdwCUT;
336 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
337 warning_error(wi, warn_buf);
341 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
343 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
345 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
346 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
348 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
350 if (!(ir->coulomb_modifier == eintmodNONE ||
351 ir->coulomb_modifier == eintmodPOTSHIFT))
353 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
354 warning_error(wi, warn_buf);
357 if (EEL_USER(ir->coulombtype))
359 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme", eel_names[ir->coulombtype]);
360 warning_error(wi, warn_buf);
363 if (ir->nstlist <= 0)
365 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
368 if (ir->nstlist < 10)
370 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.");
373 rc_max = std::max(ir->rvdw, ir->rcoulomb);
375 if (ir->verletbuf_tol <= 0)
377 if (ir->verletbuf_tol == 0)
379 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
382 if (ir->rlist < rc_max)
384 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
387 if (ir->rlist == rc_max && ir->nstlist > 1)
389 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.");
394 if (ir->rlist > rc_max)
396 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.");
399 if (ir->nstlist == 1)
401 /* No buffer required */
406 if (EI_DYNAMICS(ir->eI))
408 if (inputrec2nboundeddim(ir) < 3)
410 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.");
412 /* Set rlist temporarily so we can continue processing */
417 /* Set the buffer to 5% of the cut-off */
418 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
424 /* GENERAL INTEGRATOR STUFF */
427 if (ir->etc != etcNO)
429 if (EI_RANDOM(ir->eI))
431 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]);
435 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
437 warning_note(wi, warn_buf);
441 if (ir->eI == eiVVAK)
443 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]);
444 warning_note(wi, warn_buf);
446 if (!EI_DYNAMICS(ir->eI))
448 if (ir->epc != epcNO)
450 sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
451 warning_note(wi, warn_buf);
455 if (EI_DYNAMICS(ir->eI))
457 if (ir->nstcalcenergy < 0)
459 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
460 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
462 /* nstcalcenergy larger than nstener does not make sense.
463 * We ideally want nstcalcenergy=nstener.
467 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
471 ir->nstcalcenergy = ir->nstenergy;
475 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
476 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
477 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
480 const char *nsten = "nstenergy";
481 const char *nstdh = "nstdhdl";
482 const char *min_name = nsten;
483 int min_nst = ir->nstenergy;
485 /* find the smallest of ( nstenergy, nstdhdl ) */
486 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
487 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
489 min_nst = ir->fepvals->nstdhdl;
492 /* If the user sets nstenergy small, we should respect that */
494 "Setting nstcalcenergy (%d) equal to %s (%d)",
495 ir->nstcalcenergy, min_name, min_nst);
496 warning_note(wi, warn_buf);
497 ir->nstcalcenergy = min_nst;
500 if (ir->epc != epcNO)
502 if (ir->nstpcouple < 0)
504 ir->nstpcouple = ir_optimal_nstpcouple(ir);
508 if (ir->nstcalcenergy > 0)
510 if (ir->efep != efepNO)
512 /* nstdhdl should be a multiple of nstcalcenergy */
513 check_nst("nstcalcenergy", ir->nstcalcenergy,
514 "nstdhdl", &ir->fepvals->nstdhdl, wi);
518 /* nstexpanded should be a multiple of nstcalcenergy */
519 check_nst("nstcalcenergy", ir->nstcalcenergy,
520 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
522 /* for storing exact averages nstenergy should be
523 * a multiple of nstcalcenergy
525 check_nst("nstcalcenergy", ir->nstcalcenergy,
526 "nstenergy", &ir->nstenergy, wi);
530 if (ir->nsteps == 0 && !ir->bContinuation)
532 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
536 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
537 ir->bContinuation && ir->ld_seed != -1)
539 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)");
545 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
546 CHECK(ir->ePBC != epbcXYZ);
547 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
548 CHECK(ir->ns_type != ensGRID);
549 sprintf(err_buf, "with TPI nstlist should be larger than zero");
550 CHECK(ir->nstlist <= 0);
551 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
552 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
553 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
554 CHECK(ir->cutoff_scheme == ecutsVERLET);
558 if ( (opts->nshake > 0) && (opts->bMorse) )
561 "Using morse bond-potentials while constraining bonds is useless");
562 warning(wi, warn_buf);
565 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
566 ir->bContinuation && ir->ld_seed != -1)
568 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)");
570 /* verify simulated tempering options */
574 bool bAllTempZero = TRUE;
575 for (i = 0; i < fep->n_lambda; i++)
577 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]);
578 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
579 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
581 bAllTempZero = FALSE;
584 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
585 CHECK(bAllTempZero == TRUE);
587 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
588 CHECK(ir->eI != eiVV);
590 /* check compatability of the temperature coupling with simulated tempering */
592 if (ir->etc == etcNOSEHOOVER)
594 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
595 warning_note(wi, warn_buf);
598 /* check that the temperatures make sense */
600 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);
601 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
603 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
604 CHECK(ir->simtempvals->simtemp_high <= 0);
606 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
607 CHECK(ir->simtempvals->simtemp_low <= 0);
610 /* verify free energy options */
612 if (ir->efep != efepNO)
615 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
617 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
619 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
620 static_cast<int>(fep->sc_r_power));
621 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
623 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
624 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
626 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
627 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
629 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
630 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
632 sprintf(err_buf, "Free-energy not implemented for Ewald");
633 CHECK(ir->coulombtype == eelEWALD);
635 /* check validty of lambda inputs */
636 if (fep->n_lambda == 0)
638 /* Clear output in case of no states:*/
639 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
640 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
644 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
645 CHECK((fep->init_fep_state >= fep->n_lambda));
648 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
649 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
651 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",
652 fep->init_lambda, fep->init_fep_state);
653 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
657 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
661 for (i = 0; i < efptNR; i++)
663 if (fep->separate_dvdl[i])
668 if (n_lambda_terms > 1)
670 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.");
671 warning(wi, warn_buf);
674 if (n_lambda_terms < 2 && fep->n_lambda > 0)
677 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
681 for (j = 0; j < efptNR; j++)
683 for (i = 0; i < fep->n_lambda; i++)
685 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]);
686 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
690 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
692 for (i = 0; i < fep->n_lambda; i++)
694 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],
695 fep->all_lambda[efptCOUL][i]);
696 CHECK((fep->sc_alpha > 0) &&
697 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
698 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
699 ((fep->all_lambda[efptVDW][i] > 0.0) &&
700 (fep->all_lambda[efptVDW][i] < 1.0))));
704 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
706 real sigma, lambda, r_sc;
709 /* Maximum estimate for A and B charges equal with lambda power 1 */
711 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);
712 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.",
714 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
715 warning_note(wi, warn_buf);
718 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
719 be treated differently, but that's the next step */
721 for (i = 0; i < efptNR; i++)
723 for (j = 0; j < fep->n_lambda; j++)
725 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
726 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
731 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
735 /* checking equilibration of weights inputs for validity */
737 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
738 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
739 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
741 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
742 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
743 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
745 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
746 expand->equil_steps, elmceq_names[elmceqSTEPS]);
747 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
749 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
750 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
751 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
753 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
754 expand->equil_ratio, elmceq_names[elmceqRATIO]);
755 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
757 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
758 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
759 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
761 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
762 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
763 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
765 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
766 expand->equil_steps, elmceq_names[elmceqSTEPS]);
767 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
769 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
770 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
771 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
773 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
774 expand->equil_ratio, elmceq_names[elmceqRATIO]);
775 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
777 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
778 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
779 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
781 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
782 CHECK((expand->lmc_repeats <= 0));
783 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
784 CHECK((expand->minvarmin <= 0));
785 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
786 CHECK((expand->c_range < 0));
787 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
788 fep->init_fep_state, expand->lmc_forced_nstart);
789 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
790 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
791 CHECK((expand->lmc_forced_nstart < 0));
792 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
793 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
795 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
796 CHECK((expand->init_wl_delta < 0));
797 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
798 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
799 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
800 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
802 /* if there is no temperature control, we need to specify an MC temperature */
803 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO) && (expand->mc_temp <= 0.0))
805 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set to a positive number");
806 warning_error(wi, err_buf);
808 if (expand->nstTij > 0)
810 sprintf(err_buf, "nstlog must be non-zero");
811 CHECK(ir->nstlog == 0);
812 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
813 expand->nstTij, ir->nstlog);
814 CHECK((expand->nstTij % ir->nstlog) != 0);
819 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
820 CHECK(ir->nwall && ir->ePBC != epbcXY);
823 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
825 if (ir->ePBC == epbcNONE)
827 if (ir->epc != epcNO)
829 warning(wi, "Turning off pressure coupling for vacuum system");
835 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
836 epbc_names[ir->ePBC]);
837 CHECK(ir->epc != epcNO);
839 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
840 CHECK(EEL_FULL(ir->coulombtype));
842 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
843 epbc_names[ir->ePBC]);
844 CHECK(ir->eDispCorr != edispcNO);
847 if (ir->rlist == 0.0)
849 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
850 "with coulombtype = %s or coulombtype = %s\n"
851 "without periodic boundary conditions (pbc = %s) and\n"
852 "rcoulomb and rvdw set to zero",
853 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
854 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
855 (ir->ePBC != epbcNONE) ||
856 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
860 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
865 if (ir->nstcomm == 0)
867 // TODO Change this behaviour. There should be exactly one way
868 // to turn off an algorithm.
869 ir->comm_mode = ecmNO;
871 if (ir->comm_mode != ecmNO)
875 // TODO Such input was once valid. Now that we've been
876 // helpful for a few years, we should reject such input,
877 // lest we have to support every historical decision
879 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");
880 ir->nstcomm = abs(ir->nstcomm);
883 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
885 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
886 ir->nstcomm = ir->nstcalcenergy;
889 if (ir->comm_mode == ecmANGULAR)
891 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
892 CHECK(ir->bPeriodicMols);
893 if (ir->ePBC != epbcNONE)
895 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.");
900 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
902 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]);
903 warning_note(wi, warn_buf);
906 /* TEMPERATURE COUPLING */
907 if (ir->etc == etcYES)
909 ir->etc = etcBERENDSEN;
910 warning_note(wi, "Old option for temperature coupling given: "
911 "changing \"yes\" to \"Berendsen\"\n");
914 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
916 if (ir->opts.nhchainlength < 1)
918 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
919 ir->opts.nhchainlength = 1;
920 warning(wi, warn_buf);
923 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
925 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
926 ir->opts.nhchainlength = 1;
931 ir->opts.nhchainlength = 0;
934 if (ir->eI == eiVVAK)
936 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
938 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
941 if (ETC_ANDERSEN(ir->etc))
943 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
944 CHECK(!(EI_VV(ir->eI)));
946 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
948 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]);
949 warning_note(wi, warn_buf);
952 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]);
953 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
956 if (ir->etc == etcBERENDSEN)
958 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
959 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
960 warning_note(wi, warn_buf);
963 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
964 && ir->epc == epcBERENDSEN)
966 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
967 "true ensemble for the thermostat");
968 warning(wi, warn_buf);
971 /* PRESSURE COUPLING */
972 if (ir->epc == epcISOTROPIC)
974 ir->epc = epcBERENDSEN;
975 warning_note(wi, "Old option for pressure coupling given: "
976 "changing \"Isotropic\" to \"Berendsen\"\n");
979 if (ir->epc != epcNO)
981 dt_pcoupl = ir->nstpcouple*ir->delta_t;
983 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
984 CHECK(ir->tau_p <= 0);
986 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
988 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
989 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
990 warning(wi, warn_buf);
993 sprintf(err_buf, "compressibility must be > 0 when using pressure"
994 " coupling %s\n", EPCOUPLTYPE(ir->epc));
995 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
996 ir->compress[ZZ][ZZ] < 0 ||
997 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
998 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1000 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1003 "You are generating velocities so I am assuming you "
1004 "are equilibrating a system. You are using "
1005 "%s pressure coupling, but this can be "
1006 "unstable for equilibration. If your system crashes, try "
1007 "equilibrating first with Berendsen pressure coupling. If "
1008 "you are not equilibrating the system, you can probably "
1009 "ignore this warning.",
1010 epcoupl_names[ir->epc]);
1011 warning(wi, warn_buf);
1017 if (ir->epc > epcNO)
1019 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1021 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.");
1027 if (ir->epc == epcMTTK)
1029 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1033 /* ELECTROSTATICS */
1034 /* More checks are in triple check (grompp.c) */
1036 if (ir->coulombtype == eelSWITCH)
1038 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1039 "artifacts, advice: use coulombtype = %s",
1040 eel_names[ir->coulombtype],
1041 eel_names[eelRF_ZERO]);
1042 warning(wi, warn_buf);
1045 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1047 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);
1048 warning(wi, warn_buf);
1049 ir->epsilon_rf = ir->epsilon_r;
1050 ir->epsilon_r = 1.0;
1053 if (ir->epsilon_r == 0)
1056 "It is pointless to use long-range electrostatics with infinite relative permittivity."
1057 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1058 CHECK(EEL_FULL(ir->coulombtype));
1061 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1063 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1064 CHECK(ir->epsilon_r < 0);
1067 if (EEL_RF(ir->coulombtype))
1069 /* reaction field (at the cut-off) */
1071 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1073 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1074 eel_names[ir->coulombtype]);
1075 warning(wi, warn_buf);
1076 ir->epsilon_rf = 0.0;
1079 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1080 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1081 (ir->epsilon_r == 0));
1082 if (ir->epsilon_rf == ir->epsilon_r)
1084 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1085 eel_names[ir->coulombtype]);
1086 warning(wi, warn_buf);
1089 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1090 * means the interaction is zero outside rcoulomb, but it helps to
1091 * provide accurate energy conservation.
1093 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1095 if (ir_coulomb_switched(ir))
1098 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1099 eel_names[ir->coulombtype]);
1100 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1104 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1107 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1108 CHECK( ir->coulomb_modifier != eintmodNONE);
1110 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1113 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1114 CHECK( ir->vdw_modifier != eintmodNONE);
1117 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1118 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1121 "The switch/shift interaction settings are just for compatibility; you will get better "
1122 "performance from applying potential modifiers to your interactions!\n");
1123 warning_note(wi, warn_buf);
1126 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1128 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1130 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1131 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.",
1132 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1133 warning(wi, warn_buf);
1137 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1139 if (ir->rvdw_switch == 0)
1141 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.");
1142 warning(wi, warn_buf);
1146 if (EEL_FULL(ir->coulombtype))
1148 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1149 ir->coulombtype == eelPMEUSERSWITCH)
1151 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1152 eel_names[ir->coulombtype]);
1153 CHECK(ir->rcoulomb > ir->rlist);
1157 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1159 // TODO: Move these checks into the ewald module with the options class
1161 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1163 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1165 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
1166 warning_error(wi, warn_buf);
1170 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1172 if (ir->ewald_geometry == eewg3D)
1174 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1175 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1176 warning(wi, warn_buf);
1178 /* This check avoids extra pbc coding for exclusion corrections */
1179 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1180 CHECK(ir->wall_ewald_zfac < 2);
1182 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1183 EEL_FULL(ir->coulombtype))
1185 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1186 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1187 warning(wi, warn_buf);
1189 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1191 if (ir->cutoff_scheme == ecutsVERLET)
1193 sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
1194 eel_names[ir->coulombtype]);
1195 warning(wi, warn_buf);
1199 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",
1200 eel_names[ir->coulombtype]);
1201 warning_note(wi, warn_buf);
1205 if (ir_vdw_switched(ir))
1207 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1208 CHECK(ir->rvdw_switch >= ir->rvdw);
1210 if (ir->rvdw_switch < 0.5*ir->rvdw)
1212 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.",
1213 ir->rvdw_switch, ir->rvdw);
1214 warning_note(wi, warn_buf);
1218 if (ir->vdwtype == evdwPME)
1220 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1222 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1223 evdw_names[ir->vdwtype],
1224 eintmod_names[eintmodPOTSHIFT],
1225 eintmod_names[eintmodNONE]);
1226 warning_error(wi, err_buf);
1230 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1232 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.");
1235 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1238 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1241 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1243 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1246 /* IMPLICIT SOLVENT */
1247 if (ir->coulombtype == eelGB_NOTUSED)
1249 sprintf(warn_buf, "Invalid option %s for coulombtype",
1250 eel_names[ir->coulombtype]);
1251 warning_error(wi, warn_buf);
1256 if (ir->cutoff_scheme != ecutsGROUP)
1258 warning_error(wi, "QMMM is currently only supported with cutoff-scheme=group");
1260 if (!EI_DYNAMICS(ir->eI))
1263 sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
1264 warning_error(wi, buf);
1270 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1274 /* interpret a number of doubles from a string and put them in an array,
1275 after allocating space for them.
1276 str = the input string
1277 n = the (pre-allocated) number of doubles read
1278 r = the output array of doubles. */
1279 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1281 auto values = gmx::splitString(str);
1285 for (int i = 0; i < *n; i++)
1289 (*r)[i] = gmx::fromString<real>(values[i]);
1291 catch (gmx::GromacsException &)
1293 warning_error(wi, "Invalid value " + values[i] + " in string in mdp file. Expected a real number.");
1299 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1302 int i, j, max_n_lambda, nweights, nfep[efptNR];
1303 t_lambda *fep = ir->fepvals;
1304 t_expanded *expand = ir->expandedvals;
1305 real **count_fep_lambdas;
1306 bool bOneLambda = TRUE;
1308 snew(count_fep_lambdas, efptNR);
1310 /* FEP input processing */
1311 /* first, identify the number of lambda values for each type.
1312 All that are nonzero must have the same number */
1314 for (i = 0; i < efptNR; i++)
1316 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1319 /* now, determine the number of components. All must be either zero, or equal. */
1322 for (i = 0; i < efptNR; i++)
1324 if (nfep[i] > max_n_lambda)
1326 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1327 must have the same number if its not zero.*/
1332 for (i = 0; i < efptNR; i++)
1336 ir->fepvals->separate_dvdl[i] = FALSE;
1338 else if (nfep[i] == max_n_lambda)
1340 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1341 respect to the temperature currently */
1343 ir->fepvals->separate_dvdl[i] = TRUE;
1348 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1349 nfep[i], efpt_names[i], max_n_lambda);
1352 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1353 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1355 /* the number of lambdas is the number we've read in, which is either zero
1356 or the same for all */
1357 fep->n_lambda = max_n_lambda;
1359 /* allocate space for the array of lambda values */
1360 snew(fep->all_lambda, efptNR);
1361 /* if init_lambda is defined, we need to set lambda */
1362 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1364 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1366 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1367 for (i = 0; i < efptNR; i++)
1369 snew(fep->all_lambda[i], fep->n_lambda);
1370 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1373 for (j = 0; j < fep->n_lambda; j++)
1375 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1377 sfree(count_fep_lambdas[i]);
1380 sfree(count_fep_lambdas);
1382 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1383 bookkeeping -- for now, init_lambda */
1385 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1387 for (i = 0; i < fep->n_lambda; i++)
1389 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1393 /* check to see if only a single component lambda is defined, and soft core is defined.
1394 In this case, turn on coulomb soft core */
1396 if (max_n_lambda == 0)
1402 for (i = 0; i < efptNR; i++)
1404 if ((nfep[i] != 0) && (i != efptFEP))
1410 if ((bOneLambda) && (fep->sc_alpha > 0))
1412 fep->bScCoul = TRUE;
1415 /* Fill in the others with the efptFEP if they are not explicitly
1416 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1417 they are all zero. */
1419 for (i = 0; i < efptNR; i++)
1421 if ((nfep[i] == 0) && (i != efptFEP))
1423 for (j = 0; j < fep->n_lambda; j++)
1425 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1431 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1432 if (fep->sc_r_power == 48)
1434 if (fep->sc_alpha > 0.1)
1436 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1440 /* now read in the weights */
1441 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1444 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1446 else if (nweights != fep->n_lambda)
1448 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1449 nweights, fep->n_lambda);
1451 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1453 expand->nstexpanded = fep->nstdhdl;
1454 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1459 static void do_simtemp_params(t_inputrec *ir)
1462 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1463 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1467 convertYesNos(warninp_t /*wi*/, gmx::ArrayRef<const std::string> inputs, const char * /*name*/, gmx_bool *outputs)
1470 for (const auto &input : inputs)
1472 outputs[i] = gmx::equalCaseInsensitive(input, "Y", 1);
1477 template <typename T> void
1478 convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, T *outputs)
1481 for (const auto &input : inputs)
1485 outputs[i] = gmx::fromStdString<T>(input);
1487 catch (gmx::GromacsException &)
1489 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of integers separated by spaces.",
1491 warning_error(wi, message);
1498 convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, real *outputs)
1501 for (const auto &input : inputs)
1505 outputs[i] = gmx::fromString<real>(input);
1507 catch (gmx::GromacsException &)
1509 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of real numbers separated by spaces.",
1511 warning_error(wi, message);
1518 convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, rvec *outputs)
1521 for (const auto &input : inputs)
1525 outputs[i][d] = gmx::fromString<real>(input);
1527 catch (gmx::GromacsException &)
1529 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of real numbers separated by spaces.",
1531 warning_error(wi, message);
1542 static void do_wall_params(t_inputrec *ir,
1543 char *wall_atomtype, char *wall_density,
1547 opts->wall_atomtype[0] = nullptr;
1548 opts->wall_atomtype[1] = nullptr;
1550 ir->wall_atomtype[0] = -1;
1551 ir->wall_atomtype[1] = -1;
1552 ir->wall_density[0] = 0;
1553 ir->wall_density[1] = 0;
1557 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1558 if (wallAtomTypes.size() != size_t(ir->nwall))
1560 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu",
1561 ir->nwall, wallAtomTypes.size());
1563 for (int i = 0; i < ir->nwall; i++)
1565 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1568 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1570 auto wallDensity = gmx::splitString(wall_density);
1571 if (wallDensity.size() != size_t(ir->nwall))
1573 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall, wallDensity.size());
1575 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1576 for (int i = 0; i < ir->nwall; i++)
1578 if (ir->wall_density[i] <= 0)
1580 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1587 static void add_wall_energrps(SimulationGroups *groups, int nwall, t_symtab *symtab)
1591 AtomGroupIndices *grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1592 for (int i = 0; i < nwall; i++)
1594 groups->groupNames.emplace_back(
1597 gmx::formatString("wall%d", i).c_str()));
1598 grps->emplace_back(groups->groupNames.size() - 1);
1603 static void read_expandedparams(std::vector<t_inpfile> *inp,
1604 t_expanded *expand, warninp_t wi)
1606 /* read expanded ensemble parameters */
1607 printStringNewline(inp, "expanded ensemble variables");
1608 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1609 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1610 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1611 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1612 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1613 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1614 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1615 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1616 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1617 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1618 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1619 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1620 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1621 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1622 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1623 expand->bSymmetrizedTMatrix = (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1624 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1625 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1626 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1627 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1628 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1629 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1630 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1633 /*! \brief Return whether an end state with the given coupling-lambda
1634 * value describes fully-interacting VDW.
1636 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1637 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1639 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1641 return (couple_lambda_value == ecouplamVDW ||
1642 couple_lambda_value == ecouplamVDWQ);
1648 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1651 explicit MdpErrorHandler(warninp_t wi)
1652 : wi_(wi), mapping_(nullptr)
1656 void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
1658 mapping_ = &mapping;
1661 bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context) override
1663 ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
1664 getOptionName(context).c_str()));
1665 std::string message = gmx::formatExceptionMessageToString(*ex);
1666 warning_error(wi_, message.c_str());
1671 std::string getOptionName(const gmx::KeyValueTreePath &context)
1673 if (mapping_ != nullptr)
1675 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1676 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1679 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1684 const gmx::IKeyValueTreeBackMapping *mapping_;
1689 void get_ir(const char *mdparin, const char *mdparout,
1690 gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
1691 WriteMdpHeader writeMdpHeader, warninp_t wi)
1694 double dumdub[2][6];
1696 char warn_buf[STRLEN];
1697 t_lambda *fep = ir->fepvals;
1698 t_expanded *expand = ir->expandedvals;
1700 const char *no_names[] = { "no", nullptr };
1702 init_inputrec_strings();
1703 gmx::TextInputFile stream(mdparin);
1704 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1706 snew(dumstr[0], STRLEN);
1707 snew(dumstr[1], STRLEN);
1709 if (-1 == search_einp(inp, "cutoff-scheme"))
1712 "%s did not specify a value for the .mdp option "
1713 "\"cutoff-scheme\". As of GROMACS 2020, the Verlet scheme "
1714 "is the only cutoff scheme supported. This may affect your "
1715 "simulation if you are using an old mdp file that assumes use "
1716 "of the (removed) group cutoff scheme.", mdparin);
1717 warning_note(wi, warn_buf);
1720 /* ignore the following deprecated commands */
1721 replace_inp_entry(inp, "title", nullptr);
1722 replace_inp_entry(inp, "cpp", nullptr);
1723 replace_inp_entry(inp, "domain-decomposition", nullptr);
1724 replace_inp_entry(inp, "andersen-seed", nullptr);
1725 replace_inp_entry(inp, "dihre", nullptr);
1726 replace_inp_entry(inp, "dihre-fc", nullptr);
1727 replace_inp_entry(inp, "dihre-tau", nullptr);
1728 replace_inp_entry(inp, "nstdihreout", nullptr);
1729 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1730 replace_inp_entry(inp, "optimize-fft", nullptr);
1731 replace_inp_entry(inp, "adress_type", nullptr);
1732 replace_inp_entry(inp, "adress_const_wf", nullptr);
1733 replace_inp_entry(inp, "adress_ex_width", nullptr);
1734 replace_inp_entry(inp, "adress_hy_width", nullptr);
1735 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1736 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1737 replace_inp_entry(inp, "adress_site", nullptr);
1738 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1739 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1740 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1741 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1742 replace_inp_entry(inp, "rlistlong", nullptr);
1743 replace_inp_entry(inp, "nstcalclr", nullptr);
1744 replace_inp_entry(inp, "pull-print-com2", nullptr);
1745 replace_inp_entry(inp, "gb-algorithm", nullptr);
1746 replace_inp_entry(inp, "nstgbradii", nullptr);
1747 replace_inp_entry(inp, "rgbradii", nullptr);
1748 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1749 replace_inp_entry(inp, "gb-saltconc", nullptr);
1750 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1751 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1752 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1753 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1754 replace_inp_entry(inp, "sa-algorithm", nullptr);
1755 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1757 /* replace the following commands with the clearer new versions*/
1758 replace_inp_entry(inp, "unconstrained-start", "continuation");
1759 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1760 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1761 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1762 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1763 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1764 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1766 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1767 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1768 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1769 setStringEntry(&inp, "include", opts->include, nullptr);
1770 printStringNoNewline(&inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1771 setStringEntry(&inp, "define", opts->define, nullptr);
1773 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1774 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1775 printStringNoNewline(&inp, "Start time and timestep in ps");
1776 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1777 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1778 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1779 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1780 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1781 printStringNoNewline(&inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1782 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1783 printStringNoNewline(&inp, "mode for center of mass motion removal");
1784 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
1785 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
1786 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
1787 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
1788 setStringEntry(&inp, "comm-grps", is->vcm, nullptr);
1790 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
1791 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
1792 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
1793 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
1796 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
1797 printStringNoNewline(&inp, "Force tolerance and initial step-size");
1798 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
1799 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
1800 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
1801 ir->niter = get_eint(&inp, "niter", 20, wi);
1802 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
1803 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
1804 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
1805 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
1806 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
1808 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
1809 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
1811 /* Output options */
1812 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
1813 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
1814 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
1815 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
1816 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
1817 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
1818 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
1819 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
1820 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
1821 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
1822 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
1823 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
1824 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
1825 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
1826 printStringNoNewline(&inp, "default, all atoms will be written.");
1827 setStringEntry(&inp, "compressed-x-grps", is->x_compressed_groups, nullptr);
1828 printStringNoNewline(&inp, "Selection of energy groups");
1829 setStringEntry(&inp, "energygrps", is->energy, nullptr);
1831 /* Neighbor searching */
1832 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
1833 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1834 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
1835 printStringNoNewline(&inp, "nblist update frequency");
1836 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
1837 printStringNoNewline(&inp, "ns algorithm (simple or grid)");
1838 ir->ns_type = get_eeenum(&inp, "ns-type", ens_names, wi);
1839 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
1840 ir->ePBC = get_eeenum(&inp, "pbc", epbc_names, wi);
1841 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
1842 printStringNoNewline(&inp, "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1843 printStringNoNewline(&inp, "a value of -1 means: use rlist");
1844 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
1845 printStringNoNewline(&inp, "nblist cut-off");
1846 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
1847 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
1849 /* Electrostatics */
1850 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
1851 printStringNoNewline(&inp, "Method for doing electrostatics");
1852 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
1853 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
1854 printStringNoNewline(&inp, "cut-off lengths");
1855 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
1856 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
1857 printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
1858 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
1859 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
1860 printStringNoNewline(&inp, "Method for doing Van der Waals");
1861 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
1862 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
1863 printStringNoNewline(&inp, "cut-off lengths");
1864 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
1865 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
1866 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
1867 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
1868 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
1869 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
1870 printStringNoNewline(&inp, "Separate tables between energy group pairs");
1871 setStringEntry(&inp, "energygrp-table", is->egptable, nullptr);
1872 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
1873 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
1874 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
1875 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
1876 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
1877 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
1878 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
1879 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
1880 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
1881 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
1882 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
1883 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
1884 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
1886 /* Implicit solvation is no longer supported, but we need grompp
1887 to be able to refuse old .mdp files that would have built a tpr
1888 to run it. Thus, only "no" is accepted. */
1889 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
1891 /* Coupling stuff */
1892 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
1893 printStringNoNewline(&inp, "Temperature coupling");
1894 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
1895 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
1896 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
1897 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
1898 printStringNoNewline(&inp, "Groups to couple separately");
1899 setStringEntry(&inp, "tc-grps", is->tcgrps, nullptr);
1900 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
1901 setStringEntry(&inp, "tau-t", is->tau_t, nullptr);
1902 setStringEntry(&inp, "ref-t", is->ref_t, nullptr);
1903 printStringNoNewline(&inp, "pressure coupling");
1904 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
1905 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
1906 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
1907 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
1908 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
1909 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
1910 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
1911 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
1912 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
1915 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
1916 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
1917 printStringNoNewline(&inp, "Groups treated Quantum Mechanically");
1918 setStringEntry(&inp, "QMMM-grps", is->QMMM, nullptr);
1919 printStringNoNewline(&inp, "QM method");
1920 setStringEntry(&inp, "QMmethod", is->QMmethod, nullptr);
1921 printStringNoNewline(&inp, "QMMM scheme");
1922 ir->QMMMscheme = get_eeenum(&inp, "QMMMscheme", eQMMMscheme_names, wi);
1923 printStringNoNewline(&inp, "QM basisset");
1924 setStringEntry(&inp, "QMbasis", is->QMbasis, nullptr);
1925 printStringNoNewline(&inp, "QM charge");
1926 setStringEntry(&inp, "QMcharge", is->QMcharge, nullptr);
1927 printStringNoNewline(&inp, "QM multiplicity");
1928 setStringEntry(&inp, "QMmult", is->QMmult, nullptr);
1929 printStringNoNewline(&inp, "Surface Hopping");
1930 setStringEntry(&inp, "SH", is->bSH, nullptr);
1931 printStringNoNewline(&inp, "CAS space options");
1932 setStringEntry(&inp, "CASorbitals", is->CASorbitals, nullptr);
1933 setStringEntry(&inp, "CASelectrons", is->CASelectrons, nullptr);
1934 setStringEntry(&inp, "SAon", is->SAon, nullptr);
1935 setStringEntry(&inp, "SAoff", is->SAoff, nullptr);
1936 setStringEntry(&inp, "SAsteps", is->SAsteps, nullptr);
1937 printStringNoNewline(&inp, "Scale factor for MM charges");
1938 ir->scalefactor = get_ereal(&inp, "MMChargeScaleFactor", 1.0, wi);
1940 /* Simulated annealing */
1941 printStringNewline(&inp, "SIMULATED ANNEALING");
1942 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
1943 setStringEntry(&inp, "annealing", is->anneal, nullptr);
1944 printStringNoNewline(&inp, "Number of time points to use for specifying annealing in each group");
1945 setStringEntry(&inp, "annealing-npoints", is->anneal_npoints, nullptr);
1946 printStringNoNewline(&inp, "List of times at the annealing points for each group");
1947 setStringEntry(&inp, "annealing-time", is->anneal_time, nullptr);
1948 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
1949 setStringEntry(&inp, "annealing-temp", is->anneal_temp, nullptr);
1952 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
1953 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
1954 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
1955 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
1958 printStringNewline(&inp, "OPTIONS FOR BONDS");
1959 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
1960 printStringNoNewline(&inp, "Type of constraint algorithm");
1961 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
1962 printStringNoNewline(&inp, "Do not constrain the start configuration");
1963 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
1964 printStringNoNewline(&inp, "Use successive overrelaxation to reduce the number of shake iterations");
1965 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
1966 printStringNoNewline(&inp, "Relative tolerance of shake");
1967 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
1968 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
1969 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
1970 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
1971 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
1972 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
1973 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
1974 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
1975 printStringNoNewline(&inp, "rotates over more degrees than");
1976 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
1977 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
1978 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
1980 /* Energy group exclusions */
1981 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
1982 printStringNoNewline(&inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
1983 setStringEntry(&inp, "energygrp-excl", is->egpexcl, nullptr);
1986 printStringNewline(&inp, "WALLS");
1987 printStringNoNewline(&inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1988 ir->nwall = get_eint(&inp, "nwall", 0, wi);
1989 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
1990 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
1991 setStringEntry(&inp, "wall-atomtype", is->wall_atomtype, nullptr);
1992 setStringEntry(&inp, "wall-density", is->wall_density, nullptr);
1993 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
1996 printStringNewline(&inp, "COM PULLING");
1997 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2001 is->pull_grp = read_pullparams(&inp, ir->pull, wi);
2005 NOTE: needs COM pulling input */
2006 printStringNewline(&inp, "AWH biasing");
2007 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2012 ir->awhParams = gmx::readAndCheckAwhParams(&inp, ir, wi);
2016 gmx_fatal(FARGS, "AWH biasing is only compatible with COM pulling turned on");
2020 /* Enforced rotation */
2021 printStringNewline(&inp, "ENFORCED ROTATION");
2022 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2023 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2027 is->rot_grp = read_rotparams(&inp, ir->rot, wi);
2030 /* Interactive MD */
2032 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2033 setStringEntry(&inp, "IMD-group", is->imd_grp, nullptr);
2034 if (is->imd_grp[0] != '\0')
2041 printStringNewline(&inp, "NMR refinement stuff");
2042 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2043 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2044 printStringNoNewline(&inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2045 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2046 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2047 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2048 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2049 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2050 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2051 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2052 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2053 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2054 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2055 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2056 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2057 setStringEntry(&inp, "orire-fitgrp", is->orirefitgrp, nullptr);
2058 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2059 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2061 /* free energy variables */
2062 printStringNewline(&inp, "Free energy variables");
2063 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2064 setStringEntry(&inp, "couple-moltype", is->couple_moltype, nullptr);
2065 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2066 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2067 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2069 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2071 it was not entered */
2072 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2073 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2074 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2075 setStringEntry(&inp, "fep-lambdas", is->fep_lambda[efptFEP], nullptr);
2076 setStringEntry(&inp, "mass-lambdas", is->fep_lambda[efptMASS], nullptr);
2077 setStringEntry(&inp, "coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
2078 setStringEntry(&inp, "vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
2079 setStringEntry(&inp, "bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
2080 setStringEntry(&inp, "restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
2081 setStringEntry(&inp, "temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
2082 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2083 setStringEntry(&inp, "init-lambda-weights", is->lambda_weights, nullptr);
2084 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2085 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2086 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2087 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2088 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2089 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2090 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2091 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2092 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2093 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2094 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2095 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2097 /* Non-equilibrium MD stuff */
2098 printStringNewline(&inp, "Non-equilibrium MD stuff");
2099 setStringEntry(&inp, "acc-grps", is->accgrps, nullptr);
2100 setStringEntry(&inp, "accelerate", is->acc, nullptr);
2101 setStringEntry(&inp, "freezegrps", is->freeze, nullptr);
2102 setStringEntry(&inp, "freezedim", is->frdim, nullptr);
2103 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2104 setStringEntry(&inp, "deform", is->deform, nullptr);
2106 /* simulated tempering variables */
2107 printStringNewline(&inp, "simulated tempering variables");
2108 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2109 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2110 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2111 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2113 /* expanded ensemble variables */
2114 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2116 read_expandedparams(&inp, expand, wi);
2119 /* Electric fields */
2121 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2122 gmx::KeyValueTreeTransformer transform;
2123 transform.rules()->addRule()
2124 .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2125 mdModules->initMdpTransform(transform.rules());
2126 for (const auto &path : transform.mappedPaths())
2128 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2129 mark_einp_set(inp, path[0].c_str());
2131 MdpErrorHandler errorHandler(wi);
2133 = transform.transform(convertedValues, &errorHandler);
2134 ir->params = new gmx::KeyValueTreeObject(result.object());
2135 mdModules->adjustInputrecBasedOnModules(ir);
2136 errorHandler.setBackMapping(result.backMapping());
2137 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2140 /* Ion/water position swapping ("computational electrophysiology") */
2141 printStringNewline(&inp, "Ion/water position swapping for computational electrophysiology setups");
2142 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2143 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2144 if (ir->eSwapCoords != eswapNO)
2151 printStringNoNewline(&inp, "Swap attempt frequency");
2152 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2153 printStringNoNewline(&inp, "Number of ion types to be controlled");
2154 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2157 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2159 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2160 snew(ir->swap->grp, ir->swap->ngrp);
2161 for (i = 0; i < ir->swap->ngrp; i++)
2163 snew(ir->swap->grp[i].molname, STRLEN);
2165 printStringNoNewline(&inp, "Two index groups that contain the compartment-partitioning atoms");
2166 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2167 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2168 printStringNoNewline(&inp, "Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2169 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2170 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2172 printStringNoNewline(&inp, "Name of solvent molecules");
2173 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2175 printStringNoNewline(&inp, "Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2176 printStringNoNewline(&inp, "Note that the split cylinder settings do not have an influence on the swapping protocol,");
2177 printStringNoNewline(&inp, "however, if correctly defined, the permeation events are recorded per channel");
2178 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2179 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2180 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2181 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2182 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2183 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2185 printStringNoNewline(&inp, "Average the number of ions per compartment over these many swap attempt steps");
2186 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2188 printStringNoNewline(&inp, "Names of the ion types that can be exchanged with solvent molecules,");
2189 printStringNoNewline(&inp, "and the requested number of ions of this type in compartments A and B");
2190 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2191 for (i = 0; i < nIonTypes; i++)
2193 int ig = eSwapFixedGrpNR + i;
2195 sprintf(buf, "iontype%d-name", i);
2196 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2197 sprintf(buf, "iontype%d-in-A", i);
2198 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2199 sprintf(buf, "iontype%d-in-B", i);
2200 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2203 printStringNoNewline(&inp, "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2204 printStringNoNewline(&inp, "at maximum distance (= bulk concentration) to the split group layers. However,");
2205 printStringNoNewline(&inp, "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2206 printStringNoNewline(&inp, "towards one of the compartment-partitioning layers (at +/- 1.0).");
2207 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2208 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2209 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2210 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2212 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2215 printStringNoNewline(&inp, "Start to swap ions if threshold difference to requested count is reached");
2216 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2219 /* AdResS is no longer supported, but we need grompp to be able to
2220 refuse to process old .mdp files that used it. */
2221 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2223 /* User defined thingies */
2224 printStringNewline(&inp, "User defined thingies");
2225 setStringEntry(&inp, "user1-grps", is->user1, nullptr);
2226 setStringEntry(&inp, "user2-grps", is->user2, nullptr);
2227 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2228 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2229 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2230 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2231 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2232 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2233 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2234 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2238 gmx::TextOutputFile stream(mdparout);
2239 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2241 // Transform module data into a flat key-value tree for output.
2242 gmx::KeyValueTreeBuilder builder;
2243 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2244 mdModules->buildMdpOutput(&builderObject);
2246 gmx::TextWriter writer(&stream);
2247 writeKeyValueTreeAsMdp(&writer, builder.build());
2252 /* Process options if necessary */
2253 for (m = 0; m < 2; m++)
2255 for (i = 0; i < 2*DIM; i++)
2264 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2266 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
2268 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2270 case epctSEMIISOTROPIC:
2271 case epctSURFACETENSION:
2272 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2274 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
2276 dumdub[m][YY] = dumdub[m][XX];
2278 case epctANISOTROPIC:
2279 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2280 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2281 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2283 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
2287 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2288 epcoupltype_names[ir->epct]);
2292 clear_mat(ir->ref_p);
2293 clear_mat(ir->compress);
2294 for (i = 0; i < DIM; i++)
2296 ir->ref_p[i][i] = dumdub[1][i];
2297 ir->compress[i][i] = dumdub[0][i];
2299 if (ir->epct == epctANISOTROPIC)
2301 ir->ref_p[XX][YY] = dumdub[1][3];
2302 ir->ref_p[XX][ZZ] = dumdub[1][4];
2303 ir->ref_p[YY][ZZ] = dumdub[1][5];
2304 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2306 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2308 ir->compress[XX][YY] = dumdub[0][3];
2309 ir->compress[XX][ZZ] = dumdub[0][4];
2310 ir->compress[YY][ZZ] = dumdub[0][5];
2311 for (i = 0; i < DIM; i++)
2313 for (m = 0; m < i; m++)
2315 ir->ref_p[i][m] = ir->ref_p[m][i];
2316 ir->compress[i][m] = ir->compress[m][i];
2321 if (ir->comm_mode == ecmNO)
2326 opts->couple_moltype = nullptr;
2327 if (strlen(is->couple_moltype) > 0)
2329 if (ir->efep != efepNO)
2331 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2332 if (opts->couple_lam0 == opts->couple_lam1)
2334 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2336 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2337 opts->couple_lam1 == ecouplamNONE))
2339 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2344 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2347 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2348 if (ir->efep != efepNO)
2350 if (fep->delta_lambda > 0)
2352 ir->efep = efepSLOWGROWTH;
2356 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2358 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2359 warning_note(wi, "Old option for dhdl-print-energy given: "
2360 "changing \"yes\" to \"total\"\n");
2363 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2365 /* always print out the energy to dhdl if we are doing
2366 expanded ensemble, since we need the total energy for
2367 analysis if the temperature is changing. In some
2368 conditions one may only want the potential energy, so
2369 we will allow that if the appropriate mdp setting has
2370 been enabled. Otherwise, total it is:
2372 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2375 if ((ir->efep != efepNO) || ir->bSimTemp)
2377 ir->bExpanded = FALSE;
2378 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2380 ir->bExpanded = TRUE;
2382 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2383 if (ir->bSimTemp) /* done after fep params */
2385 do_simtemp_params(ir);
2388 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2389 * setup and not on the old way of specifying the free-energy setup,
2390 * we should check for using soft-core when not needed, since that
2391 * can complicate the sampling significantly.
2392 * Note that we only check for the automated coupling setup.
2393 * If the (advanced) user does FEP through manual topology changes,
2394 * this check will not be triggered.
2396 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2397 ir->fepvals->sc_alpha != 0 &&
2398 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2399 couple_lambda_has_vdw_on(opts->couple_lam1)))
2401 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.");
2406 ir->fepvals->n_lambda = 0;
2409 /* WALL PARAMETERS */
2411 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts, wi);
2413 /* ORIENTATION RESTRAINT PARAMETERS */
2415 if (opts->bOrire && gmx::splitString(is->orirefitgrp).size() != 1)
2417 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2420 /* DEFORMATION PARAMETERS */
2422 clear_mat(ir->deform);
2423 for (i = 0; i < 6; i++)
2428 double gmx_unused canary;
2429 int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
2430 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2431 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
2433 if (strlen(is->deform) > 0 && ndeform != 6)
2435 warning_error(wi, gmx::formatString("Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform).c_str());
2437 for (i = 0; i < 3; i++)
2439 ir->deform[i][i] = dumdub[0][i];
2441 ir->deform[YY][XX] = dumdub[0][3];
2442 ir->deform[ZZ][XX] = dumdub[0][4];
2443 ir->deform[ZZ][YY] = dumdub[0][5];
2444 if (ir->epc != epcNO)
2446 for (i = 0; i < 3; i++)
2448 for (j = 0; j <= i; j++)
2450 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2452 warning_error(wi, "A box element has deform set and compressibility > 0");
2456 for (i = 0; i < 3; i++)
2458 for (j = 0; j < i; j++)
2460 if (ir->deform[i][j] != 0)
2462 for (m = j; m < DIM; m++)
2464 if (ir->compress[m][j] != 0)
2466 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.");
2467 warning(wi, warn_buf);
2475 /* Ion/water position swapping checks */
2476 if (ir->eSwapCoords != eswapNO)
2478 if (ir->swap->nstswap < 1)
2480 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2482 if (ir->swap->nAverage < 1)
2484 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2486 if (ir->swap->threshold < 1.0)
2488 warning_error(wi, "Ion count threshold must be at least 1.\n");
2496 static int search_QMstring(const char *s, int ng, const char *gn[])
2498 /* same as normal search_string, but this one searches QM strings */
2501 for (i = 0; (i < ng); i++)
2503 if (gmx_strcasecmp(s, gn[i]) == 0)
2509 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2510 } /* search_QMstring */
2512 /* We would like gn to be const as well, but C doesn't allow this */
2513 /* TODO this is utility functionality (search for the index of a
2514 string in a collection), so should be refactored and located more
2516 int search_string(const char *s, int ng, char *gn[])
2520 for (i = 0; (i < ng); i++)
2522 if (gmx_strcasecmp(s, gn[i]) == 0)
2529 "Group %s referenced in the .mdp file was not found in the index file.\n"
2530 "Group names must match either [moleculetype] names or custom index group\n"
2531 "names, in which case you must supply an index file to the '-n' option\n"
2536 static bool do_numbering(int natoms, SimulationGroups *groups,
2537 gmx::ArrayRef<std::string> groupsFromMdpFile,
2538 t_blocka *block, char *gnames[],
2539 SimulationAtomGroupType gtype, int restnm,
2540 int grptp, bool bVerbose,
2543 unsigned short *cbuf;
2544 AtomGroupIndices *grps = &(groups->groups[gtype]);
2545 int j, gid, aj, ognr, ntot = 0;
2548 char warn_buf[STRLEN];
2550 title = shortName(gtype);
2553 /* Mark all id's as not set */
2554 for (int i = 0; (i < natoms); i++)
2559 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2561 /* Lookup the group name in the block structure */
2562 gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2563 if ((grptp != egrptpONE) || (i == 0))
2565 grps->emplace_back(gid);
2568 /* Now go over the atoms in the group */
2569 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2574 /* Range checking */
2575 if ((aj < 0) || (aj >= natoms))
2577 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2579 /* Lookup up the old group number */
2583 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2584 aj+1, title, ognr+1, i+1);
2588 /* Store the group number in buffer */
2589 if (grptp == egrptpONE)
2602 /* Now check whether we have done all atoms */
2606 if (grptp == egrptpALL)
2608 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2609 natoms-ntot, title);
2611 else if (grptp == egrptpPART)
2613 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2614 natoms-ntot, title);
2615 warning_note(wi, warn_buf);
2617 /* Assign all atoms currently unassigned to a rest group */
2618 for (j = 0; (j < natoms); j++)
2620 if (cbuf[j] == NOGID)
2622 cbuf[j] = grps->size();
2626 if (grptp != egrptpPART)
2631 "Making dummy/rest group for %s containing %d elements\n",
2632 title, natoms-ntot);
2634 /* Add group name "rest" */
2635 grps->emplace_back(restnm);
2637 /* Assign the rest name to all atoms not currently assigned to a group */
2638 for (j = 0; (j < natoms); j++)
2640 if (cbuf[j] == NOGID)
2642 // group size was not updated before this here, so need to use -1.
2643 cbuf[j] = grps->size() - 1;
2649 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2651 /* All atoms are part of one (or no) group, no index required */
2652 groups->groupNumbers[gtype].clear();
2656 for (int j = 0; (j < natoms); j++)
2658 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2664 return (bRest && grptp == egrptpPART);
2667 static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2670 pull_params_t *pull;
2671 int natoms, imin, jmin;
2672 int *nrdf2, *na_vcm, na_tot;
2673 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2678 * First calc 3xnr-atoms for each group
2679 * then subtract half a degree of freedom for each constraint
2681 * Only atoms and nuclei contribute to the degrees of freedom...
2686 const SimulationGroups &groups = mtop->groups;
2687 natoms = mtop->natoms;
2689 /* Allocate one more for a possible rest group */
2690 /* We need to sum degrees of freedom into doubles,
2691 * since floats give too low nrdf's above 3 million atoms.
2693 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size()+1);
2694 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2695 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2696 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2697 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2699 for (int i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2703 for (int i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; i++)
2706 clear_ivec(dof_vcm[i]);
2708 nrdf_vcm_sub[i] = 0;
2710 snew(nrdf2, natoms);
2711 for (const AtomProxy atomP : AtomRange(*mtop))
2713 const t_atom &local = atomP.atom();
2714 int i = atomP.globalAtomNumber();
2716 if (local.ptype == eptAtom || local.ptype == eptNucleus)
2718 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
2719 for (int d = 0; d < DIM; d++)
2721 if (opts->nFreeze[g][d] == 0)
2723 /* Add one DOF for particle i (counted as 2*1) */
2725 /* VCM group i has dim d as a DOF */
2726 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] = 1;
2729 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] += 0.5*nrdf2[i];
2730 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] += 0.5*nrdf2[i];
2735 for (const gmx_molblock_t &molb : mtop->molblock)
2737 const gmx_moltype_t &molt = mtop->moltype[molb.type];
2738 const t_atom *atom = molt.atoms.atom;
2739 for (int mol = 0; mol < molb.nmol; mol++)
2741 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2743 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
2744 for (int i = 0; i < molt.ilist[ftype].size(); )
2746 /* Subtract degrees of freedom for the constraints,
2747 * if the particles still have degrees of freedom left.
2748 * If one of the particles is a vsite or a shell, then all
2749 * constraint motion will go there, but since they do not
2750 * contribute to the constraints the degrees of freedom do not
2753 int ai = as + ia[i + 1];
2754 int aj = as + ia[i + 2];
2755 if (((atom[ia[i + 1]].ptype == eptNucleus) ||
2756 (atom[ia[i + 1]].ptype == eptAtom)) &&
2757 ((atom[ia[i + 2]].ptype == eptNucleus) ||
2758 (atom[ia[i + 2]].ptype == eptAtom)))
2776 imin = std::min(imin, nrdf2[ai]);
2777 jmin = std::min(jmin, nrdf2[aj]);
2780 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2781 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -= 0.5*jmin;
2782 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2783 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -= 0.5*jmin;
2785 i += interaction_function[ftype].nratoms+1;
2788 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
2789 for (int i = 0; i < molt.ilist[F_SETTLE].size(); )
2791 /* Subtract 1 dof from every atom in the SETTLE */
2792 for (int j = 0; j < 3; j++)
2794 int ai = as + ia[i + 1 + j];
2795 imin = std::min(2, nrdf2[ai]);
2797 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2798 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2802 as += molt.atoms.nr;
2808 /* Correct nrdf for the COM constraints.
2809 * We correct using the TC and VCM group of the first atom
2810 * in the reference and pull group. If atoms in one pull group
2811 * belong to different TC or VCM groups it is anyhow difficult
2812 * to determine the optimal nrdf assignment.
2816 for (int i = 0; i < pull->ncoord; i++)
2818 if (pull->coord[i].eType != epullCONSTRAINT)
2825 for (int j = 0; j < 2; j++)
2827 const t_pull_group *pgrp;
2829 pgrp = &pull->group[pull->coord[i].group[j]];
2833 /* Subtract 1/2 dof from each group */
2834 int ai = pgrp->ind[0];
2835 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2836 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2837 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
2839 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)]]);
2844 /* We need to subtract the whole DOF from group j=1 */
2851 if (ir->nstcomm != 0)
2855 /* We remove COM motion up to dim ndof_com() */
2856 ndim_rm_vcm = ndof_com(ir);
2858 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
2859 * the number of degrees of freedom in each vcm group when COM
2860 * translation is removed and 6 when rotation is removed as well.
2862 for (int j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2864 switch (ir->comm_mode)
2867 case ecmLINEAR_ACCELERATION_CORRECTION:
2868 nrdf_vcm_sub[j] = 0;
2869 for (int d = 0; d < ndim_rm_vcm; d++)
2878 nrdf_vcm_sub[j] = 6;
2881 gmx_incons("Checking comm_mode");
2885 for (int i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2887 /* Count the number of atoms of TC group i for every VCM group */
2888 for (int j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2893 for (int ai = 0; ai < natoms; ai++)
2895 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
2897 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
2901 /* Correct for VCM removal according to the fraction of each VCM
2902 * group present in this TC group.
2904 nrdf_uc = nrdf_tc[i];
2906 for (int j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2908 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
2910 nrdf_tc[i] += nrdf_uc*(static_cast<double>(na_vcm[j])/static_cast<double>(na_tot))*
2911 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
2916 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
2918 opts->nrdf[i] = nrdf_tc[i];
2919 if (opts->nrdf[i] < 0)
2924 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2925 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
2933 sfree(nrdf_vcm_sub);
2936 static bool do_egp_flag(t_inputrec *ir, SimulationGroups *groups,
2937 const char *option, const char *val, int flag)
2939 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2940 * But since this is much larger than STRLEN, such a line can not be parsed.
2941 * The real maximum is the number of names that fit in a string: STRLEN/2.
2943 #define EGP_MAX (STRLEN/2)
2947 auto names = gmx::splitString(val);
2948 if (names.size() % 2 != 0)
2950 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2952 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
2954 for (size_t i = 0; i < names.size() / 2; i++)
2956 // TODO this needs to be replaced by a solution using std::find_if
2959 gmx_strcasecmp(names[2*i].c_str(), *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
2965 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2966 names[2*i].c_str(), option);
2970 gmx_strcasecmp(names[2*i+1].c_str(), *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
2976 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2977 names[2*i+1].c_str(), option);
2979 if ((j < nr) && (k < nr))
2981 ir->opts.egp_flags[nr*j+k] |= flag;
2982 ir->opts.egp_flags[nr*k+j] |= flag;
2991 static void make_swap_groups(
2996 int ig = -1, i = 0, gind;
3000 /* Just a quick check here, more thorough checks are in mdrun */
3001 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3003 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3006 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3007 for (ig = 0; ig < swap->ngrp; ig++)
3009 swapg = &swap->grp[ig];
3010 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3011 swapg->nat = grps->index[gind+1] - grps->index[gind];
3015 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3016 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3017 swap->grp[ig].molname, swapg->nat);
3018 snew(swapg->ind, swapg->nat);
3019 for (i = 0; i < swapg->nat; i++)
3021 swapg->ind[i] = grps->a[grps->index[gind]+i];
3026 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3032 static void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3037 ig = search_string(IMDgname, grps->nr, gnames);
3038 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3040 if (IMDgroup->nat > 0)
3042 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3043 IMDgname, IMDgroup->nat);
3044 snew(IMDgroup->ind, IMDgroup->nat);
3045 for (i = 0; i < IMDgroup->nat; i++)
3047 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3052 void do_index(const char* mdparin, const char *ndx,
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 auto accelerations = gmx::splitString(is->acc);
3402 auto accelerationGroupNames = gmx::splitString(is->accgrps);
3403 if (accelerationGroupNames.size() * DIM != accelerations.size())
3405 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3406 accelerationGroupNames.size(), accelerations.size());
3408 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3409 SimulationAtomGroupType::Acceleration,
3410 restnm, egrptpALL_GENREST, bVerbose, wi);
3411 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3412 snew(ir->opts.acc, nr);
3413 ir->opts.ngacc = nr;
3415 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3417 auto freezeDims = gmx::splitString(is->frdim);
3418 auto freezeGroupNames = gmx::splitString(is->freeze);
3419 if (freezeDims.size() != DIM * freezeGroupNames.size())
3421 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3422 freezeGroupNames.size(), freezeDims.size());
3424 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3425 SimulationAtomGroupType::Freeze,
3426 restnm, egrptpALL_GENREST, bVerbose, wi);
3427 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3428 ir->opts.ngfrz = nr;
3429 snew(ir->opts.nFreeze, nr);
3430 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3432 for (j = 0; (j < DIM); j++, k++)
3434 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3435 if (!ir->opts.nFreeze[i][j])
3437 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3439 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3440 "(not %s)", freezeDims[k].c_str());
3441 warning(wi, warn_buf);
3446 for (; (i < nr); i++)
3448 for (j = 0; (j < DIM); j++)
3450 ir->opts.nFreeze[i][j] = 0;
3454 auto energyGroupNames = gmx::splitString(is->energy);
3455 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3456 SimulationAtomGroupType::EnergyOutput,
3457 restnm, egrptpALL_GENREST, bVerbose, wi);
3458 add_wall_energrps(groups, ir->nwall, symtab);
3459 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3460 auto vcmGroupNames = gmx::splitString(is->vcm);
3462 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3463 SimulationAtomGroupType::MassCenterVelocityRemoval,
3464 restnm, vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3467 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3468 "This may lead to artifacts.\n"
3469 "In most cases one should use one group for the whole system.");
3472 /* Now we have filled the freeze struct, so we can calculate NRDF */
3473 calc_nrdf(mtop, ir, gnames);
3475 auto user1GroupNames = gmx::splitString(is->user1);
3476 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3477 SimulationAtomGroupType::User1,
3478 restnm, egrptpALL_GENREST, bVerbose, wi);
3479 auto user2GroupNames = gmx::splitString(is->user2);
3480 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3481 SimulationAtomGroupType::User2,
3482 restnm, egrptpALL_GENREST, bVerbose, wi);
3483 auto compressedXGroupNames = gmx::splitString(is->x_compressed_groups);
3484 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3485 SimulationAtomGroupType::CompressedPositionOutput,
3486 restnm, egrptpONE, bVerbose, wi);
3487 auto orirefFitGroupNames = gmx::splitString(is->orirefitgrp);
3488 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3489 SimulationAtomGroupType::OrientationRestraintsFit,
3490 restnm, egrptpALL_GENREST, bVerbose, wi);
3492 /* QMMM input processing */
3493 auto qmGroupNames = gmx::splitString(is->QMMM);
3494 auto qmMethods = gmx::splitString(is->QMmethod);
3495 auto qmBasisSets = gmx::splitString(is->QMbasis);
3496 if (ir->eI != eiMimic)
3498 if (qmMethods.size() != qmGroupNames.size() ||
3499 qmBasisSets.size() != qmGroupNames.size())
3501 gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets"
3502 " and %zu methods\n", qmGroupNames.size(),
3503 qmBasisSets.size(), qmMethods.size());
3505 /* group rest, if any, is always MM! */
3506 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3507 SimulationAtomGroupType::QuantumMechanics,
3508 restnm, egrptpALL_GENREST, bVerbose, wi);
3509 nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
3510 ir->opts.ngQM = qmGroupNames.size();
3511 snew(ir->opts.QMmethod, nr);
3512 snew(ir->opts.QMbasis, nr);
3513 for (i = 0; i < nr; i++)
3515 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3516 * converted to the corresponding enum in names.c
3518 ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(),
3521 ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(),
3526 auto qmMultiplicities = gmx::splitString(is->QMmult);
3527 auto qmCharges = gmx::splitString(is->QMcharge);
3528 auto qmbSH = gmx::splitString(is->bSH);
3529 snew(ir->opts.QMmult, nr);
3530 snew(ir->opts.QMcharge, nr);
3531 snew(ir->opts.bSH, nr);
3532 convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult);
3533 convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge);
3534 convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH);
3536 auto CASelectrons = gmx::splitString(is->CASelectrons);
3537 auto CASorbitals = gmx::splitString(is->CASorbitals);
3538 snew(ir->opts.CASelectrons, nr);
3539 snew(ir->opts.CASorbitals, nr);
3540 convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons);
3541 convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals);
3543 auto SAon = gmx::splitString(is->SAon);
3544 auto SAoff = gmx::splitString(is->SAoff);
3545 auto SAsteps = gmx::splitString(is->SAsteps);
3546 snew(ir->opts.SAon, nr);
3547 snew(ir->opts.SAoff, nr);
3548 snew(ir->opts.SAsteps, nr);
3549 convertInts(wi, SAon, "SAon", ir->opts.SAon);
3550 convertInts(wi, SAoff, "SAoff", ir->opts.SAoff);
3551 convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps);
3556 if (qmGroupNames.size() > 1)
3558 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3560 /* group rest, if any, is always MM! */
3561 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3562 SimulationAtomGroupType::QuantumMechanics,
3563 restnm, egrptpALL_GENREST, bVerbose, wi);
3565 ir->opts.ngQM = qmGroupNames.size();
3568 /* end of QMMM input */
3572 for (auto group : gmx::keysOf(groups->groups))
3574 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3575 for (const auto &entry : groups->groups[group])
3577 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3579 fprintf(stderr, "\n");
3583 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3584 snew(ir->opts.egp_flags, nr*nr);
3586 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3587 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3589 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3591 if (bExcl && EEL_FULL(ir->coulombtype))
3593 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3596 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3597 if (bTable && !(ir->vdwtype == evdwUSER) &&
3598 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3599 !(ir->coulombtype == eelPMEUSERSWITCH))
3601 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3604 /* final check before going out of scope if simulated tempering variables
3605 * need to be set to default values.
3607 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3609 ir->expandedvals->nstexpanded = 2*static_cast<int>(ir->opts.tau_t[0]/ir->delta_t);
3610 warning(wi, gmx::formatString("the value for nstexpanded was not specified for "
3611 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3612 "by default, but it is recommended to set it to an explicit value!",
3613 ir->expandedvals->nstexpanded));
3615 for (i = 0; (i < defaultIndexGroups->nr); i++)
3620 done_blocka(defaultIndexGroups);
3621 sfree(defaultIndexGroups);
3627 static void check_disre(const gmx_mtop_t *mtop)
3629 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3631 const gmx_ffparams_t &ffparams = mtop->ffparams;
3634 for (int i = 0; i < ffparams.numTypes(); i++)
3636 int ftype = ffparams.functype[i];
3637 if (ftype == F_DISRES)
3639 int label = ffparams.iparams[i].disres.label;
3640 if (label == old_label)
3642 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3650 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3651 "probably the parameters for multiple pairs in one restraint "
3652 "are not identical\n", ndouble);
3657 static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3662 gmx_mtop_ilistloop_t iloop;
3671 for (d = 0; d < DIM; d++)
3673 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3675 /* Check for freeze groups */
3676 for (g = 0; g < ir->opts.ngfrz; g++)
3678 for (d = 0; d < DIM; d++)
3680 if (ir->opts.nFreeze[g][d] != 0)
3688 /* Check for position restraints */
3689 iloop = gmx_mtop_ilistloop_init(sys);
3690 while (const InteractionLists *ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
3693 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3695 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
3697 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
3698 for (d = 0; d < DIM; d++)
3700 if (pr->posres.fcA[d] != 0)
3706 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
3708 /* Check for flat-bottom posres */
3709 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
3710 if (pr->fbposres.k != 0)
3712 switch (pr->fbposres.geom)
3714 case efbposresSPHERE:
3715 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3717 case efbposresCYLINDERX:
3718 AbsRef[YY] = AbsRef[ZZ] = 1;
3720 case efbposresCYLINDERY:
3721 AbsRef[XX] = AbsRef[ZZ] = 1;
3723 case efbposresCYLINDER:
3724 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3725 case efbposresCYLINDERZ:
3726 AbsRef[XX] = AbsRef[YY] = 1;
3728 case efbposresX: /* d=XX */
3729 case efbposresY: /* d=YY */
3730 case efbposresZ: /* d=ZZ */
3731 d = pr->fbposres.geom - efbposresX;
3735 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3736 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3744 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3748 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3749 bool *bC6ParametersWorkWithGeometricRules,
3750 bool *bC6ParametersWorkWithLBRules,
3751 bool *bLBRulesPossible)
3753 int ntypes, tpi, tpj;
3756 double c6i, c6j, c12i, c12j;
3757 double c6, c6_geometric, c6_LB;
3758 double sigmai, sigmaj, epsi, epsj;
3759 bool bCanDoLBRules, bCanDoGeometricRules;
3762 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3763 * force-field floating point parameters.
3766 ptr = getenv("GMX_LJCOMB_TOL");
3770 double gmx_unused canary;
3772 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3774 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
3779 *bC6ParametersWorkWithLBRules = TRUE;
3780 *bC6ParametersWorkWithGeometricRules = TRUE;
3781 bCanDoLBRules = TRUE;
3782 ntypes = mtop->ffparams.atnr;
3783 snew(typecount, ntypes);
3784 gmx_mtop_count_atomtypes(mtop, state, typecount);
3785 *bLBRulesPossible = TRUE;
3786 for (tpi = 0; tpi < ntypes; ++tpi)
3788 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3789 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3790 for (tpj = tpi; tpj < ntypes; ++tpj)
3792 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3793 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3794 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3795 c6_geometric = std::sqrt(c6i * c6j);
3796 if (!gmx_numzero(c6_geometric))
3798 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3800 sigmai = gmx::sixthroot(c12i / c6i);
3801 sigmaj = gmx::sixthroot(c12j / c6j);
3802 epsi = c6i * c6i /(4.0 * c12i);
3803 epsj = c6j * c6j /(4.0 * c12j);
3804 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
3808 *bLBRulesPossible = FALSE;
3809 c6_LB = c6_geometric;
3811 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3816 *bC6ParametersWorkWithLBRules = FALSE;
3819 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3821 if (!bCanDoGeometricRules)
3823 *bC6ParametersWorkWithGeometricRules = FALSE;
3831 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3834 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3836 check_combination_rule_differences(mtop, 0,
3837 &bC6ParametersWorkWithGeometricRules,
3838 &bC6ParametersWorkWithLBRules,
3840 if (ir->ljpme_combination_rule == eljpmeLB)
3842 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
3844 warning(wi, "You are using arithmetic-geometric combination rules "
3845 "in LJ-PME, but your non-bonded C6 parameters do not "
3846 "follow these rules.");
3851 if (!bC6ParametersWorkWithGeometricRules)
3853 if (ir->eDispCorr != edispcNO)
3855 warning_note(wi, "You are using geometric combination rules in "
3856 "LJ-PME, but your non-bonded C6 parameters do "
3857 "not follow these rules. "
3858 "This will introduce very small errors in the forces and energies in "
3859 "your simulations. Dispersion correction will correct total energy "
3860 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3864 warning_note(wi, "You are using geometric combination rules in "
3865 "LJ-PME, but your non-bonded C6 parameters do "
3866 "not follow these rules. "
3867 "This will introduce very small errors in the forces and energies in "
3868 "your simulations. If your system is homogeneous, consider using dispersion correction "
3869 "for the total energy and pressure.");
3875 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3878 char err_buf[STRLEN];
3883 gmx_mtop_atomloop_block_t aloopb;
3885 char warn_buf[STRLEN];
3887 set_warning_line(wi, mdparin, -1);
3889 if (ir->cutoff_scheme == ecutsVERLET &&
3890 ir->verletbuf_tol > 0 &&
3892 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
3893 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
3895 /* Check if a too small Verlet buffer might potentially
3896 * cause more drift than the thermostat can couple off.
3898 /* Temperature error fraction for warning and suggestion */
3899 const real T_error_warn = 0.002;
3900 const real T_error_suggest = 0.001;
3901 /* For safety: 2 DOF per atom (typical with constraints) */
3902 const real nrdf_at = 2;
3903 real T, tau, max_T_error;
3908 for (i = 0; i < ir->opts.ngtc; i++)
3910 T = std::max(T, ir->opts.ref_t[i]);
3911 tau = std::max(tau, ir->opts.tau_t[i]);
3915 /* This is a worst case estimate of the temperature error,
3916 * assuming perfect buffer estimation and no cancelation
3917 * of errors. The factor 0.5 is because energy distributes
3918 * equally over Ekin and Epot.
3920 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
3921 if (max_T_error > T_error_warn)
3923 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.",
3924 ir->verletbuf_tol, T, tau,
3926 100*T_error_suggest,
3927 ir->verletbuf_tol*T_error_suggest/max_T_error);
3928 warning(wi, warn_buf);
3933 if (ETC_ANDERSEN(ir->etc))
3937 for (i = 0; i < ir->opts.ngtc; i++)
3939 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
3940 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
3941 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
3942 i, ir->opts.tau_t[i]);
3943 CHECK(ir->opts.tau_t[i] < 0);
3946 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
3948 for (i = 0; i < ir->opts.ngtc; i++)
3950 int nsteps = gmx::roundToInt(ir->opts.tau_t[i]/ir->delta_t);
3951 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);
3952 CHECK(nsteps % ir->nstcomm != 0);
3957 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3958 ir->comm_mode == ecmNO &&
3959 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
3960 !ETC_ANDERSEN(ir->etc))
3962 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");
3965 if (ir->epc == epcPARRINELLORAHMAN &&
3966 ir->etc == etcNOSEHOOVER)
3969 for (int g = 0; g < ir->opts.ngtc; g++)
3971 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
3973 if (ir->tau_p < 1.9*tau_t_max)
3975 std::string message =
3976 gmx::formatString("With %s T-coupling and %s p-coupling, "
3977 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
3978 etcoupl_names[ir->etc],
3979 epcoupl_names[ir->epc],
3981 "tau-t", tau_t_max);
3982 warning(wi, message.c_str());
3986 /* Check for pressure coupling with absolute position restraints */
3987 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3989 absolute_reference(ir, sys, TRUE, AbsRef);
3991 for (m = 0; m < DIM; m++)
3993 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3995 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4003 aloopb = gmx_mtop_atomloop_block_init(sys);
4005 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4007 if (atom->q != 0 || atom->qB != 0)
4015 if (EEL_FULL(ir->coulombtype))
4018 "You are using full electrostatics treatment %s for a system without charges.\n"
4019 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4020 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4021 warning(wi, err_buf);
4026 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4029 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4030 "You might want to consider using %s electrostatics.\n",
4032 warning_note(wi, err_buf);
4036 /* Check if combination rules used in LJ-PME are the same as in the force field */
4037 if (EVDW_PME(ir->vdwtype))
4039 check_combination_rules(ir, sys, wi);
4042 /* Generalized reaction field */
4043 if (ir->coulombtype == eelGRF_NOTUSED)
4045 warning_error(wi, "Generalized reaction-field electrostatics is no longer supported. "
4046 "You can use normal reaction-field instead and compute the reaction-field constant by hand.");
4050 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4052 for (m = 0; (m < DIM); m++)
4054 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4063 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4064 for (const AtomProxy atomP : AtomRange(*sys))
4066 const t_atom &local = atomP.atom();
4067 int i = atomP.globalAtomNumber();
4068 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4071 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4073 for (m = 0; (m < DIM); m++)
4075 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4079 for (m = 0; (m < DIM); m++)
4081 if (fabs(acc[m]) > 1e-6)
4083 const char *dim[DIM] = { "X", "Y", "Z" };
4085 "Net Acceleration in %s direction, will %s be corrected\n",
4086 dim[m], ir->nstcomm != 0 ? "" : "not");
4087 if (ir->nstcomm != 0 && m < ndof_com(ir))
4090 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4092 ir->opts.acc[i][m] -= acc[m];
4100 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4101 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4103 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4111 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4113 if (ir->pull->coord[i].group[0] == 0 ||
4114 ir->pull->coord[i].group[1] == 0)
4116 absolute_reference(ir, sys, FALSE, AbsRef);
4117 for (m = 0; m < DIM; m++)
4119 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4121 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.");
4129 for (i = 0; i < 3; i++)
4131 for (m = 0; m <= i; m++)
4133 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4134 ir->deform[i][m] != 0)
4136 for (c = 0; c < ir->pull->ncoord; c++)
4138 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4139 ir->pull->coord[c].vec[m] != 0)
4141 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4152 void double_check(t_inputrec *ir, matrix box,
4153 bool bHasNormalConstraints,
4154 bool bHasAnyConstraints,
4158 char warn_buf[STRLEN];
4161 ptr = check_box(ir->ePBC, box);
4164 warning_error(wi, ptr);
4167 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4169 if (ir->shake_tol <= 0.0)
4171 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4173 warning_error(wi, warn_buf);
4177 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4179 /* If we have Lincs constraints: */
4180 if (ir->eI == eiMD && ir->etc == etcNO &&
4181 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4183 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4184 warning_note(wi, warn_buf);
4187 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4189 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4190 warning_note(wi, warn_buf);
4192 if (ir->epc == epcMTTK)
4194 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4198 if (bHasAnyConstraints && ir->epc == epcMTTK)
4200 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4203 if (ir->LincsWarnAngle > 90.0)
4205 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4206 warning(wi, warn_buf);
4207 ir->LincsWarnAngle = 90.0;
4210 if (ir->ePBC != epbcNONE)
4212 if (ir->nstlist == 0)
4214 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4216 if (ir->ns_type == ensGRID)
4218 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4220 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");
4221 warning_error(wi, warn_buf);
4226 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4227 if (2*ir->rlist >= min_size)
4229 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4230 warning_error(wi, warn_buf);
4233 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4240 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4244 real rvdw1, rvdw2, rcoul1, rcoul2;
4245 char warn_buf[STRLEN];
4247 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4251 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4256 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4262 if (rvdw1 + rvdw2 > ir->rlist ||
4263 rcoul1 + rcoul2 > ir->rlist)
4266 "The sum of the two largest charge group radii (%f) "
4267 "is larger than rlist (%f)\n",
4268 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4269 warning(wi, warn_buf);
4273 /* Here we do not use the zero at cut-off macro,
4274 * since user defined interactions might purposely
4275 * not be zero at the cut-off.
4277 if (ir_vdw_is_zero_at_cutoff(ir) &&
4278 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4280 sprintf(warn_buf, "The sum of the two largest charge group "
4281 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4282 "With exact cut-offs, better performance can be "
4283 "obtained with cutoff-scheme = %s, because it "
4284 "does not use charge groups at all.",
4286 ir->rlist, ir->rvdw,
4287 ecutscheme_names[ecutsVERLET]);
4290 warning(wi, warn_buf);
4294 warning_note(wi, warn_buf);
4297 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4298 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4300 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4301 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4303 ir->rlist, ir->rcoulomb,
4304 ecutscheme_names[ecutsVERLET]);
4307 warning(wi, warn_buf);
4311 warning_note(wi, warn_buf);