2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/awh/read_params.h"
50 #include "gromacs/fileio/readinp.h"
51 #include "gromacs/fileio/warninp.h"
52 #include "gromacs/gmxlib/chargegroup.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/gmxpreprocess/toputil.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/calc_verletbuf.h"
59 #include "gromacs/mdrun/mdmodules.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/mdtypes/pull_params.h"
63 #include "gromacs/options/options.h"
64 #include "gromacs/options/treesupport.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/selection/indexutil.h"
67 #include "gromacs/topology/block.h"
68 #include "gromacs/topology/ifunc.h"
69 #include "gromacs/topology/index.h"
70 #include "gromacs/topology/mtop_util.h"
71 #include "gromacs/topology/symtab.h"
72 #include "gromacs/topology/topology.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/exceptions.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/filestream.h"
77 #include "gromacs/utility/gmxassert.h"
78 #include "gromacs/utility/ikeyvaluetreeerror.h"
79 #include "gromacs/utility/keyvaluetree.h"
80 #include "gromacs/utility/keyvaluetreebuilder.h"
81 #include "gromacs/utility/keyvaluetreemdpwriter.h"
82 #include "gromacs/utility/keyvaluetreetransform.h"
83 #include "gromacs/utility/mdmodulenotification.h"
84 #include "gromacs/utility/smalloc.h"
85 #include "gromacs/utility/strconvert.h"
86 #include "gromacs/utility/stringcompare.h"
87 #include "gromacs/utility/stringutil.h"
88 #include "gromacs/utility/textwriter.h"
93 /* Resource parameters
94 * Do not change any of these until you read the instruction
95 * in readinp.h. Some cpp's do not take spaces after the backslash
96 * (like the c-shell), which will give you a very weird compiler
100 typedef struct t_inputrec_strings
102 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
103 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
104 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
105 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
106 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
108 char fep_lambda[efptNR][STRLEN];
109 char lambda_weights[STRLEN];
112 char anneal[STRLEN], anneal_npoints[STRLEN],
113 anneal_time[STRLEN], anneal_temp[STRLEN];
114 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
115 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
116 SAoff[STRLEN], SAsteps[STRLEN];
118 } gmx_inputrec_strings;
120 static gmx_inputrec_strings *is = nullptr;
122 void init_inputrec_strings()
126 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.");
131 void done_inputrec_strings()
139 egrptpALL, /* All particles have to be a member of a group. */
140 egrptpALL_GENREST, /* A rest group with name is generated for particles *
141 * that are not part of any group. */
142 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
143 * for the rest group. */
144 egrptpONE /* Merge all selected groups into one group, *
145 * make a rest group for the remaining particles. */
148 static const char *constraints[eshNR+1] = {
149 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
152 static const char *couple_lam[ecouplamNR+1] = {
153 "vdw-q", "vdw", "q", "none", nullptr
156 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
161 for (i = 0; i < ntemps; i++)
163 /* simple linear scaling -- allows more control */
164 if (simtemp->eSimTempScale == esimtempLINEAR)
166 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
168 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
170 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
172 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
174 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
179 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
180 gmx_fatal(FARGS, "%s", errorstr);
187 static void _low_check(bool b, const char *s, warninp_t wi)
191 warning_error(wi, s);
195 static void check_nst(const char *desc_nst, int nst,
196 const char *desc_p, int *p,
201 if (*p > 0 && *p % nst != 0)
203 /* Round up to the next multiple of nst */
204 *p = ((*p)/nst + 1)*nst;
205 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
206 desc_p, desc_nst, desc_p, *p);
211 static bool ir_NVE(const t_inputrec *ir)
213 return (EI_MD(ir->eI) && ir->etc == etcNO);
216 static int lcd(int n1, int n2)
221 for (i = 2; (i <= n1 && i <= n2); i++)
223 if (n1 % i == 0 && n2 % i == 0)
232 //! Convert legacy mdp entries to modern ones.
233 static void process_interaction_modifier(int *eintmod)
235 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
237 *eintmod = eintmodPOTSHIFT;
241 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
243 /* Check internal consistency.
244 * NOTE: index groups are not set here yet, don't check things
245 * like temperature coupling group options here, but in triple_check
248 /* Strange macro: first one fills the err_buf, and then one can check
249 * the condition, which will print the message and increase the error
252 #define CHECK(b) _low_check(b, err_buf, wi)
253 char err_buf[256], warn_buf[STRLEN];
256 t_lambda *fep = ir->fepvals;
257 t_expanded *expand = ir->expandedvals;
259 set_warning_line(wi, mdparin, -1);
261 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
263 sprintf(warn_buf, "%s electrostatics is no longer supported",
264 eel_names[eelRF_NEC_UNSUPPORTED]);
265 warning_error(wi, warn_buf);
268 /* BASIC CUT-OFF STUFF */
269 if (ir->rcoulomb < 0)
271 warning_error(wi, "rcoulomb should be >= 0");
275 warning_error(wi, "rvdw should be >= 0");
278 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
280 warning_error(wi, "rlist should be >= 0");
282 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.)");
283 CHECK(ir->nstlist < 0);
285 process_interaction_modifier(&ir->coulomb_modifier);
286 process_interaction_modifier(&ir->vdw_modifier);
288 if (ir->cutoff_scheme == ecutsGROUP)
291 "The group cutoff scheme has been removed since GROMACS 2020. "
292 "Please use the Verlet cutoff scheme.");
294 if (ir->cutoff_scheme == ecutsVERLET)
298 /* Normal Verlet type neighbor-list, currently only limited feature support */
299 if (inputrec2nboundeddim(ir) < 3)
301 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
304 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
305 if (ir->rcoulomb != ir->rvdw)
307 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
308 // for PME load balancing, we can support this exception.
309 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
310 ir->vdwtype == evdwCUT &&
311 ir->rcoulomb > ir->rvdw);
312 if (!bUsesPmeTwinRangeKernel)
314 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
318 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
320 if (ir->vdw_modifier == eintmodNONE ||
321 ir->vdw_modifier == eintmodPOTSHIFT)
323 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
325 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]);
326 warning_note(wi, warn_buf);
328 ir->vdwtype = evdwCUT;
332 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
333 warning_error(wi, warn_buf);
337 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
339 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
341 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
342 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
344 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
346 if (!(ir->coulomb_modifier == eintmodNONE ||
347 ir->coulomb_modifier == eintmodPOTSHIFT))
349 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
350 warning_error(wi, warn_buf);
353 if (EEL_USER(ir->coulombtype))
355 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme", eel_names[ir->coulombtype]);
356 warning_error(wi, warn_buf);
359 if (ir->nstlist <= 0)
361 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
364 if (ir->nstlist < 10)
366 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.");
369 rc_max = std::max(ir->rvdw, ir->rcoulomb);
371 if (ir->verletbuf_tol <= 0)
373 if (ir->verletbuf_tol == 0)
375 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
378 if (ir->rlist < rc_max)
380 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
383 if (ir->rlist == rc_max && ir->nstlist > 1)
385 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.");
390 if (ir->rlist > rc_max)
392 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.");
395 if (ir->nstlist == 1)
397 /* No buffer required */
402 if (EI_DYNAMICS(ir->eI))
404 if (inputrec2nboundeddim(ir) < 3)
406 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.");
408 /* Set rlist temporarily so we can continue processing */
413 /* Set the buffer to 5% of the cut-off */
414 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
420 /* GENERAL INTEGRATOR STUFF */
423 if (ir->etc != etcNO)
425 if (EI_RANDOM(ir->eI))
427 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]);
431 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
433 warning_note(wi, warn_buf);
437 if (ir->eI == eiVVAK)
439 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]);
440 warning_note(wi, warn_buf);
442 if (!EI_DYNAMICS(ir->eI))
444 if (ir->epc != epcNO)
446 sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
447 warning_note(wi, warn_buf);
451 if (EI_DYNAMICS(ir->eI))
453 if (ir->nstcalcenergy < 0)
455 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
456 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
458 /* nstcalcenergy larger than nstener does not make sense.
459 * We ideally want nstcalcenergy=nstener.
463 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
467 ir->nstcalcenergy = ir->nstenergy;
471 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
472 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
473 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
476 const char *nsten = "nstenergy";
477 const char *nstdh = "nstdhdl";
478 const char *min_name = nsten;
479 int min_nst = ir->nstenergy;
481 /* find the smallest of ( nstenergy, nstdhdl ) */
482 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
483 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
485 min_nst = ir->fepvals->nstdhdl;
488 /* If the user sets nstenergy small, we should respect that */
490 "Setting nstcalcenergy (%d) equal to %s (%d)",
491 ir->nstcalcenergy, min_name, min_nst);
492 warning_note(wi, warn_buf);
493 ir->nstcalcenergy = min_nst;
496 if (ir->epc != epcNO)
498 if (ir->nstpcouple < 0)
500 ir->nstpcouple = ir_optimal_nstpcouple(ir);
504 if (ir->nstcalcenergy > 0)
506 if (ir->efep != efepNO)
508 /* nstdhdl should be a multiple of nstcalcenergy */
509 check_nst("nstcalcenergy", ir->nstcalcenergy,
510 "nstdhdl", &ir->fepvals->nstdhdl, wi);
514 /* nstexpanded should be a multiple of nstcalcenergy */
515 check_nst("nstcalcenergy", ir->nstcalcenergy,
516 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
518 /* for storing exact averages nstenergy should be
519 * a multiple of nstcalcenergy
521 check_nst("nstcalcenergy", ir->nstcalcenergy,
522 "nstenergy", &ir->nstenergy, wi);
526 if (ir->nsteps == 0 && !ir->bContinuation)
528 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
532 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
533 ir->bContinuation && ir->ld_seed != -1)
535 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)");
541 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
542 CHECK(ir->ePBC != epbcXYZ);
543 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
544 CHECK(ir->ns_type != ensGRID);
545 sprintf(err_buf, "with TPI nstlist should be larger than zero");
546 CHECK(ir->nstlist <= 0);
547 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
548 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
549 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
550 CHECK(ir->cutoff_scheme == ecutsVERLET);
554 if ( (opts->nshake > 0) && (opts->bMorse) )
557 "Using morse bond-potentials while constraining bonds is useless");
558 warning(wi, warn_buf);
561 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
562 ir->bContinuation && ir->ld_seed != -1)
564 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)");
566 /* verify simulated tempering options */
570 bool bAllTempZero = TRUE;
571 for (i = 0; i < fep->n_lambda; i++)
573 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]);
574 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
575 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
577 bAllTempZero = FALSE;
580 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
581 CHECK(bAllTempZero == TRUE);
583 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
584 CHECK(ir->eI != eiVV);
586 /* check compatability of the temperature coupling with simulated tempering */
588 if (ir->etc == etcNOSEHOOVER)
590 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
591 warning_note(wi, warn_buf);
594 /* check that the temperatures make sense */
596 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);
597 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
599 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
600 CHECK(ir->simtempvals->simtemp_high <= 0);
602 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
603 CHECK(ir->simtempvals->simtemp_low <= 0);
606 /* verify free energy options */
608 if (ir->efep != efepNO)
611 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
613 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
615 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
616 static_cast<int>(fep->sc_r_power));
617 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
619 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
620 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
622 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
623 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
625 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
626 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
628 sprintf(err_buf, "Free-energy not implemented for Ewald");
629 CHECK(ir->coulombtype == eelEWALD);
631 /* check validty of lambda inputs */
632 if (fep->n_lambda == 0)
634 /* Clear output in case of no states:*/
635 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
636 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
640 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
641 CHECK((fep->init_fep_state >= fep->n_lambda));
644 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
645 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
647 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",
648 fep->init_lambda, fep->init_fep_state);
649 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
653 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
657 for (i = 0; i < efptNR; i++)
659 if (fep->separate_dvdl[i])
664 if (n_lambda_terms > 1)
666 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.");
667 warning(wi, warn_buf);
670 if (n_lambda_terms < 2 && fep->n_lambda > 0)
673 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
677 for (j = 0; j < efptNR; j++)
679 for (i = 0; i < fep->n_lambda; i++)
681 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]);
682 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
686 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
688 for (i = 0; i < fep->n_lambda; i++)
690 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],
691 fep->all_lambda[efptCOUL][i]);
692 CHECK((fep->sc_alpha > 0) &&
693 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
694 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
695 ((fep->all_lambda[efptVDW][i] > 0.0) &&
696 (fep->all_lambda[efptVDW][i] < 1.0))));
700 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
702 real sigma, lambda, r_sc;
705 /* Maximum estimate for A and B charges equal with lambda power 1 */
707 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);
708 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.",
710 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
711 warning_note(wi, warn_buf);
714 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
715 be treated differently, but that's the next step */
717 for (i = 0; i < efptNR; i++)
719 for (j = 0; j < fep->n_lambda; j++)
721 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
722 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
727 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
731 /* checking equilibration of weights inputs for validity */
733 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
734 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
735 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
737 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
738 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
739 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
741 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
742 expand->equil_steps, elmceq_names[elmceqSTEPS]);
743 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
745 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
746 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
747 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
749 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
750 expand->equil_ratio, elmceq_names[elmceqRATIO]);
751 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
753 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
754 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
755 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
757 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
758 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
759 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
761 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
762 expand->equil_steps, elmceq_names[elmceqSTEPS]);
763 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
765 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
766 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
767 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
769 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
770 expand->equil_ratio, elmceq_names[elmceqRATIO]);
771 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
773 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
774 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
775 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
777 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
778 CHECK((expand->lmc_repeats <= 0));
779 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
780 CHECK((expand->minvarmin <= 0));
781 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
782 CHECK((expand->c_range < 0));
783 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
784 fep->init_fep_state, expand->lmc_forced_nstart);
785 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
786 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
787 CHECK((expand->lmc_forced_nstart < 0));
788 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
789 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
791 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
792 CHECK((expand->init_wl_delta < 0));
793 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
794 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
795 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
796 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
798 /* if there is no temperature control, we need to specify an MC temperature */
799 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO) && (expand->mc_temp <= 0.0))
801 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set to a positive number");
802 warning_error(wi, err_buf);
804 if (expand->nstTij > 0)
806 sprintf(err_buf, "nstlog must be non-zero");
807 CHECK(ir->nstlog == 0);
808 // Avoid modulus by zero in the case that already triggered an error exit.
811 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
812 expand->nstTij, ir->nstlog);
813 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 (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2703 for (gmx::index 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 (gmx::index 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 (gmx::index 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 (gmx::index 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 (gmx::index 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,
3055 const gmx::MdModulesNotifier ¬ifier,
3059 t_blocka *defaultIndexGroups;
3063 char warnbuf[STRLEN], **gnames;
3067 int i, j, k, restnm;
3068 bool bExcl, bTable, bAnneal, bRest;
3069 char warn_buf[STRLEN];
3073 fprintf(stderr, "processing index file...\n");
3077 snew(defaultIndexGroups, 1);
3078 snew(defaultIndexGroups->index, 1);
3080 atoms_all = gmx_mtop_global_atoms(mtop);
3081 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3082 done_atom(&atoms_all);
3086 defaultIndexGroups = init_index(ndx, &gnames);
3089 SimulationGroups *groups = &mtop->groups;
3090 natoms = mtop->natoms;
3091 symtab = &mtop->symtab;
3093 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3095 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3097 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3098 restnm = groups->groupNames.size() - 1;
3099 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3100 srenew(gnames, defaultIndexGroups->nr+1);
3101 gnames[restnm] = *(groups->groupNames.back());
3103 set_warning_line(wi, mdparin, -1);
3105 auto temperatureCouplingTauValues = gmx::splitString(is->tau_t);
3106 auto temperatureCouplingReferenceValues = gmx::splitString(is->ref_t);
3107 auto temperatureCouplingGroupNames = gmx::splitString(is->tcgrps);
3108 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size() ||
3109 temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3111 gmx_fatal(FARGS, "Invalid T coupling input: %zu groups, %zu ref-t values and "
3113 temperatureCouplingGroupNames.size(),
3114 temperatureCouplingReferenceValues.size(),
3115 temperatureCouplingTauValues.size());
3118 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3119 do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
3120 SimulationAtomGroupType::TemperatureCoupling,
3121 restnm, useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3122 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3124 snew(ir->opts.nrdf, nr);
3125 snew(ir->opts.tau_t, nr);
3126 snew(ir->opts.ref_t, nr);
3127 if (ir->eI == eiBD && ir->bd_fric == 0)
3129 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3132 if (useReferenceTemperature)
3134 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3136 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3140 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3141 for (i = 0; (i < nr); i++)
3143 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3145 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3146 warning_error(wi, warn_buf);
3149 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3151 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.");
3154 if (ir->opts.tau_t[i] >= 0)
3156 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3159 if (ir->etc != etcNO && ir->nsttcouple == -1)
3161 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3166 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3168 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");
3170 if (ir->epc == epcMTTK)
3172 if (ir->etc != etcNOSEHOOVER)
3174 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3178 if (ir->nstpcouple != ir->nsttcouple)
3180 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3181 ir->nstpcouple = ir->nsttcouple = mincouple;
3182 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);
3183 warning_note(wi, warn_buf);
3188 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3189 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3191 if (ETC_ANDERSEN(ir->etc))
3193 if (ir->nsttcouple != 1)
3196 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");
3197 warning_note(wi, warn_buf);
3200 nstcmin = tcouple_min_integration_steps(ir->etc);
3203 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3205 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3206 ETCOUPLTYPE(ir->etc),
3208 ir->nsttcouple*ir->delta_t);
3209 warning(wi, warn_buf);
3212 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3213 for (i = 0; (i < nr); i++)
3215 if (ir->opts.ref_t[i] < 0)
3217 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3220 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3221 if we are in this conditional) if mc_temp is negative */
3222 if (ir->expandedvals->mc_temp < 0)
3224 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3228 /* Simulated annealing for each group. There are nr groups */
3229 auto simulatedAnnealingGroupNames = gmx::splitString(is->anneal);
3230 if (simulatedAnnealingGroupNames.size() == 1 &&
3231 gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3233 simulatedAnnealingGroupNames.resize(0);
3235 if (!simulatedAnnealingGroupNames.empty() &&
3236 gmx::ssize(simulatedAnnealingGroupNames) != nr)
3238 gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
3239 simulatedAnnealingGroupNames.size(), nr);
3243 snew(ir->opts.annealing, nr);
3244 snew(ir->opts.anneal_npoints, nr);
3245 snew(ir->opts.anneal_time, nr);
3246 snew(ir->opts.anneal_temp, nr);
3247 for (i = 0; i < nr; i++)
3249 ir->opts.annealing[i] = eannNO;
3250 ir->opts.anneal_npoints[i] = 0;
3251 ir->opts.anneal_time[i] = nullptr;
3252 ir->opts.anneal_temp[i] = nullptr;
3254 if (!simulatedAnnealingGroupNames.empty())
3257 for (i = 0; i < nr; i++)
3259 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3261 ir->opts.annealing[i] = eannNO;
3263 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3265 ir->opts.annealing[i] = eannSINGLE;
3268 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3270 ir->opts.annealing[i] = eannPERIODIC;
3276 /* Read the other fields too */
3277 auto simulatedAnnealingPoints = gmx::splitString(is->anneal_npoints);
3278 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3280 gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
3281 simulatedAnnealingPoints.size(), simulatedAnnealingGroupNames.size());
3283 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3284 size_t numSimulatedAnnealingFields = 0;
3285 for (i = 0; i < nr; i++)
3287 if (ir->opts.anneal_npoints[i] == 1)
3289 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3291 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3292 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3293 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3296 auto simulatedAnnealingTimes = gmx::splitString(is->anneal_time);
3298 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3300 gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
3301 simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
3303 auto simulatedAnnealingTemperatures = gmx::splitString(is->anneal_temp);
3304 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3306 gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
3307 simulatedAnnealingTemperatures.size(), numSimulatedAnnealingFields);
3310 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3311 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3312 convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
3313 convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp", allSimulatedAnnealingTemperatures.data());
3314 for (i = 0, k = 0; i < nr; i++)
3316 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3318 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3319 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3322 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3324 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3330 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3332 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3333 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3336 if (ir->opts.anneal_temp[i][j] < 0)
3338 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3343 /* Print out some summary information, to make sure we got it right */
3344 for (i = 0; i < nr; i++)
3346 if (ir->opts.annealing[i] != eannNO)
3348 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3349 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3350 *(groups->groupNames[j]), eann_names[ir->opts.annealing[i]],
3351 ir->opts.anneal_npoints[i]);
3352 fprintf(stderr, "Time (ps) Temperature (K)\n");
3353 /* All terms except the last one */
3354 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3356 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3359 /* Finally the last one */
3360 j = ir->opts.anneal_npoints[i]-1;
3361 if (ir->opts.annealing[i] == eannSINGLE)
3363 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3367 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3368 if (std::fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3370 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3381 make_pull_groups(ir->pull, is->pull_grp, defaultIndexGroups, gnames);
3383 make_pull_coords(ir->pull);
3388 make_rotation_groups(ir->rot, is->rot_grp, defaultIndexGroups, gnames);
3391 if (ir->eSwapCoords != eswapNO)
3393 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3396 /* Make indices for IMD session */
3399 make_IMD_group(ir->imd, is->imd_grp, defaultIndexGroups, gnames);
3402 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3403 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3404 notifier.notifier_.notify(defaultIndexGroupsAndNames);
3406 auto accelerations = gmx::splitString(is->acc);
3407 auto accelerationGroupNames = gmx::splitString(is->accgrps);
3408 if (accelerationGroupNames.size() * DIM != accelerations.size())
3410 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3411 accelerationGroupNames.size(), accelerations.size());
3413 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3414 SimulationAtomGroupType::Acceleration,
3415 restnm, egrptpALL_GENREST, bVerbose, wi);
3416 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3417 snew(ir->opts.acc, nr);
3418 ir->opts.ngacc = nr;
3420 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3422 auto freezeDims = gmx::splitString(is->frdim);
3423 auto freezeGroupNames = gmx::splitString(is->freeze);
3424 if (freezeDims.size() != DIM * freezeGroupNames.size())
3426 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3427 freezeGroupNames.size(), freezeDims.size());
3429 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3430 SimulationAtomGroupType::Freeze,
3431 restnm, egrptpALL_GENREST, bVerbose, wi);
3432 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3433 ir->opts.ngfrz = nr;
3434 snew(ir->opts.nFreeze, nr);
3435 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3437 for (j = 0; (j < DIM); j++, k++)
3439 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3440 if (!ir->opts.nFreeze[i][j])
3442 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3444 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3445 "(not %s)", freezeDims[k].c_str());
3446 warning(wi, warn_buf);
3451 for (; (i < nr); i++)
3453 for (j = 0; (j < DIM); j++)
3455 ir->opts.nFreeze[i][j] = 0;
3459 auto energyGroupNames = gmx::splitString(is->energy);
3460 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3461 SimulationAtomGroupType::EnergyOutput,
3462 restnm, egrptpALL_GENREST, bVerbose, wi);
3463 add_wall_energrps(groups, ir->nwall, symtab);
3464 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3465 auto vcmGroupNames = gmx::splitString(is->vcm);
3467 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3468 SimulationAtomGroupType::MassCenterVelocityRemoval,
3469 restnm, vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3472 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3473 "This may lead to artifacts.\n"
3474 "In most cases one should use one group for the whole system.");
3477 /* Now we have filled the freeze struct, so we can calculate NRDF */
3478 calc_nrdf(mtop, ir, gnames);
3480 auto user1GroupNames = gmx::splitString(is->user1);
3481 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3482 SimulationAtomGroupType::User1,
3483 restnm, egrptpALL_GENREST, bVerbose, wi);
3484 auto user2GroupNames = gmx::splitString(is->user2);
3485 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3486 SimulationAtomGroupType::User2,
3487 restnm, egrptpALL_GENREST, bVerbose, wi);
3488 auto compressedXGroupNames = gmx::splitString(is->x_compressed_groups);
3489 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3490 SimulationAtomGroupType::CompressedPositionOutput,
3491 restnm, egrptpONE, bVerbose, wi);
3492 auto orirefFitGroupNames = gmx::splitString(is->orirefitgrp);
3493 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3494 SimulationAtomGroupType::OrientationRestraintsFit,
3495 restnm, egrptpALL_GENREST, bVerbose, wi);
3497 /* QMMM input processing */
3498 auto qmGroupNames = gmx::splitString(is->QMMM);
3499 auto qmMethods = gmx::splitString(is->QMmethod);
3500 auto qmBasisSets = gmx::splitString(is->QMbasis);
3501 if (ir->eI != eiMimic)
3503 if (qmMethods.size() != qmGroupNames.size() ||
3504 qmBasisSets.size() != qmGroupNames.size())
3506 gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets"
3507 " and %zu methods\n", qmGroupNames.size(),
3508 qmBasisSets.size(), qmMethods.size());
3510 /* group rest, if any, is always MM! */
3511 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3512 SimulationAtomGroupType::QuantumMechanics,
3513 restnm, egrptpALL_GENREST, bVerbose, wi);
3514 nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
3515 ir->opts.ngQM = qmGroupNames.size();
3516 snew(ir->opts.QMmethod, nr);
3517 snew(ir->opts.QMbasis, nr);
3518 for (i = 0; i < nr; i++)
3520 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3521 * converted to the corresponding enum in names.c
3523 ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(),
3526 ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(),
3531 auto qmMultiplicities = gmx::splitString(is->QMmult);
3532 auto qmCharges = gmx::splitString(is->QMcharge);
3533 auto qmbSH = gmx::splitString(is->bSH);
3534 snew(ir->opts.QMmult, nr);
3535 snew(ir->opts.QMcharge, nr);
3536 snew(ir->opts.bSH, nr);
3537 convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult);
3538 convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge);
3539 convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH);
3541 auto CASelectrons = gmx::splitString(is->CASelectrons);
3542 auto CASorbitals = gmx::splitString(is->CASorbitals);
3543 snew(ir->opts.CASelectrons, nr);
3544 snew(ir->opts.CASorbitals, nr);
3545 convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons);
3546 convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals);
3548 auto SAon = gmx::splitString(is->SAon);
3549 auto SAoff = gmx::splitString(is->SAoff);
3550 auto SAsteps = gmx::splitString(is->SAsteps);
3551 snew(ir->opts.SAon, nr);
3552 snew(ir->opts.SAoff, nr);
3553 snew(ir->opts.SAsteps, nr);
3554 convertInts(wi, SAon, "SAon", ir->opts.SAon);
3555 convertInts(wi, SAoff, "SAoff", ir->opts.SAoff);
3556 convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps);
3561 if (qmGroupNames.size() > 1)
3563 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3565 /* group rest, if any, is always MM! */
3566 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3567 SimulationAtomGroupType::QuantumMechanics,
3568 restnm, egrptpALL_GENREST, bVerbose, wi);
3570 ir->opts.ngQM = qmGroupNames.size();
3573 /* end of QMMM input */
3577 for (auto group : gmx::keysOf(groups->groups))
3579 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3580 for (const auto &entry : groups->groups[group])
3582 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3584 fprintf(stderr, "\n");
3588 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3589 snew(ir->opts.egp_flags, nr*nr);
3591 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3592 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3594 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3596 if (bExcl && EEL_FULL(ir->coulombtype))
3598 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3601 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3602 if (bTable && !(ir->vdwtype == evdwUSER) &&
3603 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3604 !(ir->coulombtype == eelPMEUSERSWITCH))
3606 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3609 /* final check before going out of scope if simulated tempering variables
3610 * need to be set to default values.
3612 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3614 ir->expandedvals->nstexpanded = 2*static_cast<int>(ir->opts.tau_t[0]/ir->delta_t);
3615 warning(wi, gmx::formatString("the value for nstexpanded was not specified for "
3616 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3617 "by default, but it is recommended to set it to an explicit value!",
3618 ir->expandedvals->nstexpanded));
3620 for (i = 0; (i < defaultIndexGroups->nr); i++)
3625 done_blocka(defaultIndexGroups);
3626 sfree(defaultIndexGroups);
3632 static void check_disre(const gmx_mtop_t *mtop)
3634 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3636 const gmx_ffparams_t &ffparams = mtop->ffparams;
3639 for (int i = 0; i < ffparams.numTypes(); i++)
3641 int ftype = ffparams.functype[i];
3642 if (ftype == F_DISRES)
3644 int label = ffparams.iparams[i].disres.label;
3645 if (label == old_label)
3647 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3655 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3656 "probably the parameters for multiple pairs in one restraint "
3657 "are not identical\n", ndouble);
3662 static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3667 gmx_mtop_ilistloop_t iloop;
3676 for (d = 0; d < DIM; d++)
3678 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3680 /* Check for freeze groups */
3681 for (g = 0; g < ir->opts.ngfrz; g++)
3683 for (d = 0; d < DIM; d++)
3685 if (ir->opts.nFreeze[g][d] != 0)
3693 /* Check for position restraints */
3694 iloop = gmx_mtop_ilistloop_init(sys);
3695 while (const InteractionLists *ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
3698 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3700 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
3702 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
3703 for (d = 0; d < DIM; d++)
3705 if (pr->posres.fcA[d] != 0)
3711 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
3713 /* Check for flat-bottom posres */
3714 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
3715 if (pr->fbposres.k != 0)
3717 switch (pr->fbposres.geom)
3719 case efbposresSPHERE:
3720 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3722 case efbposresCYLINDERX:
3723 AbsRef[YY] = AbsRef[ZZ] = 1;
3725 case efbposresCYLINDERY:
3726 AbsRef[XX] = AbsRef[ZZ] = 1;
3728 case efbposresCYLINDER:
3729 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3730 case efbposresCYLINDERZ:
3731 AbsRef[XX] = AbsRef[YY] = 1;
3733 case efbposresX: /* d=XX */
3734 case efbposresY: /* d=YY */
3735 case efbposresZ: /* d=ZZ */
3736 d = pr->fbposres.geom - efbposresX;
3740 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3741 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3749 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3753 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3754 bool *bC6ParametersWorkWithGeometricRules,
3755 bool *bC6ParametersWorkWithLBRules,
3756 bool *bLBRulesPossible)
3758 int ntypes, tpi, tpj;
3761 double c6i, c6j, c12i, c12j;
3762 double c6, c6_geometric, c6_LB;
3763 double sigmai, sigmaj, epsi, epsj;
3764 bool bCanDoLBRules, bCanDoGeometricRules;
3767 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3768 * force-field floating point parameters.
3771 ptr = getenv("GMX_LJCOMB_TOL");
3775 double gmx_unused canary;
3777 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3779 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
3784 *bC6ParametersWorkWithLBRules = TRUE;
3785 *bC6ParametersWorkWithGeometricRules = TRUE;
3786 bCanDoLBRules = TRUE;
3787 ntypes = mtop->ffparams.atnr;
3788 snew(typecount, ntypes);
3789 gmx_mtop_count_atomtypes(mtop, state, typecount);
3790 *bLBRulesPossible = TRUE;
3791 for (tpi = 0; tpi < ntypes; ++tpi)
3793 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3794 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3795 for (tpj = tpi; tpj < ntypes; ++tpj)
3797 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3798 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3799 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3800 c6_geometric = std::sqrt(c6i * c6j);
3801 if (!gmx_numzero(c6_geometric))
3803 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3805 sigmai = gmx::sixthroot(c12i / c6i);
3806 sigmaj = gmx::sixthroot(c12j / c6j);
3807 epsi = c6i * c6i /(4.0 * c12i);
3808 epsj = c6j * c6j /(4.0 * c12j);
3809 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
3813 *bLBRulesPossible = FALSE;
3814 c6_LB = c6_geometric;
3816 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3821 *bC6ParametersWorkWithLBRules = FALSE;
3824 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3826 if (!bCanDoGeometricRules)
3828 *bC6ParametersWorkWithGeometricRules = FALSE;
3836 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3839 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3841 check_combination_rule_differences(mtop, 0,
3842 &bC6ParametersWorkWithGeometricRules,
3843 &bC6ParametersWorkWithLBRules,
3845 if (ir->ljpme_combination_rule == eljpmeLB)
3847 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
3849 warning(wi, "You are using arithmetic-geometric combination rules "
3850 "in LJ-PME, but your non-bonded C6 parameters do not "
3851 "follow these rules.");
3856 if (!bC6ParametersWorkWithGeometricRules)
3858 if (ir->eDispCorr != edispcNO)
3860 warning_note(wi, "You are using geometric combination rules in "
3861 "LJ-PME, but your non-bonded C6 parameters do "
3862 "not follow these rules. "
3863 "This will introduce very small errors in the forces and energies in "
3864 "your simulations. Dispersion correction will correct total energy "
3865 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3869 warning_note(wi, "You are using geometric combination rules in "
3870 "LJ-PME, but your non-bonded C6 parameters do "
3871 "not follow these rules. "
3872 "This will introduce very small errors in the forces and energies in "
3873 "your simulations. If your system is homogeneous, consider using dispersion correction "
3874 "for the total energy and pressure.");
3880 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3883 char err_buf[STRLEN];
3888 gmx_mtop_atomloop_block_t aloopb;
3890 char warn_buf[STRLEN];
3892 set_warning_line(wi, mdparin, -1);
3894 if (ir->cutoff_scheme == ecutsVERLET &&
3895 ir->verletbuf_tol > 0 &&
3897 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
3898 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
3900 /* Check if a too small Verlet buffer might potentially
3901 * cause more drift than the thermostat can couple off.
3903 /* Temperature error fraction for warning and suggestion */
3904 const real T_error_warn = 0.002;
3905 const real T_error_suggest = 0.001;
3906 /* For safety: 2 DOF per atom (typical with constraints) */
3907 const real nrdf_at = 2;
3908 real T, tau, max_T_error;
3913 for (i = 0; i < ir->opts.ngtc; i++)
3915 T = std::max(T, ir->opts.ref_t[i]);
3916 tau = std::max(tau, ir->opts.tau_t[i]);
3920 /* This is a worst case estimate of the temperature error,
3921 * assuming perfect buffer estimation and no cancelation
3922 * of errors. The factor 0.5 is because energy distributes
3923 * equally over Ekin and Epot.
3925 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
3926 if (max_T_error > T_error_warn)
3928 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.",
3929 ir->verletbuf_tol, T, tau,
3931 100*T_error_suggest,
3932 ir->verletbuf_tol*T_error_suggest/max_T_error);
3933 warning(wi, warn_buf);
3938 if (ETC_ANDERSEN(ir->etc))
3942 for (i = 0; i < ir->opts.ngtc; i++)
3944 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
3945 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
3946 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
3947 i, ir->opts.tau_t[i]);
3948 CHECK(ir->opts.tau_t[i] < 0);
3951 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
3953 for (i = 0; i < ir->opts.ngtc; i++)
3955 int nsteps = gmx::roundToInt(ir->opts.tau_t[i]/ir->delta_t);
3956 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);
3957 CHECK(nsteps % ir->nstcomm != 0);
3962 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3963 ir->comm_mode == ecmNO &&
3964 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
3965 !ETC_ANDERSEN(ir->etc))
3967 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");
3970 if (ir->epc == epcPARRINELLORAHMAN &&
3971 ir->etc == etcNOSEHOOVER)
3974 for (int g = 0; g < ir->opts.ngtc; g++)
3976 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
3978 if (ir->tau_p < 1.9*tau_t_max)
3980 std::string message =
3981 gmx::formatString("With %s T-coupling and %s p-coupling, "
3982 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
3983 etcoupl_names[ir->etc],
3984 epcoupl_names[ir->epc],
3986 "tau-t", tau_t_max);
3987 warning(wi, message.c_str());
3991 /* Check for pressure coupling with absolute position restraints */
3992 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3994 absolute_reference(ir, sys, TRUE, AbsRef);
3996 for (m = 0; m < DIM; m++)
3998 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4000 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4008 aloopb = gmx_mtop_atomloop_block_init(sys);
4010 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4012 if (atom->q != 0 || atom->qB != 0)
4020 if (EEL_FULL(ir->coulombtype))
4023 "You are using full electrostatics treatment %s for a system without charges.\n"
4024 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4025 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4026 warning(wi, err_buf);
4031 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4034 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4035 "You might want to consider using %s electrostatics.\n",
4037 warning_note(wi, err_buf);
4041 /* Check if combination rules used in LJ-PME are the same as in the force field */
4042 if (EVDW_PME(ir->vdwtype))
4044 check_combination_rules(ir, sys, wi);
4047 /* Generalized reaction field */
4048 if (ir->coulombtype == eelGRF_NOTUSED)
4050 warning_error(wi, "Generalized reaction-field electrostatics is no longer supported. "
4051 "You can use normal reaction-field instead and compute the reaction-field constant by hand.");
4055 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4057 for (m = 0; (m < DIM); m++)
4059 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4068 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4069 for (const AtomProxy atomP : AtomRange(*sys))
4071 const t_atom &local = atomP.atom();
4072 int i = atomP.globalAtomNumber();
4073 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4076 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4078 for (m = 0; (m < DIM); m++)
4080 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4084 for (m = 0; (m < DIM); m++)
4086 if (fabs(acc[m]) > 1e-6)
4088 const char *dim[DIM] = { "X", "Y", "Z" };
4090 "Net Acceleration in %s direction, will %s be corrected\n",
4091 dim[m], ir->nstcomm != 0 ? "" : "not");
4092 if (ir->nstcomm != 0 && m < ndof_com(ir))
4095 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4097 ir->opts.acc[i][m] -= acc[m];
4105 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4106 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4108 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4116 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4118 if (ir->pull->coord[i].group[0] == 0 ||
4119 ir->pull->coord[i].group[1] == 0)
4121 absolute_reference(ir, sys, FALSE, AbsRef);
4122 for (m = 0; m < DIM; m++)
4124 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4126 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.");
4134 for (i = 0; i < 3; i++)
4136 for (m = 0; m <= i; m++)
4138 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4139 ir->deform[i][m] != 0)
4141 for (c = 0; c < ir->pull->ncoord; c++)
4143 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4144 ir->pull->coord[c].vec[m] != 0)
4146 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4157 void double_check(t_inputrec *ir, matrix box,
4158 bool bHasNormalConstraints,
4159 bool bHasAnyConstraints,
4163 char warn_buf[STRLEN];
4166 ptr = check_box(ir->ePBC, box);
4169 warning_error(wi, ptr);
4172 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4174 if (ir->shake_tol <= 0.0)
4176 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4178 warning_error(wi, warn_buf);
4182 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4184 /* If we have Lincs constraints: */
4185 if (ir->eI == eiMD && ir->etc == etcNO &&
4186 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4188 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4189 warning_note(wi, warn_buf);
4192 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4194 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4195 warning_note(wi, warn_buf);
4197 if (ir->epc == epcMTTK)
4199 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4203 if (bHasAnyConstraints && ir->epc == epcMTTK)
4205 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4208 if (ir->LincsWarnAngle > 90.0)
4210 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4211 warning(wi, warn_buf);
4212 ir->LincsWarnAngle = 90.0;
4215 if (ir->ePBC != epbcNONE)
4217 if (ir->nstlist == 0)
4219 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4221 if (ir->ns_type == ensGRID)
4223 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4225 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");
4226 warning_error(wi, warn_buf);
4231 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4232 if (2*ir->rlist >= min_size)
4234 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4235 warning_error(wi, warn_buf);
4238 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4245 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4249 real rvdw1, rvdw2, rcoul1, rcoul2;
4250 char warn_buf[STRLEN];
4252 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4256 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4261 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4267 if (rvdw1 + rvdw2 > ir->rlist ||
4268 rcoul1 + rcoul2 > ir->rlist)
4271 "The sum of the two largest charge group radii (%f) "
4272 "is larger than rlist (%f)\n",
4273 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4274 warning(wi, warn_buf);
4278 /* Here we do not use the zero at cut-off macro,
4279 * since user defined interactions might purposely
4280 * not be zero at the cut-off.
4282 if (ir_vdw_is_zero_at_cutoff(ir) &&
4283 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4285 sprintf(warn_buf, "The sum of the two largest charge group "
4286 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4287 "With exact cut-offs, better performance can be "
4288 "obtained with cutoff-scheme = %s, because it "
4289 "does not use charge groups at all.",
4291 ir->rlist, ir->rvdw,
4292 ecutscheme_names[ecutsVERLET]);
4295 warning(wi, warn_buf);
4299 warning_note(wi, warn_buf);
4302 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4303 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4305 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4306 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4308 ir->rlist, ir->rcoulomb,
4309 ecutscheme_names[ecutsVERLET]);
4312 warning(wi, warn_buf);
4316 warning_note(wi, warn_buf);