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/network.h"
53 #include "gromacs/gmxpreprocess/toputil.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdlib/calc_verletbuf.h"
58 #include "gromacs/mdrun/mdmodules.h"
59 #include "gromacs/mdtypes/inputrec.h"
60 #include "gromacs/mdtypes/md_enums.h"
61 #include "gromacs/mdtypes/pull_params.h"
62 #include "gromacs/options/options.h"
63 #include "gromacs/options/treesupport.h"
64 #include "gromacs/pbcutil/pbc.h"
65 #include "gromacs/selection/indexutil.h"
66 #include "gromacs/topology/block.h"
67 #include "gromacs/topology/ifunc.h"
68 #include "gromacs/topology/index.h"
69 #include "gromacs/topology/mtop_util.h"
70 #include "gromacs/topology/symtab.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/utility/cstringutil.h"
73 #include "gromacs/utility/exceptions.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/filestream.h"
76 #include "gromacs/utility/gmxassert.h"
77 #include "gromacs/utility/ikeyvaluetreeerror.h"
78 #include "gromacs/utility/keyvaluetree.h"
79 #include "gromacs/utility/keyvaluetreebuilder.h"
80 #include "gromacs/utility/keyvaluetreemdpwriter.h"
81 #include "gromacs/utility/keyvaluetreetransform.h"
82 #include "gromacs/utility/mdmodulenotification.h"
83 #include "gromacs/utility/smalloc.h"
84 #include "gromacs/utility/strconvert.h"
85 #include "gromacs/utility/stringcompare.h"
86 #include "gromacs/utility/stringutil.h"
87 #include "gromacs/utility/textwriter.h"
92 /* Resource parameters
93 * Do not change any of these until you read the instruction
94 * in readinp.h. Some cpp's do not take spaces after the backslash
95 * (like the c-shell), which will give you a very weird compiler
99 typedef struct t_inputrec_strings
101 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
102 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
103 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
104 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
105 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
107 char fep_lambda[efptNR][STRLEN];
108 char lambda_weights[STRLEN];
111 char anneal[STRLEN], anneal_npoints[STRLEN],
112 anneal_time[STRLEN], anneal_temp[STRLEN];
113 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
114 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
115 SAoff[STRLEN], SAsteps[STRLEN];
117 } gmx_inputrec_strings;
119 static gmx_inputrec_strings *is = nullptr;
121 void init_inputrec_strings()
125 gmx_incons("Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time.");
130 void done_inputrec_strings()
138 egrptpALL, /* All particles have to be a member of a group. */
139 egrptpALL_GENREST, /* A rest group with name is generated for particles *
140 * that are not part of any group. */
141 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
142 * for the rest group. */
143 egrptpONE /* Merge all selected groups into one group, *
144 * make a rest group for the remaining particles. */
147 static const char *constraints[eshNR+1] = {
148 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
151 static const char *couple_lam[ecouplamNR+1] = {
152 "vdw-q", "vdw", "q", "none", nullptr
155 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
160 for (i = 0; i < ntemps; i++)
162 /* simple linear scaling -- allows more control */
163 if (simtemp->eSimTempScale == esimtempLINEAR)
165 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
167 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
169 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
171 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
173 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
178 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
179 gmx_fatal(FARGS, "%s", errorstr);
186 static void _low_check(bool b, const char *s, warninp_t wi)
190 warning_error(wi, s);
194 static void check_nst(const char *desc_nst, int nst,
195 const char *desc_p, int *p,
200 if (*p > 0 && *p % nst != 0)
202 /* Round up to the next multiple of nst */
203 *p = ((*p)/nst + 1)*nst;
204 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
205 desc_p, desc_nst, desc_p, *p);
210 static int lcd(int n1, int n2)
215 for (i = 2; (i <= n1 && i <= n2); i++)
217 if (n1 % i == 0 && n2 % i == 0)
226 //! Convert legacy mdp entries to modern ones.
227 static void process_interaction_modifier(int *eintmod)
229 if (*eintmod == eintmodPOTSHIFT_VERLET_UNSUPPORTED)
231 *eintmod = eintmodPOTSHIFT;
235 void check_ir(const char *mdparin, const gmx::MdModulesNotifier &mdModulesNotifier,
236 t_inputrec *ir, t_gromppopts *opts, warninp_t wi)
237 /* Check internal consistency.
238 * NOTE: index groups are not set here yet, don't check things
239 * like temperature coupling group options here, but in triple_check
242 /* Strange macro: first one fills the err_buf, and then one can check
243 * the condition, which will print the message and increase the error
246 #define CHECK(b) _low_check(b, err_buf, wi)
247 char err_buf[256], warn_buf[STRLEN];
250 t_lambda *fep = ir->fepvals;
251 t_expanded *expand = ir->expandedvals;
253 set_warning_line(wi, mdparin, -1);
255 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
257 sprintf(warn_buf, "%s electrostatics is no longer supported",
258 eel_names[eelRF_NEC_UNSUPPORTED]);
259 warning_error(wi, warn_buf);
262 /* BASIC CUT-OFF STUFF */
263 if (ir->rcoulomb < 0)
265 warning_error(wi, "rcoulomb should be >= 0");
269 warning_error(wi, "rvdw should be >= 0");
272 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
274 warning_error(wi, "rlist should be >= 0");
276 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.)");
277 CHECK(ir->nstlist < 0);
279 process_interaction_modifier(&ir->coulomb_modifier);
280 process_interaction_modifier(&ir->vdw_modifier);
282 if (ir->cutoff_scheme == ecutsGROUP)
285 "The group cutoff scheme has been removed since GROMACS 2020. "
286 "Please use the Verlet cutoff scheme.");
288 if (ir->cutoff_scheme == ecutsVERLET)
292 /* Normal Verlet type neighbor-list, currently only limited feature support */
293 if (inputrec2nboundeddim(ir) < 3)
295 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
298 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
299 if (ir->rcoulomb != ir->rvdw)
301 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
302 // for PME load balancing, we can support this exception.
303 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
304 ir->vdwtype == evdwCUT &&
305 ir->rcoulomb > ir->rvdw);
306 if (!bUsesPmeTwinRangeKernel)
308 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
312 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
314 if (ir->vdw_modifier == eintmodNONE ||
315 ir->vdw_modifier == eintmodPOTSHIFT)
317 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
319 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]);
320 warning_note(wi, warn_buf);
322 ir->vdwtype = evdwCUT;
326 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
327 warning_error(wi, warn_buf);
331 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
333 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
335 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
336 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
338 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
340 if (!(ir->coulomb_modifier == eintmodNONE ||
341 ir->coulomb_modifier == eintmodPOTSHIFT))
343 sprintf(warn_buf, "coulomb_modifier=%s is not supported", eintmod_names[ir->coulomb_modifier]);
344 warning_error(wi, warn_buf);
347 if (EEL_USER(ir->coulombtype))
349 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme", eel_names[ir->coulombtype]);
350 warning_error(wi, warn_buf);
353 if (ir->nstlist <= 0)
355 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
358 if (ir->nstlist < 10)
360 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.");
363 rc_max = std::max(ir->rvdw, ir->rcoulomb);
367 /* With TPI we set the pairlist cut-off later using the radius of the insterted molecule */
368 ir->verletbuf_tol = 0;
371 else 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);
525 // Inquire all MdModules, if their parameters match with the energy
526 // calculation frequency
527 gmx::EnergyCalculationFrequencyErrors energyCalculationFrequencyErrors(ir->nstcalcenergy);
528 mdModulesNotifier.notifier_.notify(&energyCalculationFrequencyErrors);
530 // Emit all errors from the energy calculation frequency checks
531 for (const std::string &energyFrequencyErrorMessage : energyCalculationFrequencyErrors.errorMessages())
533 warning_error(wi, energyFrequencyErrorMessage);
537 if (ir->nsteps == 0 && !ir->bContinuation)
539 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
543 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
544 ir->bContinuation && ir->ld_seed != -1)
546 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)");
552 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
553 CHECK(ir->ePBC != epbcXYZ);
554 sprintf(err_buf, "with TPI nstlist should be larger than zero");
555 CHECK(ir->nstlist <= 0);
556 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
557 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
561 if ( (opts->nshake > 0) && (opts->bMorse) )
564 "Using morse bond-potentials while constraining bonds is useless");
565 warning(wi, warn_buf);
568 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
569 ir->bContinuation && ir->ld_seed != -1)
571 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)");
573 /* verify simulated tempering options */
577 bool bAllTempZero = TRUE;
578 for (i = 0; i < fep->n_lambda; i++)
580 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]);
581 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
582 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
584 bAllTempZero = FALSE;
587 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
588 CHECK(bAllTempZero == TRUE);
590 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
591 CHECK(ir->eI != eiVV);
593 /* check compatability of the temperature coupling with simulated tempering */
595 if (ir->etc == etcNOSEHOOVER)
597 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
598 warning_note(wi, warn_buf);
601 /* check that the temperatures make sense */
603 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);
604 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
606 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
607 CHECK(ir->simtempvals->simtemp_high <= 0);
609 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
610 CHECK(ir->simtempvals->simtemp_low <= 0);
613 /* verify free energy options */
615 if (ir->efep != efepNO)
618 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
620 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
622 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
623 static_cast<int>(fep->sc_r_power));
624 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
626 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
627 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
629 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
630 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
632 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
633 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
635 sprintf(err_buf, "Free-energy not implemented for Ewald");
636 CHECK(ir->coulombtype == eelEWALD);
638 /* check validty of lambda inputs */
639 if (fep->n_lambda == 0)
641 /* Clear output in case of no states:*/
642 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
643 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
647 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
648 CHECK((fep->init_fep_state >= fep->n_lambda));
651 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
652 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
654 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",
655 fep->init_lambda, fep->init_fep_state);
656 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
660 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
664 for (i = 0; i < efptNR; i++)
666 if (fep->separate_dvdl[i])
671 if (n_lambda_terms > 1)
673 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.");
674 warning(wi, warn_buf);
677 if (n_lambda_terms < 2 && fep->n_lambda > 0)
680 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
684 for (j = 0; j < efptNR; j++)
686 for (i = 0; i < fep->n_lambda; i++)
688 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]);
689 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
693 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
695 for (i = 0; i < fep->n_lambda; i++)
697 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],
698 fep->all_lambda[efptCOUL][i]);
699 CHECK((fep->sc_alpha > 0) &&
700 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
701 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
702 ((fep->all_lambda[efptVDW][i] > 0.0) &&
703 (fep->all_lambda[efptVDW][i] < 1.0))));
707 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
709 real sigma, lambda, r_sc;
712 /* Maximum estimate for A and B charges equal with lambda power 1 */
714 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);
715 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.",
717 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
718 warning_note(wi, warn_buf);
721 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
722 be treated differently, but that's the next step */
724 for (i = 0; i < efptNR; i++)
726 for (j = 0; j < fep->n_lambda; j++)
728 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
729 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
734 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
738 /* checking equilibration of weights inputs for validity */
740 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
741 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
742 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
744 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
745 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
746 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
748 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
749 expand->equil_steps, elmceq_names[elmceqSTEPS]);
750 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
752 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
753 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
754 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
756 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
757 expand->equil_ratio, elmceq_names[elmceqRATIO]);
758 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
760 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
761 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
762 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
764 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
765 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
766 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
768 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
769 expand->equil_steps, elmceq_names[elmceqSTEPS]);
770 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
772 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
773 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
774 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
776 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
777 expand->equil_ratio, elmceq_names[elmceqRATIO]);
778 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
780 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
781 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
782 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
784 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
785 CHECK((expand->lmc_repeats <= 0));
786 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
787 CHECK((expand->minvarmin <= 0));
788 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
789 CHECK((expand->c_range < 0));
790 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
791 fep->init_fep_state, expand->lmc_forced_nstart);
792 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
793 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
794 CHECK((expand->lmc_forced_nstart < 0));
795 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
796 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
798 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
799 CHECK((expand->init_wl_delta < 0));
800 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
801 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
802 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
803 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
805 /* if there is no temperature control, we need to specify an MC temperature */
806 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO) && (expand->mc_temp <= 0.0))
808 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set to a positive number");
809 warning_error(wi, err_buf);
811 if (expand->nstTij > 0)
813 sprintf(err_buf, "nstlog must be non-zero");
814 CHECK(ir->nstlog == 0);
815 // Avoid modulus by zero in the case that already triggered an error exit.
818 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
819 expand->nstTij, ir->nstlog);
820 CHECK((expand->nstTij % ir->nstlog) != 0);
826 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
827 CHECK(ir->nwall && ir->ePBC != epbcXY);
830 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
832 if (ir->ePBC == epbcNONE)
834 if (ir->epc != epcNO)
836 warning(wi, "Turning off pressure coupling for vacuum system");
842 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
843 epbc_names[ir->ePBC]);
844 CHECK(ir->epc != epcNO);
846 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
847 CHECK(EEL_FULL(ir->coulombtype));
849 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
850 epbc_names[ir->ePBC]);
851 CHECK(ir->eDispCorr != edispcNO);
854 if (ir->rlist == 0.0)
856 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
857 "with coulombtype = %s or coulombtype = %s\n"
858 "without periodic boundary conditions (pbc = %s) and\n"
859 "rcoulomb and rvdw set to zero",
860 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
861 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
862 (ir->ePBC != epbcNONE) ||
863 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
867 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
872 if (ir->nstcomm == 0)
874 // TODO Change this behaviour. There should be exactly one way
875 // to turn off an algorithm.
876 ir->comm_mode = ecmNO;
878 if (ir->comm_mode != ecmNO)
882 // TODO Such input was once valid. Now that we've been
883 // helpful for a few years, we should reject such input,
884 // lest we have to support every historical decision
886 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");
887 ir->nstcomm = abs(ir->nstcomm);
890 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
892 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
893 ir->nstcomm = ir->nstcalcenergy;
896 if (ir->comm_mode == ecmANGULAR)
898 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
899 CHECK(ir->bPeriodicMols);
900 if (ir->ePBC != epbcNONE)
902 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.");
907 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
909 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]);
910 warning_note(wi, warn_buf);
913 /* TEMPERATURE COUPLING */
914 if (ir->etc == etcYES)
916 ir->etc = etcBERENDSEN;
917 warning_note(wi, "Old option for temperature coupling given: "
918 "changing \"yes\" to \"Berendsen\"\n");
921 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
923 if (ir->opts.nhchainlength < 1)
925 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
926 ir->opts.nhchainlength = 1;
927 warning(wi, warn_buf);
930 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
932 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
933 ir->opts.nhchainlength = 1;
938 ir->opts.nhchainlength = 0;
941 if (ir->eI == eiVVAK)
943 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
945 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
948 if (ETC_ANDERSEN(ir->etc))
950 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
951 CHECK(!(EI_VV(ir->eI)));
953 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
955 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]);
956 warning_note(wi, warn_buf);
959 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]);
960 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
963 if (ir->etc == etcBERENDSEN)
965 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
966 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
967 warning_note(wi, warn_buf);
970 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
971 && ir->epc == epcBERENDSEN)
973 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
974 "true ensemble for the thermostat");
975 warning(wi, warn_buf);
978 /* PRESSURE COUPLING */
979 if (ir->epc == epcISOTROPIC)
981 ir->epc = epcBERENDSEN;
982 warning_note(wi, "Old option for pressure coupling given: "
983 "changing \"Isotropic\" to \"Berendsen\"\n");
986 if (ir->epc != epcNO)
988 dt_pcoupl = ir->nstpcouple*ir->delta_t;
990 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
991 CHECK(ir->tau_p <= 0);
993 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
995 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
996 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
997 warning(wi, warn_buf);
1000 sprintf(err_buf, "compressibility must be > 0 when using pressure"
1001 " coupling %s\n", EPCOUPLTYPE(ir->epc));
1002 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
1003 ir->compress[ZZ][ZZ] < 0 ||
1004 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
1005 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1007 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1010 "You are generating velocities so I am assuming you "
1011 "are equilibrating a system. You are using "
1012 "%s pressure coupling, but this can be "
1013 "unstable for equilibration. If your system crashes, try "
1014 "equilibrating first with Berendsen pressure coupling. If "
1015 "you are not equilibrating the system, you can probably "
1016 "ignore this warning.",
1017 epcoupl_names[ir->epc]);
1018 warning(wi, warn_buf);
1024 if (ir->epc == epcMTTK)
1026 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1030 /* ELECTROSTATICS */
1031 /* More checks are in triple check (grompp.c) */
1033 if (ir->coulombtype == eelSWITCH)
1035 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1036 "artifacts, advice: use coulombtype = %s",
1037 eel_names[ir->coulombtype],
1038 eel_names[eelRF_ZERO]);
1039 warning(wi, warn_buf);
1042 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1044 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);
1045 warning(wi, warn_buf);
1046 ir->epsilon_rf = ir->epsilon_r;
1047 ir->epsilon_r = 1.0;
1050 if (ir->epsilon_r == 0)
1053 "It is pointless to use long-range electrostatics with infinite relative permittivity."
1054 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1055 CHECK(EEL_FULL(ir->coulombtype));
1058 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1060 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1061 CHECK(ir->epsilon_r < 0);
1064 if (EEL_RF(ir->coulombtype))
1066 /* reaction field (at the cut-off) */
1068 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1070 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1071 eel_names[ir->coulombtype]);
1072 warning(wi, warn_buf);
1073 ir->epsilon_rf = 0.0;
1076 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1077 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1078 (ir->epsilon_r == 0));
1079 if (ir->epsilon_rf == ir->epsilon_r)
1081 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1082 eel_names[ir->coulombtype]);
1083 warning(wi, warn_buf);
1086 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1087 * means the interaction is zero outside rcoulomb, but it helps to
1088 * provide accurate energy conservation.
1090 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1092 if (ir_coulomb_switched(ir))
1095 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1096 eel_names[ir->coulombtype]);
1097 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1101 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1104 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1105 CHECK( ir->coulomb_modifier != eintmodNONE);
1107 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1110 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1111 CHECK( ir->vdw_modifier != eintmodNONE);
1114 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1115 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1118 "The switch/shift interaction settings are just for compatibility; you will get better "
1119 "performance from applying potential modifiers to your interactions!\n");
1120 warning_note(wi, warn_buf);
1123 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1125 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1127 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1128 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.",
1129 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1130 warning(wi, warn_buf);
1134 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1136 if (ir->rvdw_switch == 0)
1138 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.");
1139 warning(wi, warn_buf);
1143 if (EEL_FULL(ir->coulombtype))
1145 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1146 ir->coulombtype == eelPMEUSERSWITCH)
1148 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1149 eel_names[ir->coulombtype]);
1150 CHECK(ir->rcoulomb > ir->rlist);
1154 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1156 // TODO: Move these checks into the ewald module with the options class
1158 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1160 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1162 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
1163 warning_error(wi, warn_buf);
1167 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1169 if (ir->ewald_geometry == eewg3D)
1171 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1172 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1173 warning(wi, warn_buf);
1175 /* This check avoids extra pbc coding for exclusion corrections */
1176 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1177 CHECK(ir->wall_ewald_zfac < 2);
1179 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1180 EEL_FULL(ir->coulombtype))
1182 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1183 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1184 warning(wi, warn_buf);
1186 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1188 sprintf(err_buf, "Cannot have periodic molecules with epsilon_surface > 0");
1189 CHECK(ir->bPeriodicMols);
1190 sprintf(warn_buf, "With epsilon_surface > 0 all molecules should be neutral.");
1191 warning_note(wi, warn_buf);
1193 "With epsilon_surface > 0 you can only use domain decomposition "
1194 "when there are only small molecules with all bonds constrained (mdrun will check for this).");
1195 warning_note(wi, warn_buf);
1198 if (ir_vdw_switched(ir))
1200 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1201 CHECK(ir->rvdw_switch >= ir->rvdw);
1203 if (ir->rvdw_switch < 0.5*ir->rvdw)
1205 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.",
1206 ir->rvdw_switch, ir->rvdw);
1207 warning_note(wi, warn_buf);
1211 if (ir->vdwtype == evdwPME)
1213 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1215 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1216 evdw_names[ir->vdwtype],
1217 eintmod_names[eintmodPOTSHIFT],
1218 eintmod_names[eintmodNONE]);
1219 warning_error(wi, err_buf);
1223 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1225 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.");
1228 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1231 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1234 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1236 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1239 /* IMPLICIT SOLVENT */
1240 if (ir->coulombtype == eelGB_NOTUSED)
1242 sprintf(warn_buf, "Invalid option %s for coulombtype",
1243 eel_names[ir->coulombtype]);
1244 warning_error(wi, warn_buf);
1249 warning_error(wi, "QMMM is currently not supported");
1250 if (!EI_DYNAMICS(ir->eI))
1253 sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
1254 warning_error(wi, buf);
1260 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1264 /* interpret a number of doubles from a string and put them in an array,
1265 after allocating space for them.
1266 str = the input string
1267 n = the (pre-allocated) number of doubles read
1268 r = the output array of doubles. */
1269 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1271 auto values = gmx::splitString(str);
1275 for (int i = 0; i < *n; i++)
1279 (*r)[i] = gmx::fromString<real>(values[i]);
1281 catch (gmx::GromacsException &)
1283 warning_error(wi, "Invalid value " + values[i] + " in string in mdp file. Expected a real number.");
1289 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1292 int i, j, max_n_lambda, nweights, nfep[efptNR];
1293 t_lambda *fep = ir->fepvals;
1294 t_expanded *expand = ir->expandedvals;
1295 real **count_fep_lambdas;
1296 bool bOneLambda = TRUE;
1298 snew(count_fep_lambdas, efptNR);
1300 /* FEP input processing */
1301 /* first, identify the number of lambda values for each type.
1302 All that are nonzero must have the same number */
1304 for (i = 0; i < efptNR; i++)
1306 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1309 /* now, determine the number of components. All must be either zero, or equal. */
1312 for (i = 0; i < efptNR; i++)
1314 if (nfep[i] > max_n_lambda)
1316 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1317 must have the same number if its not zero.*/
1322 for (i = 0; i < efptNR; i++)
1326 ir->fepvals->separate_dvdl[i] = FALSE;
1328 else if (nfep[i] == max_n_lambda)
1330 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1331 respect to the temperature currently */
1333 ir->fepvals->separate_dvdl[i] = TRUE;
1338 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1339 nfep[i], efpt_names[i], max_n_lambda);
1342 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1343 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1345 /* the number of lambdas is the number we've read in, which is either zero
1346 or the same for all */
1347 fep->n_lambda = max_n_lambda;
1349 /* allocate space for the array of lambda values */
1350 snew(fep->all_lambda, efptNR);
1351 /* if init_lambda is defined, we need to set lambda */
1352 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1354 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1356 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1357 for (i = 0; i < efptNR; i++)
1359 snew(fep->all_lambda[i], fep->n_lambda);
1360 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1363 for (j = 0; j < fep->n_lambda; j++)
1365 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1367 sfree(count_fep_lambdas[i]);
1370 sfree(count_fep_lambdas);
1372 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1373 bookkeeping -- for now, init_lambda */
1375 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1377 for (i = 0; i < fep->n_lambda; i++)
1379 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1383 /* check to see if only a single component lambda is defined, and soft core is defined.
1384 In this case, turn on coulomb soft core */
1386 if (max_n_lambda == 0)
1392 for (i = 0; i < efptNR; i++)
1394 if ((nfep[i] != 0) && (i != efptFEP))
1400 if ((bOneLambda) && (fep->sc_alpha > 0))
1402 fep->bScCoul = TRUE;
1405 /* Fill in the others with the efptFEP if they are not explicitly
1406 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1407 they are all zero. */
1409 for (i = 0; i < efptNR; i++)
1411 if ((nfep[i] == 0) && (i != efptFEP))
1413 for (j = 0; j < fep->n_lambda; j++)
1415 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1421 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1422 if (fep->sc_r_power == 48)
1424 if (fep->sc_alpha > 0.1)
1426 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1430 /* now read in the weights */
1431 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1434 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1436 else if (nweights != fep->n_lambda)
1438 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1439 nweights, fep->n_lambda);
1441 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1443 expand->nstexpanded = fep->nstdhdl;
1444 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1449 static void do_simtemp_params(t_inputrec *ir)
1452 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1453 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1457 convertYesNos(warninp_t /*wi*/, gmx::ArrayRef<const std::string> inputs, const char * /*name*/, gmx_bool *outputs)
1460 for (const auto &input : inputs)
1462 outputs[i] = gmx::equalCaseInsensitive(input, "Y", 1);
1467 template <typename T> void
1468 convertInts(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, T *outputs)
1471 for (const auto &input : inputs)
1475 outputs[i] = gmx::fromStdString<T>(input);
1477 catch (gmx::GromacsException &)
1479 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of integers separated by spaces.",
1481 warning_error(wi, message);
1488 convertReals(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, real *outputs)
1491 for (const auto &input : inputs)
1495 outputs[i] = gmx::fromString<real>(input);
1497 catch (gmx::GromacsException &)
1499 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of real numbers separated by spaces.",
1501 warning_error(wi, message);
1508 convertRvecs(warninp_t wi, gmx::ArrayRef<const std::string> inputs, const char *name, rvec *outputs)
1511 for (const auto &input : inputs)
1515 outputs[i][d] = gmx::fromString<real>(input);
1517 catch (gmx::GromacsException &)
1519 auto message = gmx::formatString("Invalid value for mdp option %s. %s should only consist of real numbers separated by spaces.",
1521 warning_error(wi, message);
1532 static void do_wall_params(t_inputrec *ir,
1533 char *wall_atomtype, char *wall_density,
1537 opts->wall_atomtype[0] = nullptr;
1538 opts->wall_atomtype[1] = nullptr;
1540 ir->wall_atomtype[0] = -1;
1541 ir->wall_atomtype[1] = -1;
1542 ir->wall_density[0] = 0;
1543 ir->wall_density[1] = 0;
1547 auto wallAtomTypes = gmx::splitString(wall_atomtype);
1548 if (wallAtomTypes.size() != size_t(ir->nwall))
1550 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %zu",
1551 ir->nwall, wallAtomTypes.size());
1553 for (int i = 0; i < ir->nwall; i++)
1555 opts->wall_atomtype[i] = gmx_strdup(wallAtomTypes[i].c_str());
1558 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1560 auto wallDensity = gmx::splitString(wall_density);
1561 if (wallDensity.size() != size_t(ir->nwall))
1563 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %zu", ir->nwall, wallDensity.size());
1565 convertReals(wi, wallDensity, "wall-density", ir->wall_density);
1566 for (int i = 0; i < ir->nwall; i++)
1568 if (ir->wall_density[i] <= 0)
1570 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, ir->wall_density[i]);
1577 static void add_wall_energrps(SimulationGroups *groups, int nwall, t_symtab *symtab)
1581 AtomGroupIndices *grps = &(groups->groups[SimulationAtomGroupType::EnergyOutput]);
1582 for (int i = 0; i < nwall; i++)
1584 groups->groupNames.emplace_back(
1587 gmx::formatString("wall%d", i).c_str()));
1588 grps->emplace_back(groups->groupNames.size() - 1);
1593 static void read_expandedparams(std::vector<t_inpfile> *inp,
1594 t_expanded *expand, warninp_t wi)
1596 /* read expanded ensemble parameters */
1597 printStringNewline(inp, "expanded ensemble variables");
1598 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1599 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1600 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1601 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1602 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1603 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1604 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1605 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1606 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1607 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1608 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1609 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1610 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1611 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1612 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1613 expand->bSymmetrizedTMatrix = (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1614 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1615 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1616 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1617 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1618 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1619 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1620 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1623 /*! \brief Return whether an end state with the given coupling-lambda
1624 * value describes fully-interacting VDW.
1626 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1627 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1629 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1631 return (couple_lambda_value == ecouplamVDW ||
1632 couple_lambda_value == ecouplamVDWQ);
1638 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1641 explicit MdpErrorHandler(warninp_t wi)
1642 : wi_(wi), mapping_(nullptr)
1646 void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
1648 mapping_ = &mapping;
1651 bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context) override
1653 ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
1654 getOptionName(context).c_str()));
1655 std::string message = gmx::formatExceptionMessageToString(*ex);
1656 warning_error(wi_, message.c_str());
1661 std::string getOptionName(const gmx::KeyValueTreePath &context)
1663 if (mapping_ != nullptr)
1665 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1666 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1669 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1674 const gmx::IKeyValueTreeBackMapping *mapping_;
1679 void get_ir(const char *mdparin, const char *mdparout,
1680 gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
1681 WriteMdpHeader writeMdpHeader, warninp_t wi)
1684 double dumdub[2][6];
1686 char warn_buf[STRLEN];
1687 t_lambda *fep = ir->fepvals;
1688 t_expanded *expand = ir->expandedvals;
1690 const char *no_names[] = { "no", nullptr };
1692 init_inputrec_strings();
1693 gmx::TextInputFile stream(mdparin);
1694 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1696 snew(dumstr[0], STRLEN);
1697 snew(dumstr[1], STRLEN);
1699 /* ignore the following deprecated commands */
1700 replace_inp_entry(inp, "title", nullptr);
1701 replace_inp_entry(inp, "cpp", nullptr);
1702 replace_inp_entry(inp, "domain-decomposition", nullptr);
1703 replace_inp_entry(inp, "andersen-seed", nullptr);
1704 replace_inp_entry(inp, "dihre", nullptr);
1705 replace_inp_entry(inp, "dihre-fc", nullptr);
1706 replace_inp_entry(inp, "dihre-tau", nullptr);
1707 replace_inp_entry(inp, "nstdihreout", nullptr);
1708 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1709 replace_inp_entry(inp, "optimize-fft", nullptr);
1710 replace_inp_entry(inp, "adress_type", nullptr);
1711 replace_inp_entry(inp, "adress_const_wf", nullptr);
1712 replace_inp_entry(inp, "adress_ex_width", nullptr);
1713 replace_inp_entry(inp, "adress_hy_width", nullptr);
1714 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1715 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1716 replace_inp_entry(inp, "adress_site", nullptr);
1717 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1718 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1719 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1720 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1721 replace_inp_entry(inp, "rlistlong", nullptr);
1722 replace_inp_entry(inp, "nstcalclr", nullptr);
1723 replace_inp_entry(inp, "pull-print-com2", nullptr);
1724 replace_inp_entry(inp, "gb-algorithm", nullptr);
1725 replace_inp_entry(inp, "nstgbradii", nullptr);
1726 replace_inp_entry(inp, "rgbradii", nullptr);
1727 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1728 replace_inp_entry(inp, "gb-saltconc", nullptr);
1729 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1730 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1731 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1732 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1733 replace_inp_entry(inp, "sa-algorithm", nullptr);
1734 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1735 replace_inp_entry(inp, "ns-type", nullptr);
1737 /* replace the following commands with the clearer new versions*/
1738 replace_inp_entry(inp, "unconstrained-start", "continuation");
1739 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1740 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1741 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1742 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1743 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1744 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1746 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1747 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1748 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1749 setStringEntry(&inp, "include", opts->include, nullptr);
1750 printStringNoNewline(&inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1751 setStringEntry(&inp, "define", opts->define, nullptr);
1753 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1754 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1755 printStringNoNewline(&inp, "Start time and timestep in ps");
1756 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1757 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1758 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1759 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1760 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1761 printStringNoNewline(&inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1762 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1763 printStringNoNewline(&inp, "mode for center of mass motion removal");
1764 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
1765 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
1766 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
1767 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
1768 setStringEntry(&inp, "comm-grps", is->vcm, nullptr);
1770 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
1771 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
1772 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
1773 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
1776 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
1777 printStringNoNewline(&inp, "Force tolerance and initial step-size");
1778 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
1779 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
1780 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
1781 ir->niter = get_eint(&inp, "niter", 20, wi);
1782 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
1783 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
1784 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
1785 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
1786 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
1788 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
1789 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
1791 /* Output options */
1792 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
1793 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
1794 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
1795 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
1796 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
1797 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
1798 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
1799 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
1800 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
1801 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
1802 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
1803 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
1804 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
1805 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
1806 printStringNoNewline(&inp, "default, all atoms will be written.");
1807 setStringEntry(&inp, "compressed-x-grps", is->x_compressed_groups, nullptr);
1808 printStringNoNewline(&inp, "Selection of energy groups");
1809 setStringEntry(&inp, "energygrps", is->energy, nullptr);
1811 /* Neighbor searching */
1812 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
1813 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1814 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
1815 printStringNoNewline(&inp, "nblist update frequency");
1816 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
1817 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
1818 ir->ePBC = get_eeenum(&inp, "pbc", epbc_names, wi);
1819 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
1820 printStringNoNewline(&inp, "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1821 printStringNoNewline(&inp, "a value of -1 means: use rlist");
1822 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
1823 printStringNoNewline(&inp, "nblist cut-off");
1824 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
1825 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
1827 /* Electrostatics */
1828 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
1829 printStringNoNewline(&inp, "Method for doing electrostatics");
1830 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
1831 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
1832 printStringNoNewline(&inp, "cut-off lengths");
1833 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
1834 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
1835 printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
1836 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
1837 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
1838 printStringNoNewline(&inp, "Method for doing Van der Waals");
1839 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
1840 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
1841 printStringNoNewline(&inp, "cut-off lengths");
1842 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
1843 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
1844 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
1845 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
1846 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
1847 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
1848 printStringNoNewline(&inp, "Separate tables between energy group pairs");
1849 setStringEntry(&inp, "energygrp-table", is->egptable, nullptr);
1850 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
1851 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
1852 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
1853 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
1854 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
1855 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
1856 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
1857 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
1858 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
1859 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
1860 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
1861 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
1862 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
1864 /* Implicit solvation is no longer supported, but we need grompp
1865 to be able to refuse old .mdp files that would have built a tpr
1866 to run it. Thus, only "no" is accepted. */
1867 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
1869 /* Coupling stuff */
1870 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
1871 printStringNoNewline(&inp, "Temperature coupling");
1872 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
1873 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
1874 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
1875 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
1876 printStringNoNewline(&inp, "Groups to couple separately");
1877 setStringEntry(&inp, "tc-grps", is->tcgrps, nullptr);
1878 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
1879 setStringEntry(&inp, "tau-t", is->tau_t, nullptr);
1880 setStringEntry(&inp, "ref-t", is->ref_t, nullptr);
1881 printStringNoNewline(&inp, "pressure coupling");
1882 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
1883 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
1884 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
1885 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
1886 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
1887 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
1888 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
1889 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
1890 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
1893 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
1894 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
1895 printStringNoNewline(&inp, "Groups treated Quantum Mechanically");
1896 setStringEntry(&inp, "QMMM-grps", is->QMMM, nullptr);
1897 printStringNoNewline(&inp, "QM method");
1898 setStringEntry(&inp, "QMmethod", is->QMmethod, nullptr);
1899 printStringNoNewline(&inp, "QMMM scheme");
1900 ir->QMMMscheme = get_eeenum(&inp, "QMMMscheme", eQMMMscheme_names, wi);
1901 printStringNoNewline(&inp, "QM basisset");
1902 setStringEntry(&inp, "QMbasis", is->QMbasis, nullptr);
1903 printStringNoNewline(&inp, "QM charge");
1904 setStringEntry(&inp, "QMcharge", is->QMcharge, nullptr);
1905 printStringNoNewline(&inp, "QM multiplicity");
1906 setStringEntry(&inp, "QMmult", is->QMmult, nullptr);
1907 printStringNoNewline(&inp, "Surface Hopping");
1908 setStringEntry(&inp, "SH", is->bSH, nullptr);
1909 printStringNoNewline(&inp, "CAS space options");
1910 setStringEntry(&inp, "CASorbitals", is->CASorbitals, nullptr);
1911 setStringEntry(&inp, "CASelectrons", is->CASelectrons, nullptr);
1912 setStringEntry(&inp, "SAon", is->SAon, nullptr);
1913 setStringEntry(&inp, "SAoff", is->SAoff, nullptr);
1914 setStringEntry(&inp, "SAsteps", is->SAsteps, nullptr);
1915 printStringNoNewline(&inp, "Scale factor for MM charges");
1916 ir->scalefactor = get_ereal(&inp, "MMChargeScaleFactor", 1.0, wi);
1918 /* Simulated annealing */
1919 printStringNewline(&inp, "SIMULATED ANNEALING");
1920 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
1921 setStringEntry(&inp, "annealing", is->anneal, nullptr);
1922 printStringNoNewline(&inp, "Number of time points to use for specifying annealing in each group");
1923 setStringEntry(&inp, "annealing-npoints", is->anneal_npoints, nullptr);
1924 printStringNoNewline(&inp, "List of times at the annealing points for each group");
1925 setStringEntry(&inp, "annealing-time", is->anneal_time, nullptr);
1926 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
1927 setStringEntry(&inp, "annealing-temp", is->anneal_temp, nullptr);
1930 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
1931 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
1932 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
1933 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
1936 printStringNewline(&inp, "OPTIONS FOR BONDS");
1937 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
1938 printStringNoNewline(&inp, "Type of constraint algorithm");
1939 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
1940 printStringNoNewline(&inp, "Do not constrain the start configuration");
1941 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
1942 printStringNoNewline(&inp, "Use successive overrelaxation to reduce the number of shake iterations");
1943 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
1944 printStringNoNewline(&inp, "Relative tolerance of shake");
1945 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
1946 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
1947 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
1948 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
1949 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
1950 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
1951 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
1952 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
1953 printStringNoNewline(&inp, "rotates over more degrees than");
1954 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
1955 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
1956 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
1958 /* Energy group exclusions */
1959 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
1960 printStringNoNewline(&inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
1961 setStringEntry(&inp, "energygrp-excl", is->egpexcl, nullptr);
1964 printStringNewline(&inp, "WALLS");
1965 printStringNoNewline(&inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1966 ir->nwall = get_eint(&inp, "nwall", 0, wi);
1967 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
1968 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
1969 setStringEntry(&inp, "wall-atomtype", is->wall_atomtype, nullptr);
1970 setStringEntry(&inp, "wall-density", is->wall_density, nullptr);
1971 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
1974 printStringNewline(&inp, "COM PULLING");
1975 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
1979 is->pull_grp = read_pullparams(&inp, ir->pull, wi);
1983 NOTE: needs COM pulling input */
1984 printStringNewline(&inp, "AWH biasing");
1985 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
1990 ir->awhParams = gmx::readAndCheckAwhParams(&inp, ir, wi);
1994 gmx_fatal(FARGS, "AWH biasing is only compatible with COM pulling turned on");
1998 /* Enforced rotation */
1999 printStringNewline(&inp, "ENFORCED ROTATION");
2000 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2001 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2005 is->rot_grp = read_rotparams(&inp, ir->rot, wi);
2008 /* Interactive MD */
2010 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2011 setStringEntry(&inp, "IMD-group", is->imd_grp, nullptr);
2012 if (is->imd_grp[0] != '\0')
2019 printStringNewline(&inp, "NMR refinement stuff");
2020 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2021 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2022 printStringNoNewline(&inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2023 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2024 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2025 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2026 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2027 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2028 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2029 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2030 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2031 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2032 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2033 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2034 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2035 setStringEntry(&inp, "orire-fitgrp", is->orirefitgrp, nullptr);
2036 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2037 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2039 /* free energy variables */
2040 printStringNewline(&inp, "Free energy variables");
2041 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2042 setStringEntry(&inp, "couple-moltype", is->couple_moltype, nullptr);
2043 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2044 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2045 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2047 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2049 it was not entered */
2050 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2051 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2052 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2053 setStringEntry(&inp, "fep-lambdas", is->fep_lambda[efptFEP], nullptr);
2054 setStringEntry(&inp, "mass-lambdas", is->fep_lambda[efptMASS], nullptr);
2055 setStringEntry(&inp, "coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
2056 setStringEntry(&inp, "vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
2057 setStringEntry(&inp, "bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
2058 setStringEntry(&inp, "restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
2059 setStringEntry(&inp, "temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
2060 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2061 setStringEntry(&inp, "init-lambda-weights", is->lambda_weights, nullptr);
2062 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2063 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2064 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2065 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2066 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2067 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2068 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2069 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2070 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2071 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2072 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2073 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2075 /* Non-equilibrium MD stuff */
2076 printStringNewline(&inp, "Non-equilibrium MD stuff");
2077 setStringEntry(&inp, "acc-grps", is->accgrps, nullptr);
2078 setStringEntry(&inp, "accelerate", is->acc, nullptr);
2079 setStringEntry(&inp, "freezegrps", is->freeze, nullptr);
2080 setStringEntry(&inp, "freezedim", is->frdim, nullptr);
2081 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2082 setStringEntry(&inp, "deform", is->deform, nullptr);
2084 /* simulated tempering variables */
2085 printStringNewline(&inp, "simulated tempering variables");
2086 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2087 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2088 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2089 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2091 /* expanded ensemble variables */
2092 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2094 read_expandedparams(&inp, expand, wi);
2097 /* Electric fields */
2099 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2100 gmx::KeyValueTreeTransformer transform;
2101 transform.rules()->addRule()
2102 .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2103 mdModules->initMdpTransform(transform.rules());
2104 for (const auto &path : transform.mappedPaths())
2106 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2107 mark_einp_set(inp, path[0].c_str());
2109 MdpErrorHandler errorHandler(wi);
2111 = transform.transform(convertedValues, &errorHandler);
2112 ir->params = new gmx::KeyValueTreeObject(result.object());
2113 mdModules->adjustInputrecBasedOnModules(ir);
2114 errorHandler.setBackMapping(result.backMapping());
2115 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2118 /* Ion/water position swapping ("computational electrophysiology") */
2119 printStringNewline(&inp, "Ion/water position swapping for computational electrophysiology setups");
2120 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2121 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2122 if (ir->eSwapCoords != eswapNO)
2129 printStringNoNewline(&inp, "Swap attempt frequency");
2130 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2131 printStringNoNewline(&inp, "Number of ion types to be controlled");
2132 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2135 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2137 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2138 snew(ir->swap->grp, ir->swap->ngrp);
2139 for (i = 0; i < ir->swap->ngrp; i++)
2141 snew(ir->swap->grp[i].molname, STRLEN);
2143 printStringNoNewline(&inp, "Two index groups that contain the compartment-partitioning atoms");
2144 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2145 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2146 printStringNoNewline(&inp, "Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2147 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2148 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2150 printStringNoNewline(&inp, "Name of solvent molecules");
2151 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2153 printStringNoNewline(&inp, "Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2154 printStringNoNewline(&inp, "Note that the split cylinder settings do not have an influence on the swapping protocol,");
2155 printStringNoNewline(&inp, "however, if correctly defined, the permeation events are recorded per channel");
2156 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2157 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2158 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2159 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2160 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2161 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2163 printStringNoNewline(&inp, "Average the number of ions per compartment over these many swap attempt steps");
2164 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2166 printStringNoNewline(&inp, "Names of the ion types that can be exchanged with solvent molecules,");
2167 printStringNoNewline(&inp, "and the requested number of ions of this type in compartments A and B");
2168 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2169 for (i = 0; i < nIonTypes; i++)
2171 int ig = eSwapFixedGrpNR + i;
2173 sprintf(buf, "iontype%d-name", i);
2174 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2175 sprintf(buf, "iontype%d-in-A", i);
2176 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2177 sprintf(buf, "iontype%d-in-B", i);
2178 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2181 printStringNoNewline(&inp, "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2182 printStringNoNewline(&inp, "at maximum distance (= bulk concentration) to the split group layers. However,");
2183 printStringNoNewline(&inp, "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2184 printStringNoNewline(&inp, "towards one of the compartment-partitioning layers (at +/- 1.0).");
2185 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2186 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2187 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2188 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2190 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2193 printStringNoNewline(&inp, "Start to swap ions if threshold difference to requested count is reached");
2194 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2197 /* AdResS is no longer supported, but we need grompp to be able to
2198 refuse to process old .mdp files that used it. */
2199 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2201 /* User defined thingies */
2202 printStringNewline(&inp, "User defined thingies");
2203 setStringEntry(&inp, "user1-grps", is->user1, nullptr);
2204 setStringEntry(&inp, "user2-grps", is->user2, nullptr);
2205 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2206 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2207 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2208 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2209 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2210 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2211 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2212 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2216 gmx::TextOutputFile stream(mdparout);
2217 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2219 // Transform module data into a flat key-value tree for output.
2220 gmx::KeyValueTreeBuilder builder;
2221 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2222 mdModules->buildMdpOutput(&builderObject);
2224 gmx::TextWriter writer(&stream);
2225 writeKeyValueTreeAsMdp(&writer, builder.build());
2230 /* Process options if necessary */
2231 for (m = 0; m < 2; m++)
2233 for (i = 0; i < 2*DIM; i++)
2242 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2244 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
2246 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2248 case epctSEMIISOTROPIC:
2249 case epctSURFACETENSION:
2250 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2252 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
2254 dumdub[m][YY] = dumdub[m][XX];
2256 case epctANISOTROPIC:
2257 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2258 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2259 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2261 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
2265 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2266 epcoupltype_names[ir->epct]);
2270 clear_mat(ir->ref_p);
2271 clear_mat(ir->compress);
2272 for (i = 0; i < DIM; i++)
2274 ir->ref_p[i][i] = dumdub[1][i];
2275 ir->compress[i][i] = dumdub[0][i];
2277 if (ir->epct == epctANISOTROPIC)
2279 ir->ref_p[XX][YY] = dumdub[1][3];
2280 ir->ref_p[XX][ZZ] = dumdub[1][4];
2281 ir->ref_p[YY][ZZ] = dumdub[1][5];
2282 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2284 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2286 ir->compress[XX][YY] = dumdub[0][3];
2287 ir->compress[XX][ZZ] = dumdub[0][4];
2288 ir->compress[YY][ZZ] = dumdub[0][5];
2289 for (i = 0; i < DIM; i++)
2291 for (m = 0; m < i; m++)
2293 ir->ref_p[i][m] = ir->ref_p[m][i];
2294 ir->compress[i][m] = ir->compress[m][i];
2299 if (ir->comm_mode == ecmNO)
2304 opts->couple_moltype = nullptr;
2305 if (strlen(is->couple_moltype) > 0)
2307 if (ir->efep != efepNO)
2309 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2310 if (opts->couple_lam0 == opts->couple_lam1)
2312 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2314 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2315 opts->couple_lam1 == ecouplamNONE))
2317 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2322 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2325 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2326 if (ir->efep != efepNO)
2328 if (fep->delta_lambda > 0)
2330 ir->efep = efepSLOWGROWTH;
2334 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2336 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2337 warning_note(wi, "Old option for dhdl-print-energy given: "
2338 "changing \"yes\" to \"total\"\n");
2341 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2343 /* always print out the energy to dhdl if we are doing
2344 expanded ensemble, since we need the total energy for
2345 analysis if the temperature is changing. In some
2346 conditions one may only want the potential energy, so
2347 we will allow that if the appropriate mdp setting has
2348 been enabled. Otherwise, total it is:
2350 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2353 if ((ir->efep != efepNO) || ir->bSimTemp)
2355 ir->bExpanded = FALSE;
2356 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2358 ir->bExpanded = TRUE;
2360 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2361 if (ir->bSimTemp) /* done after fep params */
2363 do_simtemp_params(ir);
2366 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2367 * setup and not on the old way of specifying the free-energy setup,
2368 * we should check for using soft-core when not needed, since that
2369 * can complicate the sampling significantly.
2370 * Note that we only check for the automated coupling setup.
2371 * If the (advanced) user does FEP through manual topology changes,
2372 * this check will not be triggered.
2374 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2375 ir->fepvals->sc_alpha != 0 &&
2376 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2377 couple_lambda_has_vdw_on(opts->couple_lam1)))
2379 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.");
2384 ir->fepvals->n_lambda = 0;
2387 /* WALL PARAMETERS */
2389 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts, wi);
2391 /* ORIENTATION RESTRAINT PARAMETERS */
2393 if (opts->bOrire && gmx::splitString(is->orirefitgrp).size() != 1)
2395 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2398 /* DEFORMATION PARAMETERS */
2400 clear_mat(ir->deform);
2401 for (i = 0; i < 6; i++)
2406 double gmx_unused canary;
2407 int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
2408 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2409 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
2411 if (strlen(is->deform) > 0 && ndeform != 6)
2413 warning_error(wi, gmx::formatString("Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform).c_str());
2415 for (i = 0; i < 3; i++)
2417 ir->deform[i][i] = dumdub[0][i];
2419 ir->deform[YY][XX] = dumdub[0][3];
2420 ir->deform[ZZ][XX] = dumdub[0][4];
2421 ir->deform[ZZ][YY] = dumdub[0][5];
2422 if (ir->epc != epcNO)
2424 for (i = 0; i < 3; i++)
2426 for (j = 0; j <= i; j++)
2428 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2430 warning_error(wi, "A box element has deform set and compressibility > 0");
2434 for (i = 0; i < 3; i++)
2436 for (j = 0; j < i; j++)
2438 if (ir->deform[i][j] != 0)
2440 for (m = j; m < DIM; m++)
2442 if (ir->compress[m][j] != 0)
2444 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.");
2445 warning(wi, warn_buf);
2453 /* Ion/water position swapping checks */
2454 if (ir->eSwapCoords != eswapNO)
2456 if (ir->swap->nstswap < 1)
2458 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2460 if (ir->swap->nAverage < 1)
2462 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2464 if (ir->swap->threshold < 1.0)
2466 warning_error(wi, "Ion count threshold must be at least 1.\n");
2474 static int search_QMstring(const char *s, int ng, const char *gn[])
2476 /* same as normal search_string, but this one searches QM strings */
2479 for (i = 0; (i < ng); i++)
2481 if (gmx_strcasecmp(s, gn[i]) == 0)
2487 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2488 } /* search_QMstring */
2490 /* We would like gn to be const as well, but C doesn't allow this */
2491 /* TODO this is utility functionality (search for the index of a
2492 string in a collection), so should be refactored and located more
2494 int search_string(const char *s, int ng, char *gn[])
2498 for (i = 0; (i < ng); i++)
2500 if (gmx_strcasecmp(s, gn[i]) == 0)
2507 "Group %s referenced in the .mdp file was not found in the index file.\n"
2508 "Group names must match either [moleculetype] names or custom index group\n"
2509 "names, in which case you must supply an index file to the '-n' option\n"
2514 static bool do_numbering(int natoms, SimulationGroups *groups,
2515 gmx::ArrayRef<std::string> groupsFromMdpFile,
2516 t_blocka *block, char *gnames[],
2517 SimulationAtomGroupType gtype, int restnm,
2518 int grptp, bool bVerbose,
2521 unsigned short *cbuf;
2522 AtomGroupIndices *grps = &(groups->groups[gtype]);
2523 int j, gid, aj, ognr, ntot = 0;
2526 char warn_buf[STRLEN];
2528 title = shortName(gtype);
2531 /* Mark all id's as not set */
2532 for (int i = 0; (i < natoms); i++)
2537 for (int i = 0; i != groupsFromMdpFile.ssize(); ++i)
2539 /* Lookup the group name in the block structure */
2540 gid = search_string(groupsFromMdpFile[i].c_str(), block->nr, gnames);
2541 if ((grptp != egrptpONE) || (i == 0))
2543 grps->emplace_back(gid);
2546 /* Now go over the atoms in the group */
2547 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2552 /* Range checking */
2553 if ((aj < 0) || (aj >= natoms))
2555 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj + 1);
2557 /* Lookup up the old group number */
2561 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2562 aj+1, title, ognr+1, i+1);
2566 /* Store the group number in buffer */
2567 if (grptp == egrptpONE)
2580 /* Now check whether we have done all atoms */
2584 if (grptp == egrptpALL)
2586 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2587 natoms-ntot, title);
2589 else if (grptp == egrptpPART)
2591 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2592 natoms-ntot, title);
2593 warning_note(wi, warn_buf);
2595 /* Assign all atoms currently unassigned to a rest group */
2596 for (j = 0; (j < natoms); j++)
2598 if (cbuf[j] == NOGID)
2600 cbuf[j] = grps->size();
2604 if (grptp != egrptpPART)
2609 "Making dummy/rest group for %s containing %d elements\n",
2610 title, natoms-ntot);
2612 /* Add group name "rest" */
2613 grps->emplace_back(restnm);
2615 /* Assign the rest name to all atoms not currently assigned to a group */
2616 for (j = 0; (j < natoms); j++)
2618 if (cbuf[j] == NOGID)
2620 // group size was not updated before this here, so need to use -1.
2621 cbuf[j] = grps->size() - 1;
2627 if (grps->size() == 1 && (ntot == 0 || ntot == natoms))
2629 /* All atoms are part of one (or no) group, no index required */
2630 groups->groupNumbers[gtype].clear();
2634 for (int j = 0; (j < natoms); j++)
2636 groups->groupNumbers[gtype].emplace_back(cbuf[j]);
2642 return (bRest && grptp == egrptpPART);
2645 static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2648 pull_params_t *pull;
2649 int natoms, imin, jmin;
2650 int *nrdf2, *na_vcm, na_tot;
2651 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2656 * First calc 3xnr-atoms for each group
2657 * then subtract half a degree of freedom for each constraint
2659 * Only atoms and nuclei contribute to the degrees of freedom...
2664 const SimulationGroups &groups = mtop->groups;
2665 natoms = mtop->natoms;
2667 /* Allocate one more for a possible rest group */
2668 /* We need to sum degrees of freedom into doubles,
2669 * since floats give too low nrdf's above 3 million atoms.
2671 snew(nrdf_tc, groups.groups[SimulationAtomGroupType::TemperatureCoupling].size()+1);
2672 snew(nrdf_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2673 snew(dof_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2674 snew(na_vcm, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2675 snew(nrdf_vcm_sub, groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval].size()+1);
2677 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2681 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; i++)
2684 clear_ivec(dof_vcm[i]);
2686 nrdf_vcm_sub[i] = 0;
2688 snew(nrdf2, natoms);
2689 for (const AtomProxy atomP : AtomRange(*mtop))
2691 const t_atom &local = atomP.atom();
2692 int i = atomP.globalAtomNumber();
2694 if (local.ptype == eptAtom || local.ptype == eptNucleus)
2696 int g = getGroupType(groups, SimulationAtomGroupType::Freeze, i);
2697 for (int d = 0; d < DIM; d++)
2699 if (opts->nFreeze[g][d] == 0)
2701 /* Add one DOF for particle i (counted as 2*1) */
2703 /* VCM group i has dim d as a DOF */
2704 dof_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)][d] = 1;
2707 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, i)] += 0.5*nrdf2[i];
2708 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, i)] += 0.5*nrdf2[i];
2713 for (const gmx_molblock_t &molb : mtop->molblock)
2715 const gmx_moltype_t &molt = mtop->moltype[molb.type];
2716 const t_atom *atom = molt.atoms.atom;
2717 for (int mol = 0; mol < molb.nmol; mol++)
2719 for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2721 gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
2722 for (int i = 0; i < molt.ilist[ftype].size(); )
2724 /* Subtract degrees of freedom for the constraints,
2725 * if the particles still have degrees of freedom left.
2726 * If one of the particles is a vsite or a shell, then all
2727 * constraint motion will go there, but since they do not
2728 * contribute to the constraints the degrees of freedom do not
2731 int ai = as + ia[i + 1];
2732 int aj = as + ia[i + 2];
2733 if (((atom[ia[i + 1]].ptype == eptNucleus) ||
2734 (atom[ia[i + 1]].ptype == eptAtom)) &&
2735 ((atom[ia[i + 2]].ptype == eptNucleus) ||
2736 (atom[ia[i + 2]].ptype == eptAtom)))
2754 imin = std::min(imin, nrdf2[ai]);
2755 jmin = std::min(jmin, nrdf2[aj]);
2758 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2759 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, aj)] -= 0.5*jmin;
2760 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2761 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, aj)] -= 0.5*jmin;
2763 i += interaction_function[ftype].nratoms+1;
2766 gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
2767 for (int i = 0; i < molt.ilist[F_SETTLE].size(); )
2769 /* Subtract 1 dof from every atom in the SETTLE */
2770 for (int j = 0; j < 3; j++)
2772 int ai = as + ia[i + 1 + j];
2773 imin = std::min(2, nrdf2[ai]);
2775 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2776 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2780 as += molt.atoms.nr;
2786 /* Correct nrdf for the COM constraints.
2787 * We correct using the TC and VCM group of the first atom
2788 * in the reference and pull group. If atoms in one pull group
2789 * belong to different TC or VCM groups it is anyhow difficult
2790 * to determine the optimal nrdf assignment.
2794 for (int i = 0; i < pull->ncoord; i++)
2796 if (pull->coord[i].eType != epullCONSTRAINT)
2803 for (int j = 0; j < 2; j++)
2805 const t_pull_group *pgrp;
2807 pgrp = &pull->group[pull->coord[i].group[j]];
2811 /* Subtract 1/2 dof from each group */
2812 int ai = pgrp->ind[0];
2813 nrdf_tc [getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] -= 0.5*imin;
2814 nrdf_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)] -= 0.5*imin;
2815 if (nrdf_tc[getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai)] < 0)
2817 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)]]);
2822 /* We need to subtract the whole DOF from group j=1 */
2829 if (ir->nstcomm != 0)
2833 /* We remove COM motion up to dim ndof_com() */
2834 ndim_rm_vcm = ndof_com(ir);
2836 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
2837 * the number of degrees of freedom in each vcm group when COM
2838 * translation is removed and 6 when rotation is removed as well.
2840 for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2842 switch (ir->comm_mode)
2845 case ecmLINEAR_ACCELERATION_CORRECTION:
2846 nrdf_vcm_sub[j] = 0;
2847 for (int d = 0; d < ndim_rm_vcm; d++)
2856 nrdf_vcm_sub[j] = 6;
2859 gmx_incons("Checking comm_mode");
2863 for (gmx::index i = 0; i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling]); i++)
2865 /* Count the number of atoms of TC group i for every VCM group */
2866 for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2871 for (int ai = 0; ai < natoms; ai++)
2873 if (getGroupType(groups, SimulationAtomGroupType::TemperatureCoupling, ai) == i)
2875 na_vcm[getGroupType(groups, SimulationAtomGroupType::MassCenterVelocityRemoval, ai)]++;
2879 /* Correct for VCM removal according to the fraction of each VCM
2880 * group present in this TC group.
2882 nrdf_uc = nrdf_tc[i];
2884 for (gmx::index j = 0; j < gmx::ssize(groups.groups[SimulationAtomGroupType::MassCenterVelocityRemoval])+1; j++)
2886 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
2888 nrdf_tc[i] += nrdf_uc*(static_cast<double>(na_vcm[j])/static_cast<double>(na_tot))*
2889 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
2894 for (int i = 0; (i < gmx::ssize(groups.groups[SimulationAtomGroupType::TemperatureCoupling])); i++)
2896 opts->nrdf[i] = nrdf_tc[i];
2897 if (opts->nrdf[i] < 0)
2902 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2903 gnames[groups.groups[SimulationAtomGroupType::TemperatureCoupling][i]], opts->nrdf[i]);
2911 sfree(nrdf_vcm_sub);
2914 static bool do_egp_flag(t_inputrec *ir, SimulationGroups *groups,
2915 const char *option, const char *val, int flag)
2917 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2918 * But since this is much larger than STRLEN, such a line can not be parsed.
2919 * The real maximum is the number of names that fit in a string: STRLEN/2.
2921 #define EGP_MAX (STRLEN/2)
2925 auto names = gmx::splitString(val);
2926 if (names.size() % 2 != 0)
2928 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2930 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
2932 for (size_t i = 0; i < names.size() / 2; i++)
2934 // TODO this needs to be replaced by a solution using std::find_if
2937 gmx_strcasecmp(names[2*i].c_str(), *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][j]])))
2943 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2944 names[2*i].c_str(), option);
2948 gmx_strcasecmp(names[2*i+1].c_str(), *(groups->groupNames[groups->groups[SimulationAtomGroupType::EnergyOutput][k]])))
2954 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2955 names[2*i+1].c_str(), option);
2957 if ((j < nr) && (k < nr))
2959 ir->opts.egp_flags[nr*j+k] |= flag;
2960 ir->opts.egp_flags[nr*k+j] |= flag;
2969 static void make_swap_groups(
2974 int ig = -1, i = 0, gind;
2978 /* Just a quick check here, more thorough checks are in mdrun */
2979 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
2981 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
2984 /* Get the index atoms of the split0, split1, solvent, and swap groups */
2985 for (ig = 0; ig < swap->ngrp; ig++)
2987 swapg = &swap->grp[ig];
2988 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
2989 swapg->nat = grps->index[gind+1] - grps->index[gind];
2993 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
2994 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
2995 swap->grp[ig].molname, swapg->nat);
2996 snew(swapg->ind, swapg->nat);
2997 for (i = 0; i < swapg->nat; i++)
2999 swapg->ind[i] = grps->a[grps->index[gind]+i];
3004 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3010 static void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3015 ig = search_string(IMDgname, grps->nr, gnames);
3016 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3018 if (IMDgroup->nat > 0)
3020 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3021 IMDgname, IMDgroup->nat);
3022 snew(IMDgroup->ind, IMDgroup->nat);
3023 for (i = 0; i < IMDgroup->nat; i++)
3025 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3030 void do_index(const char* mdparin, const char *ndx,
3033 const gmx::MdModulesNotifier ¬ifier,
3037 t_blocka *defaultIndexGroups;
3041 char warnbuf[STRLEN], **gnames;
3045 int i, j, k, restnm;
3046 bool bExcl, bTable, bAnneal, bRest;
3047 char warn_buf[STRLEN];
3051 fprintf(stderr, "processing index file...\n");
3055 snew(defaultIndexGroups, 1);
3056 snew(defaultIndexGroups->index, 1);
3058 atoms_all = gmx_mtop_global_atoms(mtop);
3059 analyse(&atoms_all, defaultIndexGroups, &gnames, FALSE, TRUE);
3060 done_atom(&atoms_all);
3064 defaultIndexGroups = init_index(ndx, &gnames);
3067 SimulationGroups *groups = &mtop->groups;
3068 natoms = mtop->natoms;
3069 symtab = &mtop->symtab;
3071 for (int i = 0; (i < defaultIndexGroups->nr); i++)
3073 groups->groupNames.emplace_back(put_symtab(symtab, gnames[i]));
3075 groups->groupNames.emplace_back(put_symtab(symtab, "rest"));
3076 restnm = groups->groupNames.size() - 1;
3077 GMX_RELEASE_ASSERT(restnm == defaultIndexGroups->nr, "Size of allocations must match");
3078 srenew(gnames, defaultIndexGroups->nr+1);
3079 gnames[restnm] = *(groups->groupNames.back());
3081 set_warning_line(wi, mdparin, -1);
3083 auto temperatureCouplingTauValues = gmx::splitString(is->tau_t);
3084 auto temperatureCouplingReferenceValues = gmx::splitString(is->ref_t);
3085 auto temperatureCouplingGroupNames = gmx::splitString(is->tcgrps);
3086 if (temperatureCouplingTauValues.size() != temperatureCouplingGroupNames.size() ||
3087 temperatureCouplingReferenceValues.size() != temperatureCouplingGroupNames.size())
3089 gmx_fatal(FARGS, "Invalid T coupling input: %zu groups, %zu ref-t values and "
3091 temperatureCouplingGroupNames.size(),
3092 temperatureCouplingReferenceValues.size(),
3093 temperatureCouplingTauValues.size());
3096 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3097 do_numbering(natoms, groups, temperatureCouplingGroupNames, defaultIndexGroups, gnames,
3098 SimulationAtomGroupType::TemperatureCoupling,
3099 restnm, useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3100 nr = groups->groups[SimulationAtomGroupType::TemperatureCoupling].size();
3102 snew(ir->opts.nrdf, nr);
3103 snew(ir->opts.tau_t, nr);
3104 snew(ir->opts.ref_t, nr);
3105 if (ir->eI == eiBD && ir->bd_fric == 0)
3107 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3110 if (useReferenceTemperature)
3112 if (size_t(nr) != temperatureCouplingReferenceValues.size())
3114 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3118 convertReals(wi, temperatureCouplingTauValues, "tau-t", ir->opts.tau_t);
3119 for (i = 0; (i < nr); i++)
3121 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3123 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3124 warning_error(wi, warn_buf);
3127 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3129 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.");
3132 if (ir->opts.tau_t[i] >= 0)
3134 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3137 if (ir->etc != etcNO && ir->nsttcouple == -1)
3139 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3144 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3146 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");
3148 if (ir->epc == epcMTTK)
3150 if (ir->etc != etcNOSEHOOVER)
3152 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3156 if (ir->nstpcouple != ir->nsttcouple)
3158 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3159 ir->nstpcouple = ir->nsttcouple = mincouple;
3160 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);
3161 warning_note(wi, warn_buf);
3166 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3167 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3169 if (ETC_ANDERSEN(ir->etc))
3171 if (ir->nsttcouple != 1)
3174 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");
3175 warning_note(wi, warn_buf);
3178 nstcmin = tcouple_min_integration_steps(ir->etc);
3181 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3183 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3184 ETCOUPLTYPE(ir->etc),
3186 ir->nsttcouple*ir->delta_t);
3187 warning(wi, warn_buf);
3190 convertReals(wi, temperatureCouplingReferenceValues, "ref-t", ir->opts.ref_t);
3191 for (i = 0; (i < nr); i++)
3193 if (ir->opts.ref_t[i] < 0)
3195 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3198 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3199 if we are in this conditional) if mc_temp is negative */
3200 if (ir->expandedvals->mc_temp < 0)
3202 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3206 /* Simulated annealing for each group. There are nr groups */
3207 auto simulatedAnnealingGroupNames = gmx::splitString(is->anneal);
3208 if (simulatedAnnealingGroupNames.size() == 1 &&
3209 gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[0], "N", 1))
3211 simulatedAnnealingGroupNames.resize(0);
3213 if (!simulatedAnnealingGroupNames.empty() &&
3214 gmx::ssize(simulatedAnnealingGroupNames) != nr)
3216 gmx_fatal(FARGS, "Wrong number of annealing values: %zu (for %d groups)\n",
3217 simulatedAnnealingGroupNames.size(), nr);
3221 snew(ir->opts.annealing, nr);
3222 snew(ir->opts.anneal_npoints, nr);
3223 snew(ir->opts.anneal_time, nr);
3224 snew(ir->opts.anneal_temp, nr);
3225 for (i = 0; i < nr; i++)
3227 ir->opts.annealing[i] = eannNO;
3228 ir->opts.anneal_npoints[i] = 0;
3229 ir->opts.anneal_time[i] = nullptr;
3230 ir->opts.anneal_temp[i] = nullptr;
3232 if (!simulatedAnnealingGroupNames.empty())
3235 for (i = 0; i < nr; i++)
3237 if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "N", 1))
3239 ir->opts.annealing[i] = eannNO;
3241 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "S", 1))
3243 ir->opts.annealing[i] = eannSINGLE;
3246 else if (gmx::equalCaseInsensitive(simulatedAnnealingGroupNames[i], "P", 1))
3248 ir->opts.annealing[i] = eannPERIODIC;
3254 /* Read the other fields too */
3255 auto simulatedAnnealingPoints = gmx::splitString(is->anneal_npoints);
3256 if (simulatedAnnealingPoints.size() != simulatedAnnealingGroupNames.size())
3258 gmx_fatal(FARGS, "Found %zu annealing-npoints values for %zu groups\n",
3259 simulatedAnnealingPoints.size(), simulatedAnnealingGroupNames.size());
3261 convertInts(wi, simulatedAnnealingPoints, "annealing points", ir->opts.anneal_npoints);
3262 size_t numSimulatedAnnealingFields = 0;
3263 for (i = 0; i < nr; i++)
3265 if (ir->opts.anneal_npoints[i] == 1)
3267 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3269 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3270 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3271 numSimulatedAnnealingFields += ir->opts.anneal_npoints[i];
3274 auto simulatedAnnealingTimes = gmx::splitString(is->anneal_time);
3276 if (simulatedAnnealingTimes.size() != numSimulatedAnnealingFields)
3278 gmx_fatal(FARGS, "Found %zu annealing-time values, wanted %zu\n",
3279 simulatedAnnealingTimes.size(), numSimulatedAnnealingFields);
3281 auto simulatedAnnealingTemperatures = gmx::splitString(is->anneal_temp);
3282 if (simulatedAnnealingTemperatures.size() != numSimulatedAnnealingFields)
3284 gmx_fatal(FARGS, "Found %zu annealing-temp values, wanted %zu\n",
3285 simulatedAnnealingTemperatures.size(), numSimulatedAnnealingFields);
3288 std::vector<real> allSimulatedAnnealingTimes(numSimulatedAnnealingFields);
3289 std::vector<real> allSimulatedAnnealingTemperatures(numSimulatedAnnealingFields);
3290 convertReals(wi, simulatedAnnealingTimes, "anneal-time", allSimulatedAnnealingTimes.data());
3291 convertReals(wi, simulatedAnnealingTemperatures, "anneal-temp", allSimulatedAnnealingTemperatures.data());
3292 for (i = 0, k = 0; i < nr; i++)
3294 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3296 ir->opts.anneal_time[i][j] = allSimulatedAnnealingTimes[k];
3297 ir->opts.anneal_temp[i][j] = allSimulatedAnnealingTemperatures[k];
3300 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3302 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3308 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3310 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3311 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3314 if (ir->opts.anneal_temp[i][j] < 0)
3316 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3321 /* Print out some summary information, to make sure we got it right */
3322 for (i = 0; i < nr; i++)
3324 if (ir->opts.annealing[i] != eannNO)
3326 j = groups->groups[SimulationAtomGroupType::TemperatureCoupling][i];
3327 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3328 *(groups->groupNames[j]), eann_names[ir->opts.annealing[i]],
3329 ir->opts.anneal_npoints[i]);
3330 fprintf(stderr, "Time (ps) Temperature (K)\n");
3331 /* All terms except the last one */
3332 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3334 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3337 /* Finally the last one */
3338 j = ir->opts.anneal_npoints[i]-1;
3339 if (ir->opts.annealing[i] == eannSINGLE)
3341 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3345 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3346 if (std::fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3348 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3359 make_pull_groups(ir->pull, is->pull_grp, defaultIndexGroups, gnames);
3361 make_pull_coords(ir->pull);
3366 make_rotation_groups(ir->rot, is->rot_grp, defaultIndexGroups, gnames);
3369 if (ir->eSwapCoords != eswapNO)
3371 make_swap_groups(ir->swap, defaultIndexGroups, gnames);
3374 /* Make indices for IMD session */
3377 make_IMD_group(ir->imd, is->imd_grp, defaultIndexGroups, gnames);
3380 gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
3381 *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
3382 notifier.notifier_.notify(defaultIndexGroupsAndNames);
3384 auto accelerations = gmx::splitString(is->acc);
3385 auto accelerationGroupNames = gmx::splitString(is->accgrps);
3386 if (accelerationGroupNames.size() * DIM != accelerations.size())
3388 gmx_fatal(FARGS, "Invalid Acceleration input: %zu groups and %zu acc. values",
3389 accelerationGroupNames.size(), accelerations.size());
3391 do_numbering(natoms, groups, accelerationGroupNames, defaultIndexGroups, gnames,
3392 SimulationAtomGroupType::Acceleration,
3393 restnm, egrptpALL_GENREST, bVerbose, wi);
3394 nr = groups->groups[SimulationAtomGroupType::Acceleration].size();
3395 snew(ir->opts.acc, nr);
3396 ir->opts.ngacc = nr;
3398 convertRvecs(wi, accelerations, "anneal-time", ir->opts.acc);
3400 auto freezeDims = gmx::splitString(is->frdim);
3401 auto freezeGroupNames = gmx::splitString(is->freeze);
3402 if (freezeDims.size() != DIM * freezeGroupNames.size())
3404 gmx_fatal(FARGS, "Invalid Freezing input: %zu groups and %zu freeze values",
3405 freezeGroupNames.size(), freezeDims.size());
3407 do_numbering(natoms, groups, freezeGroupNames, defaultIndexGroups, gnames,
3408 SimulationAtomGroupType::Freeze,
3409 restnm, egrptpALL_GENREST, bVerbose, wi);
3410 nr = groups->groups[SimulationAtomGroupType::Freeze].size();
3411 ir->opts.ngfrz = nr;
3412 snew(ir->opts.nFreeze, nr);
3413 for (i = k = 0; (size_t(i) < freezeGroupNames.size()); i++)
3415 for (j = 0; (j < DIM); j++, k++)
3417 ir->opts.nFreeze[i][j] = static_cast<int>(gmx::equalCaseInsensitive(freezeDims[k], "Y", 1));
3418 if (!ir->opts.nFreeze[i][j])
3420 if (!gmx::equalCaseInsensitive(freezeDims[k], "N", 1))
3422 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3423 "(not %s)", freezeDims[k].c_str());
3424 warning(wi, warn_buf);
3429 for (; (i < nr); i++)
3431 for (j = 0; (j < DIM); j++)
3433 ir->opts.nFreeze[i][j] = 0;
3437 auto energyGroupNames = gmx::splitString(is->energy);
3438 do_numbering(natoms, groups, energyGroupNames, defaultIndexGroups, gnames,
3439 SimulationAtomGroupType::EnergyOutput,
3440 restnm, egrptpALL_GENREST, bVerbose, wi);
3441 add_wall_energrps(groups, ir->nwall, symtab);
3442 ir->opts.ngener = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3443 auto vcmGroupNames = gmx::splitString(is->vcm);
3445 do_numbering(natoms, groups, vcmGroupNames, defaultIndexGroups, gnames,
3446 SimulationAtomGroupType::MassCenterVelocityRemoval,
3447 restnm, vcmGroupNames.empty() ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3450 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3451 "This may lead to artifacts.\n"
3452 "In most cases one should use one group for the whole system.");
3455 /* Now we have filled the freeze struct, so we can calculate NRDF */
3456 calc_nrdf(mtop, ir, gnames);
3458 auto user1GroupNames = gmx::splitString(is->user1);
3459 do_numbering(natoms, groups, user1GroupNames, defaultIndexGroups, gnames,
3460 SimulationAtomGroupType::User1,
3461 restnm, egrptpALL_GENREST, bVerbose, wi);
3462 auto user2GroupNames = gmx::splitString(is->user2);
3463 do_numbering(natoms, groups, user2GroupNames, defaultIndexGroups, gnames,
3464 SimulationAtomGroupType::User2,
3465 restnm, egrptpALL_GENREST, bVerbose, wi);
3466 auto compressedXGroupNames = gmx::splitString(is->x_compressed_groups);
3467 do_numbering(natoms, groups, compressedXGroupNames, defaultIndexGroups, gnames,
3468 SimulationAtomGroupType::CompressedPositionOutput,
3469 restnm, egrptpONE, bVerbose, wi);
3470 auto orirefFitGroupNames = gmx::splitString(is->orirefitgrp);
3471 do_numbering(natoms, groups, orirefFitGroupNames, defaultIndexGroups, gnames,
3472 SimulationAtomGroupType::OrientationRestraintsFit,
3473 restnm, egrptpALL_GENREST, bVerbose, wi);
3475 /* QMMM input processing */
3476 auto qmGroupNames = gmx::splitString(is->QMMM);
3477 auto qmMethods = gmx::splitString(is->QMmethod);
3478 auto qmBasisSets = gmx::splitString(is->QMbasis);
3479 if (ir->eI != eiMimic)
3481 if (qmMethods.size() != qmGroupNames.size() ||
3482 qmBasisSets.size() != qmGroupNames.size())
3484 gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets"
3485 " and %zu methods\n", qmGroupNames.size(),
3486 qmBasisSets.size(), qmMethods.size());
3488 /* group rest, if any, is always MM! */
3489 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3490 SimulationAtomGroupType::QuantumMechanics,
3491 restnm, egrptpALL_GENREST, bVerbose, wi);
3492 nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
3493 ir->opts.ngQM = qmGroupNames.size();
3494 snew(ir->opts.QMmethod, nr);
3495 snew(ir->opts.QMbasis, nr);
3496 for (i = 0; i < nr; i++)
3498 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3499 * converted to the corresponding enum in names.c
3501 ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(),
3504 ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(),
3509 auto qmMultiplicities = gmx::splitString(is->QMmult);
3510 auto qmCharges = gmx::splitString(is->QMcharge);
3511 auto qmbSH = gmx::splitString(is->bSH);
3512 snew(ir->opts.QMmult, nr);
3513 snew(ir->opts.QMcharge, nr);
3514 snew(ir->opts.bSH, nr);
3515 convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult);
3516 convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge);
3517 convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH);
3519 auto CASelectrons = gmx::splitString(is->CASelectrons);
3520 auto CASorbitals = gmx::splitString(is->CASorbitals);
3521 snew(ir->opts.CASelectrons, nr);
3522 snew(ir->opts.CASorbitals, nr);
3523 convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons);
3524 convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals);
3526 auto SAon = gmx::splitString(is->SAon);
3527 auto SAoff = gmx::splitString(is->SAoff);
3528 auto SAsteps = gmx::splitString(is->SAsteps);
3529 snew(ir->opts.SAon, nr);
3530 snew(ir->opts.SAoff, nr);
3531 snew(ir->opts.SAsteps, nr);
3532 convertInts(wi, SAon, "SAon", ir->opts.SAon);
3533 convertInts(wi, SAoff, "SAoff", ir->opts.SAoff);
3534 convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps);
3539 if (qmGroupNames.size() > 1)
3541 gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
3543 /* group rest, if any, is always MM! */
3544 do_numbering(natoms, groups, qmGroupNames, defaultIndexGroups, gnames,
3545 SimulationAtomGroupType::QuantumMechanics,
3546 restnm, egrptpALL_GENREST, bVerbose, wi);
3548 ir->opts.ngQM = qmGroupNames.size();
3551 /* end of QMMM input */
3555 for (auto group : gmx::keysOf(groups->groups))
3557 fprintf(stderr, "%-16s has %zu element(s):", shortName(group), groups->groups[group].size());
3558 for (const auto &entry : groups->groups[group])
3560 fprintf(stderr, " %s", *(groups->groupNames[entry]));
3562 fprintf(stderr, "\n");
3566 nr = groups->groups[SimulationAtomGroupType::EnergyOutput].size();
3567 snew(ir->opts.egp_flags, nr*nr);
3569 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3570 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3572 warning_error(wi, "Energy group exclusions are currently not supported");
3574 if (bExcl && EEL_FULL(ir->coulombtype))
3576 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3579 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3580 if (bTable && !(ir->vdwtype == evdwUSER) &&
3581 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3582 !(ir->coulombtype == eelPMEUSERSWITCH))
3584 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3587 /* final check before going out of scope if simulated tempering variables
3588 * need to be set to default values.
3590 if ((ir->expandedvals->nstexpanded < 0) && ir->bSimTemp)
3592 ir->expandedvals->nstexpanded = 2*static_cast<int>(ir->opts.tau_t[0]/ir->delta_t);
3593 warning(wi, gmx::formatString("the value for nstexpanded was not specified for "
3594 " expanded ensemble simulated tempering. It is set to 2*tau_t (%d) "
3595 "by default, but it is recommended to set it to an explicit value!",
3596 ir->expandedvals->nstexpanded));
3598 for (i = 0; (i < defaultIndexGroups->nr); i++)
3603 done_blocka(defaultIndexGroups);
3604 sfree(defaultIndexGroups);
3610 static void check_disre(const gmx_mtop_t *mtop)
3612 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3614 const gmx_ffparams_t &ffparams = mtop->ffparams;
3617 for (int i = 0; i < ffparams.numTypes(); i++)
3619 int ftype = ffparams.functype[i];
3620 if (ftype == F_DISRES)
3622 int label = ffparams.iparams[i].disres.label;
3623 if (label == old_label)
3625 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3633 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3634 "probably the parameters for multiple pairs in one restraint "
3635 "are not identical\n", ndouble);
3640 static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3645 gmx_mtop_ilistloop_t iloop;
3654 for (d = 0; d < DIM; d++)
3656 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3658 /* Check for freeze groups */
3659 for (g = 0; g < ir->opts.ngfrz; g++)
3661 for (d = 0; d < DIM; d++)
3663 if (ir->opts.nFreeze[g][d] != 0)
3671 /* Check for position restraints */
3672 iloop = gmx_mtop_ilistloop_init(sys);
3673 while (const InteractionLists *ilist = gmx_mtop_ilistloop_next(iloop, &nmol))
3676 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3678 for (i = 0; i < (*ilist)[F_POSRES].size(); i += 2)
3680 pr = &sys->ffparams.iparams[(*ilist)[F_POSRES].iatoms[i]];
3681 for (d = 0; d < DIM; d++)
3683 if (pr->posres.fcA[d] != 0)
3689 for (i = 0; i < (*ilist)[F_FBPOSRES].size(); i += 2)
3691 /* Check for flat-bottom posres */
3692 pr = &sys->ffparams.iparams[(*ilist)[F_FBPOSRES].iatoms[i]];
3693 if (pr->fbposres.k != 0)
3695 switch (pr->fbposres.geom)
3697 case efbposresSPHERE:
3698 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3700 case efbposresCYLINDERX:
3701 AbsRef[YY] = AbsRef[ZZ] = 1;
3703 case efbposresCYLINDERY:
3704 AbsRef[XX] = AbsRef[ZZ] = 1;
3706 case efbposresCYLINDER:
3707 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3708 case efbposresCYLINDERZ:
3709 AbsRef[XX] = AbsRef[YY] = 1;
3711 case efbposresX: /* d=XX */
3712 case efbposresY: /* d=YY */
3713 case efbposresZ: /* d=ZZ */
3714 d = pr->fbposres.geom - efbposresX;
3718 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3719 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3727 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3731 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3732 bool *bC6ParametersWorkWithGeometricRules,
3733 bool *bC6ParametersWorkWithLBRules,
3734 bool *bLBRulesPossible)
3736 int ntypes, tpi, tpj;
3739 double c6i, c6j, c12i, c12j;
3740 double c6, c6_geometric, c6_LB;
3741 double sigmai, sigmaj, epsi, epsj;
3742 bool bCanDoLBRules, bCanDoGeometricRules;
3745 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3746 * force-field floating point parameters.
3749 ptr = getenv("GMX_LJCOMB_TOL");
3753 double gmx_unused canary;
3755 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3757 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
3762 *bC6ParametersWorkWithLBRules = TRUE;
3763 *bC6ParametersWorkWithGeometricRules = TRUE;
3764 bCanDoLBRules = TRUE;
3765 ntypes = mtop->ffparams.atnr;
3766 snew(typecount, ntypes);
3767 gmx_mtop_count_atomtypes(mtop, state, typecount);
3768 *bLBRulesPossible = TRUE;
3769 for (tpi = 0; tpi < ntypes; ++tpi)
3771 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3772 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3773 for (tpj = tpi; tpj < ntypes; ++tpj)
3775 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3776 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3777 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3778 c6_geometric = std::sqrt(c6i * c6j);
3779 if (!gmx_numzero(c6_geometric))
3781 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3783 sigmai = gmx::sixthroot(c12i / c6i);
3784 sigmaj = gmx::sixthroot(c12j / c6j);
3785 epsi = c6i * c6i /(4.0 * c12i);
3786 epsj = c6j * c6j /(4.0 * c12j);
3787 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
3791 *bLBRulesPossible = FALSE;
3792 c6_LB = c6_geometric;
3794 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3799 *bC6ParametersWorkWithLBRules = FALSE;
3802 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3804 if (!bCanDoGeometricRules)
3806 *bC6ParametersWorkWithGeometricRules = FALSE;
3814 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3817 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3819 check_combination_rule_differences(mtop, 0,
3820 &bC6ParametersWorkWithGeometricRules,
3821 &bC6ParametersWorkWithLBRules,
3823 if (ir->ljpme_combination_rule == eljpmeLB)
3825 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
3827 warning(wi, "You are using arithmetic-geometric combination rules "
3828 "in LJ-PME, but your non-bonded C6 parameters do not "
3829 "follow these rules.");
3834 if (!bC6ParametersWorkWithGeometricRules)
3836 if (ir->eDispCorr != edispcNO)
3838 warning_note(wi, "You are using geometric combination rules in "
3839 "LJ-PME, but your non-bonded C6 parameters do "
3840 "not follow these rules. "
3841 "This will introduce very small errors in the forces and energies in "
3842 "your simulations. Dispersion correction will correct total energy "
3843 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3847 warning_note(wi, "You are using geometric combination rules in "
3848 "LJ-PME, but your non-bonded C6 parameters do "
3849 "not follow these rules. "
3850 "This will introduce very small errors in the forces and energies in "
3851 "your simulations. If your system is homogeneous, consider using dispersion correction "
3852 "for the total energy and pressure.");
3858 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3861 char err_buf[STRLEN];
3866 gmx_mtop_atomloop_block_t aloopb;
3868 char warn_buf[STRLEN];
3870 set_warning_line(wi, mdparin, -1);
3872 if (ir->cutoff_scheme == ecutsVERLET &&
3873 ir->verletbuf_tol > 0 &&
3875 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
3876 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
3878 /* Check if a too small Verlet buffer might potentially
3879 * cause more drift than the thermostat can couple off.
3881 /* Temperature error fraction for warning and suggestion */
3882 const real T_error_warn = 0.002;
3883 const real T_error_suggest = 0.001;
3884 /* For safety: 2 DOF per atom (typical with constraints) */
3885 const real nrdf_at = 2;
3886 real T, tau, max_T_error;
3891 for (i = 0; i < ir->opts.ngtc; i++)
3893 T = std::max(T, ir->opts.ref_t[i]);
3894 tau = std::max(tau, ir->opts.tau_t[i]);
3898 /* This is a worst case estimate of the temperature error,
3899 * assuming perfect buffer estimation and no cancelation
3900 * of errors. The factor 0.5 is because energy distributes
3901 * equally over Ekin and Epot.
3903 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
3904 if (max_T_error > T_error_warn)
3906 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.",
3907 ir->verletbuf_tol, T, tau,
3909 100*T_error_suggest,
3910 ir->verletbuf_tol*T_error_suggest/max_T_error);
3911 warning(wi, warn_buf);
3916 if (ETC_ANDERSEN(ir->etc))
3920 for (i = 0; i < ir->opts.ngtc; i++)
3922 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
3923 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
3924 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
3925 i, ir->opts.tau_t[i]);
3926 CHECK(ir->opts.tau_t[i] < 0);
3929 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
3931 for (i = 0; i < ir->opts.ngtc; i++)
3933 int nsteps = gmx::roundToInt(ir->opts.tau_t[i]/ir->delta_t);
3934 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);
3935 CHECK(nsteps % ir->nstcomm != 0);
3940 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3941 ir->comm_mode == ecmNO &&
3942 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
3943 !ETC_ANDERSEN(ir->etc))
3945 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");
3948 if (ir->epc == epcPARRINELLORAHMAN &&
3949 ir->etc == etcNOSEHOOVER)
3952 for (int g = 0; g < ir->opts.ngtc; g++)
3954 tau_t_max = std::max(tau_t_max, ir->opts.tau_t[g]);
3956 if (ir->tau_p < 1.9*tau_t_max)
3958 std::string message =
3959 gmx::formatString("With %s T-coupling and %s p-coupling, "
3960 "%s (%g) should be at least twice as large as %s (%g) to avoid resonances",
3961 etcoupl_names[ir->etc],
3962 epcoupl_names[ir->epc],
3964 "tau-t", tau_t_max);
3965 warning(wi, message.c_str());
3969 /* Check for pressure coupling with absolute position restraints */
3970 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3972 absolute_reference(ir, sys, TRUE, AbsRef);
3974 for (m = 0; m < DIM; m++)
3976 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3978 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3986 aloopb = gmx_mtop_atomloop_block_init(sys);
3988 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3990 if (atom->q != 0 || atom->qB != 0)
3998 if (EEL_FULL(ir->coulombtype))
4001 "You are using full electrostatics treatment %s for a system without charges.\n"
4002 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4003 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4004 warning(wi, err_buf);
4009 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4012 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4013 "You might want to consider using %s electrostatics.\n",
4015 warning_note(wi, err_buf);
4019 /* Check if combination rules used in LJ-PME are the same as in the force field */
4020 if (EVDW_PME(ir->vdwtype))
4022 check_combination_rules(ir, sys, wi);
4025 /* Generalized reaction field */
4026 if (ir->coulombtype == eelGRF_NOTUSED)
4028 warning_error(wi, "Generalized reaction-field electrostatics is no longer supported. "
4029 "You can use normal reaction-field instead and compute the reaction-field constant by hand.");
4033 for (int i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4035 for (m = 0; (m < DIM); m++)
4037 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4046 snew(mgrp, sys->groups.groups[SimulationAtomGroupType::Acceleration].size());
4047 for (const AtomProxy atomP : AtomRange(*sys))
4049 const t_atom &local = atomP.atom();
4050 int i = atomP.globalAtomNumber();
4051 mgrp[getGroupType(sys->groups, SimulationAtomGroupType::Acceleration, i)] += local.m;
4054 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4056 for (m = 0; (m < DIM); m++)
4058 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4062 for (m = 0; (m < DIM); m++)
4064 if (fabs(acc[m]) > 1e-6)
4066 const char *dim[DIM] = { "X", "Y", "Z" };
4068 "Net Acceleration in %s direction, will %s be corrected\n",
4069 dim[m], ir->nstcomm != 0 ? "" : "not");
4070 if (ir->nstcomm != 0 && m < ndof_com(ir))
4073 for (i = 0; (i < gmx::ssize(sys->groups.groups[SimulationAtomGroupType::Acceleration])); i++)
4075 ir->opts.acc[i][m] -= acc[m];
4083 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4084 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4086 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4094 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4096 if (ir->pull->coord[i].group[0] == 0 ||
4097 ir->pull->coord[i].group[1] == 0)
4099 absolute_reference(ir, sys, FALSE, AbsRef);
4100 for (m = 0; m < DIM; m++)
4102 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4104 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.");
4112 for (i = 0; i < 3; i++)
4114 for (m = 0; m <= i; m++)
4116 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4117 ir->deform[i][m] != 0)
4119 for (c = 0; c < ir->pull->ncoord; c++)
4121 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4122 ir->pull->coord[c].vec[m] != 0)
4124 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4135 void double_check(t_inputrec *ir, matrix box,
4136 bool bHasNormalConstraints,
4137 bool bHasAnyConstraints,
4140 char warn_buf[STRLEN];
4143 ptr = check_box(ir->ePBC, box);
4146 warning_error(wi, ptr);
4149 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4151 if (ir->shake_tol <= 0.0)
4153 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4155 warning_error(wi, warn_buf);
4159 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4161 /* If we have Lincs constraints: */
4162 if (ir->eI == eiMD && ir->etc == etcNO &&
4163 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4165 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4166 warning_note(wi, warn_buf);
4169 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4171 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4172 warning_note(wi, warn_buf);
4174 if (ir->epc == epcMTTK)
4176 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4180 if (bHasAnyConstraints && ir->epc == epcMTTK)
4182 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4185 if (ir->LincsWarnAngle > 90.0)
4187 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4188 warning(wi, warn_buf);
4189 ir->LincsWarnAngle = 90.0;
4192 if (ir->ePBC != epbcNONE)
4194 if (ir->nstlist == 0)
4196 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4198 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4200 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");
4201 warning_error(wi, warn_buf);