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, 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.
50 #include "gromacs/fileio/readinp.h"
51 #include "gromacs/fileio/warninp.h"
52 #include "gromacs/gmxlib/chargegroup.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/gmxpreprocess/toputil.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/calc_verletbuf.h"
59 #include "gromacs/mdrunutility/mdmodules.h"
60 #include "gromacs/mdtypes/inputrec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/mdtypes/pull-params.h"
63 #include "gromacs/options/options.h"
64 #include "gromacs/options/treesupport.h"
65 #include "gromacs/pbcutil/pbc.h"
66 #include "gromacs/topology/block.h"
67 #include "gromacs/topology/ifunc.h"
68 #include "gromacs/topology/index.h"
69 #include "gromacs/topology/mtop_util.h"
70 #include "gromacs/topology/symtab.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/utility/cstringutil.h"
73 #include "gromacs/utility/exceptions.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/filestream.h"
76 #include "gromacs/utility/gmxassert.h"
77 #include "gromacs/utility/ikeyvaluetreeerror.h"
78 #include "gromacs/utility/keyvaluetree.h"
79 #include "gromacs/utility/keyvaluetreetransform.h"
80 #include "gromacs/utility/smalloc.h"
81 #include "gromacs/utility/stringcompare.h"
82 #include "gromacs/utility/stringutil.h"
87 /* Resource parameters
88 * Do not change any of these until you read the instruction
89 * in readinp.h. Some cpp's do not take spaces after the backslash
90 * (like the c-shell), which will give you a very weird compiler
94 typedef struct t_inputrec_strings
96 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
97 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
98 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
99 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
100 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
102 char fep_lambda[efptNR][STRLEN];
103 char lambda_weights[STRLEN];
106 char anneal[STRLEN], anneal_npoints[STRLEN],
107 anneal_time[STRLEN], anneal_temp[STRLEN];
108 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
109 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
110 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
112 } gmx_inputrec_strings;
114 static gmx_inputrec_strings *is = nullptr;
116 void init_inputrec_strings()
120 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.");
125 void done_inputrec_strings()
133 egrptpALL, /* All particles have to be a member of a group. */
134 egrptpALL_GENREST, /* A rest group with name is generated for particles *
135 * that are not part of any group. */
136 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
137 * for the rest group. */
138 egrptpONE /* Merge all selected groups into one group, *
139 * make a rest group for the remaining particles. */
142 static const char *constraints[eshNR+1] = {
143 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
146 static const char *couple_lam[ecouplamNR+1] = {
147 "vdw-q", "vdw", "q", "none", nullptr
150 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
155 for (i = 0; i < ntemps; i++)
157 /* simple linear scaling -- allows more control */
158 if (simtemp->eSimTempScale == esimtempLINEAR)
160 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
162 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
164 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
166 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
168 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
173 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
174 gmx_fatal(FARGS, errorstr);
181 static void _low_check(gmx_bool b, const char *s, warninp_t wi)
185 warning_error(wi, s);
189 static void check_nst(const char *desc_nst, int nst,
190 const char *desc_p, int *p,
195 if (*p > 0 && *p % nst != 0)
197 /* Round up to the next multiple of nst */
198 *p = ((*p)/nst + 1)*nst;
199 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
200 desc_p, desc_nst, desc_p, *p);
205 static gmx_bool ir_NVE(const t_inputrec *ir)
207 return (EI_MD(ir->eI) && ir->etc == etcNO);
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 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
228 if (*eintmod == eintmodPOTSHIFT_VERLET)
230 if (ir->cutoff_scheme == ecutsVERLET)
232 *eintmod = eintmodPOTSHIFT;
236 *eintmod = eintmodNONE;
241 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
243 /* Check internal consistency.
244 * NOTE: index groups are not set here yet, don't check things
245 * like temperature coupling group options here, but in triple_check
248 /* Strange macro: first one fills the err_buf, and then one can check
249 * the condition, which will print the message and increase the error
252 #define CHECK(b) _low_check(b, err_buf, wi)
253 char err_buf[256], warn_buf[STRLEN];
256 t_lambda *fep = ir->fepvals;
257 t_expanded *expand = ir->expandedvals;
259 set_warning_line(wi, mdparin, -1);
261 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
263 sprintf(warn_buf, "%s electrostatics is no longer supported",
264 eel_names[eelRF_NEC_UNSUPPORTED]);
265 warning_error(wi, warn_buf);
268 /* BASIC CUT-OFF STUFF */
269 if (ir->rcoulomb < 0)
271 warning_error(wi, "rcoulomb should be >= 0");
275 warning_error(wi, "rvdw should be >= 0");
278 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
280 warning_error(wi, "rlist should be >= 0");
282 sprintf(err_buf, "nstlist can not be smaller than 0. (If you were trying to use the heuristic neighbour-list update scheme for efficient buffering for improved energy conservation, please use the Verlet cut-off scheme instead.)");
283 CHECK(ir->nstlist < 0);
285 process_interaction_modifier(ir, &ir->coulomb_modifier);
286 process_interaction_modifier(ir, &ir->vdw_modifier);
288 if (ir->cutoff_scheme == ecutsGROUP)
291 "The group cutoff scheme is deprecated since GROMACS 5.0 and will be removed in a future "
292 "release when all interaction forms are supported for the verlet scheme. The verlet "
293 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
295 if (ir->rlist > 0 && ir->rlist < ir->rcoulomb)
297 gmx_fatal(FARGS, "rcoulomb must not be greater than rlist (twin-range schemes are not supported)");
299 if (ir->rlist > 0 && ir->rlist < ir->rvdw)
301 gmx_fatal(FARGS, "rvdw must not be greater than rlist (twin-range schemes are not supported)");
304 if (ir->rlist == 0 && ir->ePBC != epbcNONE)
306 warning_error(wi, "Can not have an infinite cut-off with PBC");
310 if (ir->cutoff_scheme == ecutsVERLET)
314 /* Normal Verlet type neighbor-list, currently only limited feature support */
315 if (inputrec2nboundeddim(ir) < 3)
317 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
320 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
321 if (ir->rcoulomb != ir->rvdw)
323 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
324 // for PME load balancing, we can support this exception.
325 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
326 ir->vdwtype == evdwCUT &&
327 ir->rcoulomb > ir->rvdw);
328 if (!bUsesPmeTwinRangeKernel)
330 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
334 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
336 if (ir->vdw_modifier == eintmodNONE ||
337 ir->vdw_modifier == eintmodPOTSHIFT)
339 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
341 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]);
342 warning_note(wi, warn_buf);
344 ir->vdwtype = evdwCUT;
348 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
349 warning_error(wi, warn_buf);
353 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
355 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
357 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
358 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
360 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
362 if (!(ir->coulomb_modifier == eintmodNONE ||
363 ir->coulomb_modifier == eintmodPOTSHIFT))
365 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
366 warning_error(wi, warn_buf);
369 if (ir->implicit_solvent != eisNO)
371 warning_error(wi, "Implicit solvent is not (yet) supported with the with Verlet lists.");
374 if (ir->nstlist <= 0)
376 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
379 if (ir->nstlist < 10)
381 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.");
384 rc_max = std::max(ir->rvdw, ir->rcoulomb);
386 if (ir->verletbuf_tol <= 0)
388 if (ir->verletbuf_tol == 0)
390 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
393 if (ir->rlist < rc_max)
395 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
398 if (ir->rlist == rc_max && ir->nstlist > 1)
400 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.");
405 if (ir->rlist > rc_max)
407 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.");
410 if (ir->nstlist == 1)
412 /* No buffer required */
417 if (EI_DYNAMICS(ir->eI))
419 if (inputrec2nboundeddim(ir) < 3)
421 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.");
423 /* Set rlist temporarily so we can continue processing */
428 /* Set the buffer to 5% of the cut-off */
429 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
435 /* GENERAL INTEGRATOR STUFF */
438 if (ir->etc != etcNO)
440 if (EI_RANDOM(ir->eI))
442 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]);
446 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
448 warning_note(wi, warn_buf);
452 if (ir->eI == eiVVAK)
454 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]);
455 warning_note(wi, warn_buf);
457 if (!EI_DYNAMICS(ir->eI))
459 if (ir->epc != epcNO)
461 sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
462 warning_note(wi, warn_buf);
466 if (EI_DYNAMICS(ir->eI))
468 if (ir->nstcalcenergy < 0)
470 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
471 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
473 /* nstcalcenergy larger than nstener does not make sense.
474 * We ideally want nstcalcenergy=nstener.
478 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
482 ir->nstcalcenergy = ir->nstenergy;
486 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
487 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
488 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
491 const char *nsten = "nstenergy";
492 const char *nstdh = "nstdhdl";
493 const char *min_name = nsten;
494 int min_nst = ir->nstenergy;
496 /* find the smallest of ( nstenergy, nstdhdl ) */
497 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
498 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
500 min_nst = ir->fepvals->nstdhdl;
503 /* If the user sets nstenergy small, we should respect that */
505 "Setting nstcalcenergy (%d) equal to %s (%d)",
506 ir->nstcalcenergy, min_name, min_nst);
507 warning_note(wi, warn_buf);
508 ir->nstcalcenergy = min_nst;
511 if (ir->epc != epcNO)
513 if (ir->nstpcouple < 0)
515 ir->nstpcouple = ir_optimal_nstpcouple(ir);
519 if (ir->nstcalcenergy > 0)
521 if (ir->efep != efepNO)
523 /* nstdhdl should be a multiple of nstcalcenergy */
524 check_nst("nstcalcenergy", ir->nstcalcenergy,
525 "nstdhdl", &ir->fepvals->nstdhdl, wi);
526 /* nstexpanded should be a multiple of nstcalcenergy */
527 check_nst("nstcalcenergy", ir->nstcalcenergy,
528 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
530 /* for storing exact averages nstenergy should be
531 * a multiple of nstcalcenergy
533 check_nst("nstcalcenergy", ir->nstcalcenergy,
534 "nstenergy", &ir->nstenergy, wi);
538 if (ir->nsteps == 0 && !ir->bContinuation)
540 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
544 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
545 ir->bContinuation && ir->ld_seed != -1)
547 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)");
553 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
554 CHECK(ir->ePBC != epbcXYZ);
555 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
556 CHECK(ir->ns_type != ensGRID);
557 sprintf(err_buf, "with TPI nstlist should be larger than zero");
558 CHECK(ir->nstlist <= 0);
559 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
560 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
561 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
562 CHECK(ir->cutoff_scheme == ecutsVERLET);
566 if ( (opts->nshake > 0) && (opts->bMorse) )
569 "Using morse bond-potentials while constraining bonds is useless");
570 warning(wi, warn_buf);
573 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
574 ir->bContinuation && ir->ld_seed != -1)
576 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)");
578 /* verify simulated tempering options */
582 gmx_bool bAllTempZero = TRUE;
583 for (i = 0; i < fep->n_lambda; i++)
585 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]);
586 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
587 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
589 bAllTempZero = FALSE;
592 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
593 CHECK(bAllTempZero == TRUE);
595 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
596 CHECK(ir->eI != eiVV);
598 /* check compatability of the temperature coupling with simulated tempering */
600 if (ir->etc == etcNOSEHOOVER)
602 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
603 warning_note(wi, warn_buf);
606 /* check that the temperatures make sense */
608 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);
609 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
611 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
612 CHECK(ir->simtempvals->simtemp_high <= 0);
614 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
615 CHECK(ir->simtempvals->simtemp_low <= 0);
618 /* verify free energy options */
620 if (ir->efep != efepNO)
623 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
625 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
627 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
628 (int)fep->sc_r_power);
629 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
631 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
632 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
634 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
635 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
637 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
638 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
640 sprintf(err_buf, "Free-energy not implemented for Ewald");
641 CHECK(ir->coulombtype == eelEWALD);
643 /* check validty of lambda inputs */
644 if (fep->n_lambda == 0)
646 /* Clear output in case of no states:*/
647 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
648 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
652 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
653 CHECK((fep->init_fep_state >= fep->n_lambda));
656 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
657 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
659 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",
660 fep->init_lambda, fep->init_fep_state);
661 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
665 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
669 for (i = 0; i < efptNR; i++)
671 if (fep->separate_dvdl[i])
676 if (n_lambda_terms > 1)
678 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.");
679 warning(wi, warn_buf);
682 if (n_lambda_terms < 2 && fep->n_lambda > 0)
685 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
689 for (j = 0; j < efptNR; j++)
691 for (i = 0; i < fep->n_lambda; i++)
693 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]);
694 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
698 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
700 for (i = 0; i < fep->n_lambda; i++)
702 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],
703 fep->all_lambda[efptCOUL][i]);
704 CHECK((fep->sc_alpha > 0) &&
705 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
706 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
707 ((fep->all_lambda[efptVDW][i] > 0.0) &&
708 (fep->all_lambda[efptVDW][i] < 1.0))));
712 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
714 real sigma, lambda, r_sc;
717 /* Maximum estimate for A and B charges equal with lambda power 1 */
719 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);
720 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.",
722 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
723 warning_note(wi, warn_buf);
726 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
727 be treated differently, but that's the next step */
729 for (i = 0; i < efptNR; i++)
731 for (j = 0; j < fep->n_lambda; j++)
733 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
734 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
739 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
743 /* checking equilibration of weights inputs for validity */
745 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
746 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
747 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
749 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
750 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
751 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
753 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
754 expand->equil_steps, elmceq_names[elmceqSTEPS]);
755 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
757 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
758 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
759 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
761 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
762 expand->equil_ratio, elmceq_names[elmceqRATIO]);
763 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
765 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
766 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
767 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
769 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
770 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
771 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
773 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
774 expand->equil_steps, elmceq_names[elmceqSTEPS]);
775 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
777 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
778 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
779 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
781 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
782 expand->equil_ratio, elmceq_names[elmceqRATIO]);
783 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
785 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
786 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
787 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
789 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
790 CHECK((expand->lmc_repeats <= 0));
791 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
792 CHECK((expand->minvarmin <= 0));
793 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
794 CHECK((expand->c_range < 0));
795 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
796 fep->init_fep_state, expand->lmc_forced_nstart);
797 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
798 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
799 CHECK((expand->lmc_forced_nstart < 0));
800 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
801 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
803 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
804 CHECK((expand->init_wl_delta < 0));
805 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
806 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
807 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
808 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
810 /* if there is no temperature control, we need to specify an MC temperature */
811 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
812 if (expand->nstTij > 0)
814 sprintf(err_buf, "nstlog must be non-zero");
815 CHECK(ir->nstlog == 0);
816 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
817 expand->nstTij, ir->nstlog);
818 CHECK((expand->nstTij % ir->nstlog) != 0);
823 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
824 CHECK(ir->nwall && ir->ePBC != epbcXY);
827 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
829 if (ir->ePBC == epbcNONE)
831 if (ir->epc != epcNO)
833 warning(wi, "Turning off pressure coupling for vacuum system");
839 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
840 epbc_names[ir->ePBC]);
841 CHECK(ir->epc != epcNO);
843 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
844 CHECK(EEL_FULL(ir->coulombtype));
846 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
847 epbc_names[ir->ePBC]);
848 CHECK(ir->eDispCorr != edispcNO);
851 if (ir->rlist == 0.0)
853 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
854 "with coulombtype = %s or coulombtype = %s\n"
855 "without periodic boundary conditions (pbc = %s) and\n"
856 "rcoulomb and rvdw set to zero",
857 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
858 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
859 (ir->ePBC != epbcNONE) ||
860 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
864 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
869 if (ir->nstcomm == 0)
871 ir->comm_mode = ecmNO;
873 if (ir->comm_mode != ecmNO)
877 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");
878 ir->nstcomm = abs(ir->nstcomm);
881 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
883 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
884 ir->nstcomm = ir->nstcalcenergy;
887 if (ir->comm_mode == ecmANGULAR)
889 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
890 CHECK(ir->bPeriodicMols);
891 if (ir->ePBC != epbcNONE)
893 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.");
898 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
900 warning_note(wi, "Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
903 /* TEMPERATURE COUPLING */
904 if (ir->etc == etcYES)
906 ir->etc = etcBERENDSEN;
907 warning_note(wi, "Old option for temperature coupling given: "
908 "changing \"yes\" to \"Berendsen\"\n");
911 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
913 if (ir->opts.nhchainlength < 1)
915 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
916 ir->opts.nhchainlength = 1;
917 warning(wi, warn_buf);
920 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
922 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
923 ir->opts.nhchainlength = 1;
928 ir->opts.nhchainlength = 0;
931 if (ir->eI == eiVVAK)
933 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
935 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
938 if (ETC_ANDERSEN(ir->etc))
940 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
941 CHECK(!(EI_VV(ir->eI)));
943 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
945 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]);
946 warning_note(wi, warn_buf);
949 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]);
950 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
953 if (ir->etc == etcBERENDSEN)
955 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
956 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
957 warning_note(wi, warn_buf);
960 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
961 && ir->epc == epcBERENDSEN)
963 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
964 "true ensemble for the thermostat");
965 warning(wi, warn_buf);
968 /* PRESSURE COUPLING */
969 if (ir->epc == epcISOTROPIC)
971 ir->epc = epcBERENDSEN;
972 warning_note(wi, "Old option for pressure coupling given: "
973 "changing \"Isotropic\" to \"Berendsen\"\n");
976 if (ir->epc != epcNO)
978 dt_pcoupl = ir->nstpcouple*ir->delta_t;
980 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
981 CHECK(ir->tau_p <= 0);
983 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
985 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
986 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
987 warning(wi, warn_buf);
990 sprintf(err_buf, "compressibility must be > 0 when using pressure"
991 " coupling %s\n", EPCOUPLTYPE(ir->epc));
992 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
993 ir->compress[ZZ][ZZ] < 0 ||
994 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
995 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
997 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1000 "You are generating velocities so I am assuming you "
1001 "are equilibrating a system. You are using "
1002 "%s pressure coupling, but this can be "
1003 "unstable for equilibration. If your system crashes, try "
1004 "equilibrating first with Berendsen pressure coupling. If "
1005 "you are not equilibrating the system, you can probably "
1006 "ignore this warning.",
1007 epcoupl_names[ir->epc]);
1008 warning(wi, warn_buf);
1014 if (ir->epc > epcNO)
1016 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1018 warning_error(wi, "for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
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 (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1044 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1045 warning_note(wi, warn_buf);
1048 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1050 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);
1051 warning(wi, warn_buf);
1052 ir->epsilon_rf = ir->epsilon_r;
1053 ir->epsilon_r = 1.0;
1056 if (ir->epsilon_r == 0)
1059 "It is pointless to use long-range or Generalized Born electrostatics with infinite relative permittivity."
1060 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1061 CHECK(EEL_FULL(ir->coulombtype) || ir->implicit_solvent == eisGBSA);
1064 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1066 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1067 CHECK(ir->epsilon_r < 0);
1070 if (EEL_RF(ir->coulombtype))
1072 /* reaction field (at the cut-off) */
1074 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1076 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1077 eel_names[ir->coulombtype]);
1078 warning(wi, warn_buf);
1079 ir->epsilon_rf = 0.0;
1082 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1083 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1084 (ir->epsilon_r == 0));
1085 if (ir->epsilon_rf == ir->epsilon_r)
1087 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1088 eel_names[ir->coulombtype]);
1089 warning(wi, warn_buf);
1092 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1093 * means the interaction is zero outside rcoulomb, but it helps to
1094 * provide accurate energy conservation.
1096 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1098 if (ir_coulomb_switched(ir))
1101 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1102 eel_names[ir->coulombtype]);
1103 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1106 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1108 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1110 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1111 eel_names[ir->coulombtype]);
1112 CHECK(ir->rlist > ir->rcoulomb);
1116 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1119 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1120 CHECK( ir->coulomb_modifier != eintmodNONE);
1122 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1125 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1126 CHECK( ir->vdw_modifier != eintmodNONE);
1129 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1130 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1133 "The switch/shift interaction settings are just for compatibility; you will get better "
1134 "performance from applying potential modifiers to your interactions!\n");
1135 warning_note(wi, warn_buf);
1138 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1140 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1142 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1143 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.",
1144 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1145 warning(wi, warn_buf);
1149 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1151 if (ir->rvdw_switch == 0)
1153 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.");
1154 warning(wi, warn_buf);
1158 if (EEL_FULL(ir->coulombtype))
1160 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1161 ir->coulombtype == eelPMEUSERSWITCH)
1163 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1164 eel_names[ir->coulombtype]);
1165 CHECK(ir->rcoulomb > ir->rlist);
1167 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1169 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1172 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist.\n"
1173 "For optimal energy conservation,consider using\n"
1174 "a potential modifier.", eel_names[ir->coulombtype]);
1175 CHECK(ir->rcoulomb != ir->rlist);
1180 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1182 // TODO: Move these checks into the ewald module with the options class
1184 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1186 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1188 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
1189 warning_error(wi, warn_buf);
1193 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1195 if (ir->ewald_geometry == eewg3D)
1197 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1198 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1199 warning(wi, warn_buf);
1201 /* This check avoids extra pbc coding for exclusion corrections */
1202 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1203 CHECK(ir->wall_ewald_zfac < 2);
1205 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1206 EEL_FULL(ir->coulombtype))
1208 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1209 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1210 warning(wi, warn_buf);
1212 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1214 if (ir->cutoff_scheme == ecutsVERLET)
1216 sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
1217 eel_names[ir->coulombtype]);
1218 warning(wi, warn_buf);
1222 sprintf(warn_buf, "Dipole corrections to %s electrostatics only work if all charge groups that can cross PBC boundaries are dipoles. If this is not the case set epsilon_surface to 0",
1223 eel_names[ir->coulombtype]);
1224 warning_note(wi, warn_buf);
1228 if (ir_vdw_switched(ir))
1230 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1231 CHECK(ir->rvdw_switch >= ir->rvdw);
1233 if (ir->rvdw_switch < 0.5*ir->rvdw)
1235 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.",
1236 ir->rvdw_switch, ir->rvdw);
1237 warning_note(wi, warn_buf);
1240 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1242 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1244 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1245 CHECK(ir->rlist > ir->rvdw);
1249 if (ir->vdwtype == evdwPME)
1251 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1253 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1254 evdw_names[ir->vdwtype],
1255 eintmod_names[eintmodPOTSHIFT],
1256 eintmod_names[eintmodNONE]);
1260 if (ir->cutoff_scheme == ecutsGROUP)
1262 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1263 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)))
1265 warning_note(wi, "With exact cut-offs, rlist should be "
1266 "larger than rcoulomb and rvdw, so that there "
1267 "is a buffer region for particle motion "
1268 "between neighborsearch steps");
1271 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlist <= ir->rcoulomb)
1273 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rcoulomb.");
1274 warning_note(wi, warn_buf);
1276 if (ir_vdw_switched(ir) && (ir->rlist <= ir->rvdw))
1278 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rvdw.");
1279 warning_note(wi, warn_buf);
1283 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1285 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.");
1288 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1291 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1294 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1296 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1299 /* ENERGY CONSERVATION */
1300 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1302 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1304 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1305 evdw_names[evdwSHIFT]);
1306 warning_note(wi, warn_buf);
1308 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1310 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1311 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1312 warning_note(wi, warn_buf);
1316 /* IMPLICIT SOLVENT */
1317 if (ir->coulombtype == eelGB_NOTUSED)
1319 sprintf(warn_buf, "Invalid option %s for coulombtype",
1320 eel_names[ir->coulombtype]);
1321 warning_error(wi, warn_buf);
1324 if (ir->sa_algorithm == esaSTILL)
1326 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1327 CHECK(ir->sa_algorithm == esaSTILL);
1330 if (ir->implicit_solvent == eisGBSA)
1332 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1333 CHECK(ir->rgbradii != ir->rlist);
1335 if (ir->coulombtype != eelCUT)
1337 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1338 CHECK(ir->coulombtype != eelCUT);
1340 if (ir->vdwtype != evdwCUT)
1342 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1343 CHECK(ir->vdwtype != evdwCUT);
1345 if (ir->nstgbradii < 1)
1347 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1348 warning_note(wi, warn_buf);
1351 if (ir->sa_algorithm == esaNO)
1353 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1354 warning_note(wi, warn_buf);
1356 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1358 sprintf(warn_buf, "Value of sa_surface_tension is < 0. Changing it to 2.05016 or 2.25936 kJ/nm^2/mol for Still and HCT/OBC respectively\n");
1359 warning_note(wi, warn_buf);
1361 if (ir->gb_algorithm == egbSTILL)
1363 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1367 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1370 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1372 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1373 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1380 if (ir->cutoff_scheme != ecutsGROUP)
1382 warning_error(wi, "QMMM is currently only supported with cutoff-scheme=group");
1384 if (!EI_DYNAMICS(ir->eI))
1387 sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
1388 warning_error(wi, buf);
1394 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1398 /* count the number of text elemets separated by whitespace in a string.
1399 str = the input string
1400 maxptr = the maximum number of allowed elements
1401 ptr = the output array of pointers to the first character of each element
1402 returns: the number of elements. */
1403 int str_nelem(const char *str, int maxptr, char *ptr[])
1408 copy0 = gmx_strdup(str);
1411 while (*copy != '\0')
1415 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1423 while ((*copy != '\0') && !isspace(*copy))
1442 /* interpret a number of doubles from a string and put them in an array,
1443 after allocating space for them.
1444 str = the input string
1445 n = the (pre-allocated) number of doubles read
1446 r = the output array of doubles. */
1447 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1452 char warn_buf[STRLEN];
1454 *n = str_nelem(str, MAXPTR, ptr);
1457 for (i = 0; i < *n; i++)
1459 (*r)[i] = strtod(ptr[i], &endptr);
1462 sprintf(warn_buf, "Invalid value %s in string in mdp file. Expected a real number.", ptr[i]);
1463 warning_error(wi, warn_buf);
1468 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1471 int i, j, max_n_lambda, nweights, nfep[efptNR];
1472 t_lambda *fep = ir->fepvals;
1473 t_expanded *expand = ir->expandedvals;
1474 real **count_fep_lambdas;
1475 gmx_bool bOneLambda = TRUE;
1477 snew(count_fep_lambdas, efptNR);
1479 /* FEP input processing */
1480 /* first, identify the number of lambda values for each type.
1481 All that are nonzero must have the same number */
1483 for (i = 0; i < efptNR; i++)
1485 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1488 /* now, determine the number of components. All must be either zero, or equal. */
1491 for (i = 0; i < efptNR; i++)
1493 if (nfep[i] > max_n_lambda)
1495 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1496 must have the same number if its not zero.*/
1501 for (i = 0; i < efptNR; i++)
1505 ir->fepvals->separate_dvdl[i] = FALSE;
1507 else if (nfep[i] == max_n_lambda)
1509 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1510 respect to the temperature currently */
1512 ir->fepvals->separate_dvdl[i] = TRUE;
1517 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1518 nfep[i], efpt_names[i], max_n_lambda);
1521 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1522 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1524 /* the number of lambdas is the number we've read in, which is either zero
1525 or the same for all */
1526 fep->n_lambda = max_n_lambda;
1528 /* allocate space for the array of lambda values */
1529 snew(fep->all_lambda, efptNR);
1530 /* if init_lambda is defined, we need to set lambda */
1531 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1533 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1535 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1536 for (i = 0; i < efptNR; i++)
1538 snew(fep->all_lambda[i], fep->n_lambda);
1539 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1542 for (j = 0; j < fep->n_lambda; j++)
1544 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1546 sfree(count_fep_lambdas[i]);
1549 sfree(count_fep_lambdas);
1551 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1552 bookkeeping -- for now, init_lambda */
1554 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1556 for (i = 0; i < fep->n_lambda; i++)
1558 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1562 /* check to see if only a single component lambda is defined, and soft core is defined.
1563 In this case, turn on coulomb soft core */
1565 if (max_n_lambda == 0)
1571 for (i = 0; i < efptNR; i++)
1573 if ((nfep[i] != 0) && (i != efptFEP))
1579 if ((bOneLambda) && (fep->sc_alpha > 0))
1581 fep->bScCoul = TRUE;
1584 /* Fill in the others with the efptFEP if they are not explicitly
1585 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1586 they are all zero. */
1588 for (i = 0; i < efptNR; i++)
1590 if ((nfep[i] == 0) && (i != efptFEP))
1592 for (j = 0; j < fep->n_lambda; j++)
1594 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1600 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1601 if (fep->sc_r_power == 48)
1603 if (fep->sc_alpha > 0.1)
1605 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1609 /* now read in the weights */
1610 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1613 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1615 else if (nweights != fep->n_lambda)
1617 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1618 nweights, fep->n_lambda);
1620 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1622 expand->nstexpanded = fep->nstdhdl;
1623 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1625 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1627 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1628 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1629 2*tau_t just to be careful so it's not to frequent */
1634 static void do_simtemp_params(t_inputrec *ir)
1637 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1638 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1643 static void do_wall_params(t_inputrec *ir,
1644 char *wall_atomtype, char *wall_density,
1648 char *names[MAXPTR];
1651 opts->wall_atomtype[0] = nullptr;
1652 opts->wall_atomtype[1] = nullptr;
1654 ir->wall_atomtype[0] = -1;
1655 ir->wall_atomtype[1] = -1;
1656 ir->wall_density[0] = 0;
1657 ir->wall_density[1] = 0;
1661 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1662 if (nstr != ir->nwall)
1664 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1667 for (i = 0; i < ir->nwall; i++)
1669 opts->wall_atomtype[i] = gmx_strdup(names[i]);
1672 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1674 nstr = str_nelem(wall_density, MAXPTR, names);
1675 if (nstr != ir->nwall)
1677 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1679 for (i = 0; i < ir->nwall; i++)
1681 if (sscanf(names[i], "%lf", &dbl) != 1)
1683 gmx_fatal(FARGS, "Could not parse wall-density value from string '%s'", names[i]);
1687 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1689 ir->wall_density[i] = dbl;
1695 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1703 srenew(groups->grpname, groups->ngrpname+nwall);
1704 grps = &(groups->grps[egcENER]);
1705 srenew(grps->nm_ind, grps->nr+nwall);
1706 for (i = 0; i < nwall; i++)
1708 sprintf(str, "wall%d", i);
1709 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1710 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1715 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1716 t_expanded *expand, warninp_t wi)
1724 /* read expanded ensemble parameters */
1725 CCTYPE ("expanded ensemble variables");
1726 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1727 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1728 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1729 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1730 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1731 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1732 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1733 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1734 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1735 CCTYPE("Seed for Monte Carlo in lambda space");
1736 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1737 RTYPE ("mc-temperature", expand->mc_temp, -1);
1738 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1739 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1740 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1741 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1742 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1743 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1744 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1745 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1746 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1747 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1748 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1756 /*! \brief Return whether an end state with the given coupling-lambda
1757 * value describes fully-interacting VDW.
1759 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1760 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1762 static gmx_bool couple_lambda_has_vdw_on(int couple_lambda_value)
1764 return (couple_lambda_value == ecouplamVDW ||
1765 couple_lambda_value == ecouplamVDWQ);
1771 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1774 explicit MdpErrorHandler(warninp_t wi)
1775 : wi_(wi), mapping_(nullptr)
1779 void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
1781 mapping_ = &mapping;
1784 virtual bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context)
1786 ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
1787 getOptionName(context).c_str()));
1788 std::string message = gmx::formatExceptionMessageToString(*ex);
1789 warning_error(wi_, message.c_str());
1794 std::string getOptionName(const gmx::KeyValueTreePath &context)
1796 if (mapping_ != nullptr)
1798 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1799 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1802 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1807 const gmx::IKeyValueTreeBackMapping *mapping_;
1812 void get_ir(const char *mdparin, const char *mdparout,
1813 gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
1814 WriteMdpHeader writeMdpHeader, warninp_t wi)
1817 double dumdub[2][6];
1821 char warn_buf[STRLEN];
1822 t_lambda *fep = ir->fepvals;
1823 t_expanded *expand = ir->expandedvals;
1825 init_inputrec_strings();
1826 gmx::TextInputFile stream(mdparin);
1827 inp = read_inpfile(&stream, mdparin, &ninp, wi);
1829 snew(dumstr[0], STRLEN);
1830 snew(dumstr[1], STRLEN);
1832 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1835 "%s did not specify a value for the .mdp option "
1836 "\"cutoff-scheme\". Probably it was first intended for use "
1837 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1838 "introduced, but the group scheme was still the default. "
1839 "The default is now the Verlet scheme, so you will observe "
1840 "different behaviour.", mdparin);
1841 warning_note(wi, warn_buf);
1844 /* ignore the following deprecated commands */
1847 REM_TYPE("domain-decomposition");
1848 REM_TYPE("andersen-seed");
1850 REM_TYPE("dihre-fc");
1851 REM_TYPE("dihre-tau");
1852 REM_TYPE("nstdihreout");
1853 REM_TYPE("nstcheckpoint");
1854 REM_TYPE("optimize-fft");
1855 REM_TYPE("adress_type");
1856 REM_TYPE("adress_const_wf");
1857 REM_TYPE("adress_ex_width");
1858 REM_TYPE("adress_hy_width");
1859 REM_TYPE("adress_ex_forcecap");
1860 REM_TYPE("adress_interface_correction");
1861 REM_TYPE("adress_site");
1862 REM_TYPE("adress_reference_coords");
1863 REM_TYPE("adress_tf_grp_names");
1864 REM_TYPE("adress_cg_grp_names");
1865 REM_TYPE("adress_do_hybridpairs");
1866 REM_TYPE("rlistlong");
1867 REM_TYPE("nstcalclr");
1868 REM_TYPE("pull-print-com2");
1870 /* replace the following commands with the clearer new versions*/
1871 REPL_TYPE("unconstrained-start", "continuation");
1872 REPL_TYPE("foreign-lambda", "fep-lambdas");
1873 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1874 REPL_TYPE("nstxtcout", "nstxout-compressed");
1875 REPL_TYPE("xtc-grps", "compressed-x-grps");
1876 REPL_TYPE("xtc-precision", "compressed-x-precision");
1877 REPL_TYPE("pull-print-com1", "pull-print-com");
1879 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1880 CTYPE ("Preprocessor information: use cpp syntax.");
1881 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1882 STYPE ("include", opts->include, nullptr);
1883 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1884 STYPE ("define", opts->define, nullptr);
1886 CCTYPE ("RUN CONTROL PARAMETERS");
1887 EETYPE("integrator", ir->eI, ei_names);
1888 CTYPE ("Start time and timestep in ps");
1889 RTYPE ("tinit", ir->init_t, 0.0);
1890 RTYPE ("dt", ir->delta_t, 0.001);
1891 STEPTYPE ("nsteps", ir->nsteps, 0);
1892 CTYPE ("For exact run continuation or redoing part of a run");
1893 STEPTYPE ("init-step", ir->init_step, 0);
1894 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1895 ITYPE ("simulation-part", ir->simulation_part, 1);
1896 CTYPE ("mode for center of mass motion removal");
1897 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1898 CTYPE ("number of steps for center of mass motion removal");
1899 ITYPE ("nstcomm", ir->nstcomm, 100);
1900 CTYPE ("group(s) for center of mass motion removal");
1901 STYPE ("comm-grps", is->vcm, nullptr);
1903 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1904 CTYPE ("Friction coefficient (amu/ps) and random seed");
1905 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1906 STEPTYPE ("ld-seed", ir->ld_seed, -1);
1909 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1910 CTYPE ("Force tolerance and initial step-size");
1911 RTYPE ("emtol", ir->em_tol, 10.0);
1912 RTYPE ("emstep", ir->em_stepsize, 0.01);
1913 CTYPE ("Max number of iterations in relax-shells");
1914 ITYPE ("niter", ir->niter, 20);
1915 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1916 RTYPE ("fcstep", ir->fc_stepsize, 0);
1917 CTYPE ("Frequency of steepest descents steps when doing CG");
1918 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1919 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1921 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1922 RTYPE ("rtpi", ir->rtpi, 0.05);
1924 /* Output options */
1925 CCTYPE ("OUTPUT CONTROL OPTIONS");
1926 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1927 ITYPE ("nstxout", ir->nstxout, 0);
1928 ITYPE ("nstvout", ir->nstvout, 0);
1929 ITYPE ("nstfout", ir->nstfout, 0);
1930 CTYPE ("Output frequency for energies to log file and energy file");
1931 ITYPE ("nstlog", ir->nstlog, 1000);
1932 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1933 ITYPE ("nstenergy", ir->nstenergy, 1000);
1934 CTYPE ("Output frequency and precision for .xtc file");
1935 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
1936 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1937 CTYPE ("This selects the subset of atoms for the compressed");
1938 CTYPE ("trajectory file. You can select multiple groups. By");
1939 CTYPE ("default, all atoms will be written.");
1940 STYPE ("compressed-x-grps", is->x_compressed_groups, nullptr);
1941 CTYPE ("Selection of energy groups");
1942 STYPE ("energygrps", is->energy, nullptr);
1944 /* Neighbor searching */
1945 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1946 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1947 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1948 CTYPE ("nblist update frequency");
1949 ITYPE ("nstlist", ir->nstlist, 10);
1950 CTYPE ("ns algorithm (simple or grid)");
1951 EETYPE("ns-type", ir->ns_type, ens_names);
1952 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1953 EETYPE("pbc", ir->ePBC, epbc_names);
1954 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1955 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1956 CTYPE ("a value of -1 means: use rlist");
1957 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1958 CTYPE ("nblist cut-off");
1959 RTYPE ("rlist", ir->rlist, 1.0);
1960 CTYPE ("long-range cut-off for switched potentials");
1962 /* Electrostatics */
1963 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1964 CTYPE ("Method for doing electrostatics");
1965 EETYPE("coulombtype", ir->coulombtype, eel_names);
1966 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1967 CTYPE ("cut-off lengths");
1968 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1969 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1970 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1971 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1972 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1973 CTYPE ("Method for doing Van der Waals");
1974 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1975 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1976 CTYPE ("cut-off lengths");
1977 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1978 RTYPE ("rvdw", ir->rvdw, 1.0);
1979 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1980 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1981 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1982 RTYPE ("table-extension", ir->tabext, 1.0);
1983 CTYPE ("Separate tables between energy group pairs");
1984 STYPE ("energygrp-table", is->egptable, nullptr);
1985 CTYPE ("Spacing for the PME/PPPM FFT grid");
1986 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1987 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1988 ITYPE ("fourier-nx", ir->nkx, 0);
1989 ITYPE ("fourier-ny", ir->nky, 0);
1990 ITYPE ("fourier-nz", ir->nkz, 0);
1991 CTYPE ("EWALD/PME/PPPM parameters");
1992 ITYPE ("pme-order", ir->pme_order, 4);
1993 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1994 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1995 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1996 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1997 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1999 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
2000 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
2002 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
2003 CTYPE ("Algorithm for calculating Born radii");
2004 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
2005 CTYPE ("Frequency of calculating the Born radii inside rlist");
2006 ITYPE ("nstgbradii", ir->nstgbradii, 1);
2007 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
2008 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
2009 RTYPE ("rgbradii", ir->rgbradii, 1.0);
2010 CTYPE ("Dielectric coefficient of the implicit solvent");
2011 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
2012 CTYPE ("Salt concentration in M for Generalized Born models");
2013 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
2014 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
2015 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
2016 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
2017 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
2018 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
2019 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
2020 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
2021 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
2022 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
2024 /* Coupling stuff */
2025 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
2026 CTYPE ("Temperature coupling");
2027 EETYPE("tcoupl", ir->etc, etcoupl_names);
2028 ITYPE ("nsttcouple", ir->nsttcouple, -1);
2029 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
2030 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
2031 CTYPE ("Groups to couple separately");
2032 STYPE ("tc-grps", is->tcgrps, nullptr);
2033 CTYPE ("Time constant (ps) and reference temperature (K)");
2034 STYPE ("tau-t", is->tau_t, nullptr);
2035 STYPE ("ref-t", is->ref_t, nullptr);
2036 CTYPE ("pressure coupling");
2037 EETYPE("pcoupl", ir->epc, epcoupl_names);
2038 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
2039 ITYPE ("nstpcouple", ir->nstpcouple, -1);
2040 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
2041 RTYPE ("tau-p", ir->tau_p, 1.0);
2042 STYPE ("compressibility", dumstr[0], nullptr);
2043 STYPE ("ref-p", dumstr[1], nullptr);
2044 CTYPE ("Scaling of reference coordinates, No, All or COM");
2045 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
2048 CCTYPE ("OPTIONS FOR QMMM calculations");
2049 EETYPE("QMMM", ir->bQMMM, yesno_names);
2050 CTYPE ("Groups treated Quantum Mechanically");
2051 STYPE ("QMMM-grps", is->QMMM, nullptr);
2052 CTYPE ("QM method");
2053 STYPE("QMmethod", is->QMmethod, nullptr);
2054 CTYPE ("QMMM scheme");
2055 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
2056 CTYPE ("QM basisset");
2057 STYPE("QMbasis", is->QMbasis, nullptr);
2058 CTYPE ("QM charge");
2059 STYPE ("QMcharge", is->QMcharge, nullptr);
2060 CTYPE ("QM multiplicity");
2061 STYPE ("QMmult", is->QMmult, nullptr);
2062 CTYPE ("Surface Hopping");
2063 STYPE ("SH", is->bSH, nullptr);
2064 CTYPE ("CAS space options");
2065 STYPE ("CASorbitals", is->CASorbitals, nullptr);
2066 STYPE ("CASelectrons", is->CASelectrons, nullptr);
2067 STYPE ("SAon", is->SAon, nullptr);
2068 STYPE ("SAoff", is->SAoff, nullptr);
2069 STYPE ("SAsteps", is->SAsteps, nullptr);
2070 CTYPE ("Scale factor for MM charges");
2071 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
2072 CTYPE ("Optimization of QM subsystem");
2073 STYPE ("bOPT", is->bOPT, nullptr);
2074 STYPE ("bTS", is->bTS, nullptr);
2076 /* Simulated annealing */
2077 CCTYPE("SIMULATED ANNEALING");
2078 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
2079 STYPE ("annealing", is->anneal, nullptr);
2080 CTYPE ("Number of time points to use for specifying annealing in each group");
2081 STYPE ("annealing-npoints", is->anneal_npoints, nullptr);
2082 CTYPE ("List of times at the annealing points for each group");
2083 STYPE ("annealing-time", is->anneal_time, nullptr);
2084 CTYPE ("Temp. at each annealing point, for each group.");
2085 STYPE ("annealing-temp", is->anneal_temp, nullptr);
2088 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
2089 EETYPE("gen-vel", opts->bGenVel, yesno_names);
2090 RTYPE ("gen-temp", opts->tempi, 300.0);
2091 ITYPE ("gen-seed", opts->seed, -1);
2094 CCTYPE ("OPTIONS FOR BONDS");
2095 EETYPE("constraints", opts->nshake, constraints);
2096 CTYPE ("Type of constraint algorithm");
2097 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
2098 CTYPE ("Do not constrain the start configuration");
2099 EETYPE("continuation", ir->bContinuation, yesno_names);
2100 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2101 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2102 CTYPE ("Relative tolerance of shake");
2103 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2104 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2105 ITYPE ("lincs-order", ir->nProjOrder, 4);
2106 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2107 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2108 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2109 ITYPE ("lincs-iter", ir->nLincsIter, 1);
2110 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2111 CTYPE ("rotates over more degrees than");
2112 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2113 CTYPE ("Convert harmonic bonds to morse potentials");
2114 EETYPE("morse", opts->bMorse, yesno_names);
2116 /* Energy group exclusions */
2117 CCTYPE ("ENERGY GROUP EXCLUSIONS");
2118 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2119 STYPE ("energygrp-excl", is->egpexcl, nullptr);
2123 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2124 ITYPE ("nwall", ir->nwall, 0);
2125 EETYPE("wall-type", ir->wall_type, ewt_names);
2126 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2127 STYPE ("wall-atomtype", is->wall_atomtype, nullptr);
2128 STYPE ("wall-density", is->wall_density, nullptr);
2129 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2132 CCTYPE("COM PULLING");
2133 EETYPE("pull", ir->bPull, yesno_names);
2137 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, wi);
2140 /* Enforced rotation */
2141 CCTYPE("ENFORCED ROTATION");
2142 CTYPE("Enforced rotation: No or Yes");
2143 EETYPE("rotation", ir->bRot, yesno_names);
2147 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2150 /* Interactive MD */
2152 CCTYPE("Group to display and/or manipulate in interactive MD session");
2153 STYPE ("IMD-group", is->imd_grp, nullptr);
2154 if (is->imd_grp[0] != '\0')
2161 CCTYPE("NMR refinement stuff");
2162 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2163 EETYPE("disre", ir->eDisre, edisre_names);
2164 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2165 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2166 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2167 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2168 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2169 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2170 CTYPE ("Output frequency for pair distances to energy file");
2171 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2172 CTYPE ("Orientation restraints: No or Yes");
2173 EETYPE("orire", opts->bOrire, yesno_names);
2174 CTYPE ("Orientation restraints force constant and tau for time averaging");
2175 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2176 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2177 STYPE ("orire-fitgrp", is->orirefitgrp, nullptr);
2178 CTYPE ("Output frequency for trace(SD) and S to energy file");
2179 ITYPE ("nstorireout", ir->nstorireout, 100);
2181 /* free energy variables */
2182 CCTYPE ("Free energy variables");
2183 EETYPE("free-energy", ir->efep, efep_names);
2184 STYPE ("couple-moltype", is->couple_moltype, nullptr);
2185 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2186 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2187 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2189 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2191 it was not entered */
2192 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2193 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2194 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2195 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], nullptr);
2196 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], nullptr);
2197 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
2198 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
2199 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
2200 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
2201 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
2202 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2203 STYPE ("init-lambda-weights", is->lambda_weights, nullptr);
2204 EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
2205 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2206 ITYPE ("sc-power", fep->sc_power, 1);
2207 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2208 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2209 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2210 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2211 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2212 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2213 separate_dhdl_file_names);
2214 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2215 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2216 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2218 /* Non-equilibrium MD stuff */
2219 CCTYPE("Non-equilibrium MD stuff");
2220 STYPE ("acc-grps", is->accgrps, nullptr);
2221 STYPE ("accelerate", is->acc, nullptr);
2222 STYPE ("freezegrps", is->freeze, nullptr);
2223 STYPE ("freezedim", is->frdim, nullptr);
2224 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2225 STYPE ("deform", is->deform, nullptr);
2227 /* simulated tempering variables */
2228 CCTYPE("simulated tempering variables");
2229 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2230 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2231 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2232 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2234 /* expanded ensemble variables */
2235 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2237 read_expandedparams(&ninp, &inp, expand, wi);
2240 /* Electric fields */
2242 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(ninp, inp);
2243 gmx::KeyValueTreeTransformer transform;
2244 transform.rules()->addRule()
2245 .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2246 mdModules->initMdpTransform(transform.rules());
2247 for (const auto &path : transform.mappedPaths())
2249 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2250 mark_einp_set(ninp, inp, path[0].c_str());
2252 MdpErrorHandler errorHandler(wi);
2254 = transform.transform(convertedValues, &errorHandler);
2255 ir->params = new gmx::KeyValueTreeObject(result.object());
2256 mdModules->adjustInputrecBasedOnModules(ir);
2257 errorHandler.setBackMapping(result.backMapping());
2258 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2261 /* Ion/water position swapping ("computational electrophysiology") */
2262 CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2263 CTYPE("Swap positions along direction: no, X, Y, Z");
2264 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2265 if (ir->eSwapCoords != eswapNO)
2272 CTYPE("Swap attempt frequency");
2273 ITYPE("swap-frequency", ir->swap->nstswap, 1);
2274 CTYPE("Number of ion types to be controlled");
2275 ITYPE("iontypes", nIonTypes, 1);
2278 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2280 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2281 snew(ir->swap->grp, ir->swap->ngrp);
2282 for (i = 0; i < ir->swap->ngrp; i++)
2284 snew(ir->swap->grp[i].molname, STRLEN);
2286 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2287 STYPE("split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2288 STYPE("split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2289 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2290 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2291 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2293 CTYPE("Name of solvent molecules");
2294 STYPE("solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2296 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2297 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2298 CTYPE("however, if correctly defined, the permeation events are recorded per channel");
2299 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2300 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2301 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2302 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2303 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2304 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2306 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2307 ITYPE("coupl-steps", ir->swap->nAverage, 10);
2309 CTYPE("Names of the ion types that can be exchanged with solvent molecules,");
2310 CTYPE("and the requested number of ions of this type in compartments A and B");
2311 CTYPE("-1 means fix the numbers as found in step 0");
2312 for (i = 0; i < nIonTypes; i++)
2314 int ig = eSwapFixedGrpNR + i;
2316 sprintf(buf, "iontype%d-name", i);
2317 STYPE(buf, ir->swap->grp[ig].molname, nullptr);
2318 sprintf(buf, "iontype%d-in-A", i);
2319 ITYPE(buf, ir->swap->grp[ig].nmolReq[0], -1);
2320 sprintf(buf, "iontype%d-in-B", i);
2321 ITYPE(buf, ir->swap->grp[ig].nmolReq[1], -1);
2324 CTYPE("By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2325 CTYPE("at maximum distance (= bulk concentration) to the split group layers. However,");
2326 CTYPE("an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2327 CTYPE("towards one of the compartment-partitioning layers (at +/- 1.0).");
2328 RTYPE("bulk-offsetA", ir->swap->bulkOffset[0], 0.0);
2329 RTYPE("bulk-offsetB", ir->swap->bulkOffset[1], 0.0);
2330 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2331 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2333 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2336 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2337 RTYPE("threshold", ir->swap->threshold, 1.0);
2340 /* AdResS is no longer supported, but we need mdrun to be able to refuse to run old AdResS .tpr files */
2341 EETYPE("adress", ir->bAdress, yesno_names);
2343 /* User defined thingies */
2344 CCTYPE ("User defined thingies");
2345 STYPE ("user1-grps", is->user1, nullptr);
2346 STYPE ("user2-grps", is->user2, nullptr);
2347 ITYPE ("userint1", ir->userint1, 0);
2348 ITYPE ("userint2", ir->userint2, 0);
2349 ITYPE ("userint3", ir->userint3, 0);
2350 ITYPE ("userint4", ir->userint4, 0);
2351 RTYPE ("userreal1", ir->userreal1, 0);
2352 RTYPE ("userreal2", ir->userreal2, 0);
2353 RTYPE ("userreal3", ir->userreal3, 0);
2354 RTYPE ("userreal4", ir->userreal4, 0);
2358 gmx::TextOutputFile stream(mdparout);
2359 write_inpfile(&stream, mdparout, ninp, inp, FALSE, writeMdpHeader, wi);
2363 for (i = 0; (i < ninp); i++)
2366 sfree(inp[i].value);
2370 /* Process options if necessary */
2371 for (m = 0; m < 2; m++)
2373 for (i = 0; i < 2*DIM; i++)
2382 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2384 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
2386 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2388 case epctSEMIISOTROPIC:
2389 case epctSURFACETENSION:
2390 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2392 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
2394 dumdub[m][YY] = dumdub[m][XX];
2396 case epctANISOTROPIC:
2397 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2398 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2399 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2401 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
2405 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2406 epcoupltype_names[ir->epct]);
2410 clear_mat(ir->ref_p);
2411 clear_mat(ir->compress);
2412 for (i = 0; i < DIM; i++)
2414 ir->ref_p[i][i] = dumdub[1][i];
2415 ir->compress[i][i] = dumdub[0][i];
2417 if (ir->epct == epctANISOTROPIC)
2419 ir->ref_p[XX][YY] = dumdub[1][3];
2420 ir->ref_p[XX][ZZ] = dumdub[1][4];
2421 ir->ref_p[YY][ZZ] = dumdub[1][5];
2422 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2424 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2426 ir->compress[XX][YY] = dumdub[0][3];
2427 ir->compress[XX][ZZ] = dumdub[0][4];
2428 ir->compress[YY][ZZ] = dumdub[0][5];
2429 for (i = 0; i < DIM; i++)
2431 for (m = 0; m < i; m++)
2433 ir->ref_p[i][m] = ir->ref_p[m][i];
2434 ir->compress[i][m] = ir->compress[m][i];
2439 if (ir->comm_mode == ecmNO)
2444 opts->couple_moltype = nullptr;
2445 if (strlen(is->couple_moltype) > 0)
2447 if (ir->efep != efepNO)
2449 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2450 if (opts->couple_lam0 == opts->couple_lam1)
2452 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2454 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2455 opts->couple_lam1 == ecouplamNONE))
2457 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2462 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2465 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2466 if (ir->efep != efepNO)
2468 if (fep->delta_lambda > 0)
2470 ir->efep = efepSLOWGROWTH;
2474 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2476 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2477 warning_note(wi, "Old option for dhdl-print-energy given: "
2478 "changing \"yes\" to \"total\"\n");
2481 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2483 /* always print out the energy to dhdl if we are doing
2484 expanded ensemble, since we need the total energy for
2485 analysis if the temperature is changing. In some
2486 conditions one may only want the potential energy, so
2487 we will allow that if the appropriate mdp setting has
2488 been enabled. Otherwise, total it is:
2490 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2493 if ((ir->efep != efepNO) || ir->bSimTemp)
2495 ir->bExpanded = FALSE;
2496 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2498 ir->bExpanded = TRUE;
2500 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2501 if (ir->bSimTemp) /* done after fep params */
2503 do_simtemp_params(ir);
2506 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2507 * setup and not on the old way of specifying the free-energy setup,
2508 * we should check for using soft-core when not needed, since that
2509 * can complicate the sampling significantly.
2510 * Note that we only check for the automated coupling setup.
2511 * If the (advanced) user does FEP through manual topology changes,
2512 * this check will not be triggered.
2514 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2515 ir->fepvals->sc_alpha != 0 &&
2516 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2517 couple_lambda_has_vdw_on(opts->couple_lam1)))
2519 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.");
2524 ir->fepvals->n_lambda = 0;
2527 /* WALL PARAMETERS */
2529 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2531 /* ORIENTATION RESTRAINT PARAMETERS */
2533 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, nullptr) != 1)
2535 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2538 /* DEFORMATION PARAMETERS */
2540 clear_mat(ir->deform);
2541 for (i = 0; i < 6; i++)
2546 double gmx_unused canary;
2547 int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
2548 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2549 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
2551 if (strlen(is->deform) > 0 && ndeform != 6)
2553 warning_error(wi, gmx::formatString("Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform).c_str());
2555 for (i = 0; i < 3; i++)
2557 ir->deform[i][i] = dumdub[0][i];
2559 ir->deform[YY][XX] = dumdub[0][3];
2560 ir->deform[ZZ][XX] = dumdub[0][4];
2561 ir->deform[ZZ][YY] = dumdub[0][5];
2562 if (ir->epc != epcNO)
2564 for (i = 0; i < 3; i++)
2566 for (j = 0; j <= i; j++)
2568 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2570 warning_error(wi, "A box element has deform set and compressibility > 0");
2574 for (i = 0; i < 3; i++)
2576 for (j = 0; j < i; j++)
2578 if (ir->deform[i][j] != 0)
2580 for (m = j; m < DIM; m++)
2582 if (ir->compress[m][j] != 0)
2584 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.");
2585 warning(wi, warn_buf);
2593 /* Ion/water position swapping checks */
2594 if (ir->eSwapCoords != eswapNO)
2596 if (ir->swap->nstswap < 1)
2598 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2600 if (ir->swap->nAverage < 1)
2602 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2604 if (ir->swap->threshold < 1.0)
2606 warning_error(wi, "Ion count threshold must be at least 1.\n");
2614 static int search_QMstring(const char *s, int ng, const char *gn[])
2616 /* same as normal search_string, but this one searches QM strings */
2619 for (i = 0; (i < ng); i++)
2621 if (gmx_strcasecmp(s, gn[i]) == 0)
2627 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2631 } /* search_QMstring */
2633 /* We would like gn to be const as well, but C doesn't allow this */
2634 /* TODO this is utility functionality (search for the index of a
2635 string in a collection), so should be refactored and located more
2637 int search_string(const char *s, int ng, char *gn[])
2641 for (i = 0; (i < ng); i++)
2643 if (gmx_strcasecmp(s, gn[i]) == 0)
2650 "Group %s referenced in the .mdp file was not found in the index file.\n"
2651 "Group names must match either [moleculetype] names or custom index group\n"
2652 "names, in which case you must supply an index file to the '-n' option\n"
2659 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2660 t_blocka *block, char *gnames[],
2661 int gtype, int restnm,
2662 int grptp, gmx_bool bVerbose,
2665 unsigned short *cbuf;
2666 t_grps *grps = &(groups->grps[gtype]);
2667 int i, j, gid, aj, ognr, ntot = 0;
2670 char warn_buf[STRLEN];
2674 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2677 title = gtypes[gtype];
2680 /* Mark all id's as not set */
2681 for (i = 0; (i < natoms); i++)
2686 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2687 for (i = 0; (i < ng); i++)
2689 /* Lookup the group name in the block structure */
2690 gid = search_string(ptrs[i], block->nr, gnames);
2691 if ((grptp != egrptpONE) || (i == 0))
2693 grps->nm_ind[grps->nr++] = gid;
2697 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2700 /* Now go over the atoms in the group */
2701 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2706 /* Range checking */
2707 if ((aj < 0) || (aj >= natoms))
2709 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2711 /* Lookup up the old group number */
2715 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2716 aj+1, title, ognr+1, i+1);
2720 /* Store the group number in buffer */
2721 if (grptp == egrptpONE)
2734 /* Now check whether we have done all atoms */
2738 if (grptp == egrptpALL)
2740 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2741 natoms-ntot, title);
2743 else if (grptp == egrptpPART)
2745 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2746 natoms-ntot, title);
2747 warning_note(wi, warn_buf);
2749 /* Assign all atoms currently unassigned to a rest group */
2750 for (j = 0; (j < natoms); j++)
2752 if (cbuf[j] == NOGID)
2758 if (grptp != egrptpPART)
2763 "Making dummy/rest group for %s containing %d elements\n",
2764 title, natoms-ntot);
2766 /* Add group name "rest" */
2767 grps->nm_ind[grps->nr] = restnm;
2769 /* Assign the rest name to all atoms not currently assigned to a group */
2770 for (j = 0; (j < natoms); j++)
2772 if (cbuf[j] == NOGID)
2781 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2783 /* All atoms are part of one (or no) group, no index required */
2784 groups->ngrpnr[gtype] = 0;
2785 groups->grpnr[gtype] = nullptr;
2789 groups->ngrpnr[gtype] = natoms;
2790 snew(groups->grpnr[gtype], natoms);
2791 for (j = 0; (j < natoms); j++)
2793 groups->grpnr[gtype][j] = cbuf[j];
2799 return (bRest && grptp == egrptpPART);
2802 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2805 gmx_groups_t *groups;
2806 pull_params_t *pull;
2807 int natoms, ai, aj, i, j, d, g, imin, jmin;
2809 int *nrdf2, *na_vcm, na_tot;
2810 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2812 gmx_mtop_atomloop_all_t aloop;
2813 int mb, mol, ftype, as;
2814 gmx_molblock_t *molb;
2815 gmx_moltype_t *molt;
2818 * First calc 3xnr-atoms for each group
2819 * then subtract half a degree of freedom for each constraint
2821 * Only atoms and nuclei contribute to the degrees of freedom...
2826 groups = &mtop->groups;
2827 natoms = mtop->natoms;
2829 /* Allocate one more for a possible rest group */
2830 /* We need to sum degrees of freedom into doubles,
2831 * since floats give too low nrdf's above 3 million atoms.
2833 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2834 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2835 snew(dof_vcm, groups->grps[egcVCM].nr+1);
2836 snew(na_vcm, groups->grps[egcVCM].nr+1);
2837 snew(nrdf_vcm_sub, groups->grps[egcVCM].nr+1);
2839 for (i = 0; i < groups->grps[egcTC].nr; i++)
2843 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2846 clear_ivec(dof_vcm[i]);
2848 nrdf_vcm_sub[i] = 0;
2851 snew(nrdf2, natoms);
2852 aloop = gmx_mtop_atomloop_all_init(mtop);
2854 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2857 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2859 g = ggrpnr(groups, egcFREEZE, i);
2860 for (d = 0; d < DIM; d++)
2862 if (opts->nFreeze[g][d] == 0)
2864 /* Add one DOF for particle i (counted as 2*1) */
2866 /* VCM group i has dim d as a DOF */
2867 dof_vcm[ggrpnr(groups, egcVCM, i)][d] = 1;
2870 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2871 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2876 for (mb = 0; mb < mtop->nmolblock; mb++)
2878 molb = &mtop->molblock[mb];
2879 molt = &mtop->moltype[molb->type];
2880 atom = molt->atoms.atom;
2881 for (mol = 0; mol < molb->nmol; mol++)
2883 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2885 ia = molt->ilist[ftype].iatoms;
2886 for (i = 0; i < molt->ilist[ftype].nr; )
2888 /* Subtract degrees of freedom for the constraints,
2889 * if the particles still have degrees of freedom left.
2890 * If one of the particles is a vsite or a shell, then all
2891 * constraint motion will go there, but since they do not
2892 * contribute to the constraints the degrees of freedom do not
2897 if (((atom[ia[1]].ptype == eptNucleus) ||
2898 (atom[ia[1]].ptype == eptAtom)) &&
2899 ((atom[ia[2]].ptype == eptNucleus) ||
2900 (atom[ia[2]].ptype == eptAtom)))
2918 imin = std::min(imin, nrdf2[ai]);
2919 jmin = std::min(jmin, nrdf2[aj]);
2922 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2923 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2924 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2925 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2927 ia += interaction_function[ftype].nratoms+1;
2928 i += interaction_function[ftype].nratoms+1;
2931 ia = molt->ilist[F_SETTLE].iatoms;
2932 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2934 /* Subtract 1 dof from every atom in the SETTLE */
2935 for (j = 0; j < 3; j++)
2938 imin = std::min(2, nrdf2[ai]);
2940 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2941 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2946 as += molt->atoms.nr;
2952 /* Correct nrdf for the COM constraints.
2953 * We correct using the TC and VCM group of the first atom
2954 * in the reference and pull group. If atoms in one pull group
2955 * belong to different TC or VCM groups it is anyhow difficult
2956 * to determine the optimal nrdf assignment.
2960 for (i = 0; i < pull->ncoord; i++)
2962 if (pull->coord[i].eType != epullCONSTRAINT)
2969 for (j = 0; j < 2; j++)
2971 const t_pull_group *pgrp;
2973 pgrp = &pull->group[pull->coord[i].group[j]];
2977 /* Subtract 1/2 dof from each group */
2979 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2980 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2981 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2983 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->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]);
2988 /* We need to subtract the whole DOF from group j=1 */
2995 if (ir->nstcomm != 0)
2999 /* We remove COM motion up to dim ndof_com() */
3000 ndim_rm_vcm = ndof_com(ir);
3002 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3003 * the number of degrees of freedom in each vcm group when COM
3004 * translation is removed and 6 when rotation is removed as well.
3006 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3008 switch (ir->comm_mode)
3011 nrdf_vcm_sub[j] = 0;
3012 for (d = 0; d < ndim_rm_vcm; d++)
3021 nrdf_vcm_sub[j] = 6;
3024 gmx_incons("Checking comm_mode");
3028 for (i = 0; i < groups->grps[egcTC].nr; i++)
3030 /* Count the number of atoms of TC group i for every VCM group */
3031 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3036 for (ai = 0; ai < natoms; ai++)
3038 if (ggrpnr(groups, egcTC, ai) == i)
3040 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
3044 /* Correct for VCM removal according to the fraction of each VCM
3045 * group present in this TC group.
3047 nrdf_uc = nrdf_tc[i];
3050 fprintf(debug, "T-group[%d] nrdf_uc = %g\n", i, nrdf_uc);
3053 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3055 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3057 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
3058 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
3062 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
3063 j, nrdf_vcm[j], nrdf_tc[i]);
3068 for (i = 0; (i < groups->grps[egcTC].nr); i++)
3070 opts->nrdf[i] = nrdf_tc[i];
3071 if (opts->nrdf[i] < 0)
3076 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3077 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
3085 sfree(nrdf_vcm_sub);
3088 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
3089 const char *option, const char *val, int flag)
3091 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3092 * But since this is much larger than STRLEN, such a line can not be parsed.
3093 * The real maximum is the number of names that fit in a string: STRLEN/2.
3095 #define EGP_MAX (STRLEN/2)
3096 int nelem, i, j, k, nr;
3097 char *names[EGP_MAX];
3101 gnames = groups->grpname;
3103 nelem = str_nelem(val, EGP_MAX, names);
3106 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3108 nr = groups->grps[egcENER].nr;
3110 for (i = 0; i < nelem/2; i++)
3114 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
3120 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3121 names[2*i], option);
3125 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3131 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3132 names[2*i+1], option);
3134 if ((j < nr) && (k < nr))
3136 ir->opts.egp_flags[nr*j+k] |= flag;
3137 ir->opts.egp_flags[nr*k+j] |= flag;
3146 static void make_swap_groups(
3151 int ig = -1, i = 0, gind;
3155 /* Just a quick check here, more thorough checks are in mdrun */
3156 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3158 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3161 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3162 for (ig = 0; ig < swap->ngrp; ig++)
3164 swapg = &swap->grp[ig];
3165 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3166 swapg->nat = grps->index[gind+1] - grps->index[gind];
3170 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3171 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3172 swap->grp[ig].molname, swapg->nat);
3173 snew(swapg->ind, swapg->nat);
3174 for (i = 0; i < swapg->nat; i++)
3176 swapg->ind[i] = grps->a[grps->index[gind]+i];
3181 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3187 void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3192 ig = search_string(IMDgname, grps->nr, gnames);
3193 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3195 if (IMDgroup->nat > 0)
3197 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3198 IMDgname, IMDgroup->nat);
3199 snew(IMDgroup->ind, IMDgroup->nat);
3200 for (i = 0; i < IMDgroup->nat; i++)
3202 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3208 void do_index(const char* mdparin, const char *ndx,
3215 gmx_groups_t *groups;
3219 char warnbuf[STRLEN], **gnames;
3220 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3223 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3224 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3225 int i, j, k, restnm;
3226 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3227 int nQMmethod, nQMbasis, nQMg;
3228 char warn_buf[STRLEN];
3233 fprintf(stderr, "processing index file...\n");
3238 snew(grps->index, 1);
3240 atoms_all = gmx_mtop_global_atoms(mtop);
3241 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3242 done_atom(&atoms_all);
3246 grps = init_index(ndx, &gnames);
3249 groups = &mtop->groups;
3250 natoms = mtop->natoms;
3251 symtab = &mtop->symtab;
3253 snew(groups->grpname, grps->nr+1);
3255 for (i = 0; (i < grps->nr); i++)
3257 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3259 groups->grpname[i] = put_symtab(symtab, "rest");
3261 srenew(gnames, grps->nr+1);
3262 gnames[restnm] = *(groups->grpname[i]);
3263 groups->ngrpname = grps->nr+1;
3265 set_warning_line(wi, mdparin, -1);
3267 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3268 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3269 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3270 if ((ntau_t != ntcg) || (nref_t != ntcg))
3272 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3273 "%d tau-t values", ntcg, nref_t, ntau_t);
3276 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3277 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3278 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3279 nr = groups->grps[egcTC].nr;
3281 snew(ir->opts.nrdf, nr);
3282 snew(ir->opts.tau_t, nr);
3283 snew(ir->opts.ref_t, nr);
3284 if (ir->eI == eiBD && ir->bd_fric == 0)
3286 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3293 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3297 for (i = 0; (i < nr); i++)
3299 ir->opts.tau_t[i] = strtod(ptr1[i], &endptr);
3302 warning_error(wi, "Invalid value for mdp option tau-t. tau-t should only consist of real numbers separated by spaces.");
3304 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3306 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3307 warning_error(wi, warn_buf);
3310 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3312 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.");
3315 if (ir->opts.tau_t[i] >= 0)
3317 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3320 if (ir->etc != etcNO && ir->nsttcouple == -1)
3322 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3327 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3329 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");
3331 if (ir->epc == epcMTTK)
3333 if (ir->etc != etcNOSEHOOVER)
3335 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3339 if (ir->nstpcouple != ir->nsttcouple)
3341 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3342 ir->nstpcouple = ir->nsttcouple = mincouple;
3343 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);
3344 warning_note(wi, warn_buf);
3349 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3350 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3352 if (ETC_ANDERSEN(ir->etc))
3354 if (ir->nsttcouple != 1)
3357 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");
3358 warning_note(wi, warn_buf);
3361 nstcmin = tcouple_min_integration_steps(ir->etc);
3364 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3366 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3367 ETCOUPLTYPE(ir->etc),
3369 ir->nsttcouple*ir->delta_t);
3370 warning(wi, warn_buf);
3373 for (i = 0; (i < nr); i++)
3375 ir->opts.ref_t[i] = strtod(ptr2[i], &endptr);
3378 warning_error(wi, "Invalid value for mdp option ref-t. ref-t should only consist of real numbers separated by spaces.");
3380 if (ir->opts.ref_t[i] < 0)
3382 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3385 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3386 if we are in this conditional) if mc_temp is negative */
3387 if (ir->expandedvals->mc_temp < 0)
3389 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3393 /* Simulated annealing for each group. There are nr groups */
3394 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3395 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3399 if (nSA > 0 && nSA != nr)
3401 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3405 snew(ir->opts.annealing, nr);
3406 snew(ir->opts.anneal_npoints, nr);
3407 snew(ir->opts.anneal_time, nr);
3408 snew(ir->opts.anneal_temp, nr);
3409 for (i = 0; i < nr; i++)
3411 ir->opts.annealing[i] = eannNO;
3412 ir->opts.anneal_npoints[i] = 0;
3413 ir->opts.anneal_time[i] = nullptr;
3414 ir->opts.anneal_temp[i] = nullptr;
3419 for (i = 0; i < nr; i++)
3421 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3423 ir->opts.annealing[i] = eannNO;
3425 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3427 ir->opts.annealing[i] = eannSINGLE;
3430 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3432 ir->opts.annealing[i] = eannPERIODIC;
3438 /* Read the other fields too */
3439 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3440 if (nSA_points != nSA)
3442 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3444 for (k = 0, i = 0; i < nr; i++)
3446 ir->opts.anneal_npoints[i] = strtol(ptr1[i], &endptr, 10);
3449 warning_error(wi, "Invalid value for mdp option annealing-npoints. annealing should only consist of integers separated by spaces.");
3451 if (ir->opts.anneal_npoints[i] == 1)
3453 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3455 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3456 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3457 k += ir->opts.anneal_npoints[i];
3460 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3463 gmx_fatal(FARGS, "Found %d annealing-time values, wanted %d\n", nSA_time, k);
3465 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3468 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3471 for (i = 0, k = 0; i < nr; i++)
3474 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3476 ir->opts.anneal_time[i][j] = strtod(ptr1[k], &endptr);
3479 warning_error(wi, "Invalid value for mdp option anneal-time. anneal-time should only consist of real numbers separated by spaces.");
3481 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], &endptr);
3484 warning_error(wi, "Invalid value for anneal-temp. anneal-temp should only consist of real numbers separated by spaces.");
3488 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3490 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3496 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3498 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3499 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3502 if (ir->opts.anneal_temp[i][j] < 0)
3504 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3509 /* Print out some summary information, to make sure we got it right */
3510 for (i = 0, k = 0; i < nr; i++)
3512 if (ir->opts.annealing[i] != eannNO)
3514 j = groups->grps[egcTC].nm_ind[i];
3515 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3516 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3517 ir->opts.anneal_npoints[i]);
3518 fprintf(stderr, "Time (ps) Temperature (K)\n");
3519 /* All terms except the last one */
3520 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3522 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3525 /* Finally the last one */
3526 j = ir->opts.anneal_npoints[i]-1;
3527 if (ir->opts.annealing[i] == eannSINGLE)
3529 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3533 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3534 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3536 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3547 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3549 make_pull_coords(ir->pull);
3554 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3557 if (ir->eSwapCoords != eswapNO)
3559 make_swap_groups(ir->swap, grps, gnames);
3562 /* Make indices for IMD session */
3565 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3568 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3569 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3570 if (nacg*DIM != nacc)
3572 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3575 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3576 restnm, egrptpALL_GENREST, bVerbose, wi);
3577 nr = groups->grps[egcACC].nr;
3578 snew(ir->opts.acc, nr);
3579 ir->opts.ngacc = nr;
3581 for (i = k = 0; (i < nacg); i++)
3583 for (j = 0; (j < DIM); j++, k++)
3585 ir->opts.acc[i][j] = strtod(ptr1[k], &endptr);
3588 warning_error(wi, "Invalid value for mdp option accelerate. accelerate should only consist of real numbers separated by spaces.");
3592 for (; (i < nr); i++)
3594 for (j = 0; (j < DIM); j++)
3596 ir->opts.acc[i][j] = 0;
3600 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3601 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3602 if (nfrdim != DIM*nfreeze)
3604 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3607 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3608 restnm, egrptpALL_GENREST, bVerbose, wi);
3609 nr = groups->grps[egcFREEZE].nr;
3610 ir->opts.ngfrz = nr;
3611 snew(ir->opts.nFreeze, nr);
3612 for (i = k = 0; (i < nfreeze); i++)
3614 for (j = 0; (j < DIM); j++, k++)
3616 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3617 if (!ir->opts.nFreeze[i][j])
3619 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3621 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3622 "(not %s)", ptr1[k]);
3623 warning(wi, warn_buf);
3628 for (; (i < nr); i++)
3630 for (j = 0; (j < DIM); j++)
3632 ir->opts.nFreeze[i][j] = 0;
3636 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3637 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3638 restnm, egrptpALL_GENREST, bVerbose, wi);
3639 add_wall_energrps(groups, ir->nwall, symtab);
3640 ir->opts.ngener = groups->grps[egcENER].nr;
3641 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3643 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3644 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3647 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3648 "This may lead to artifacts.\n"
3649 "In most cases one should use one group for the whole system.");
3652 /* Now we have filled the freeze struct, so we can calculate NRDF */
3653 calc_nrdf(mtop, ir, gnames);
3655 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3656 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3657 restnm, egrptpALL_GENREST, bVerbose, wi);
3658 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3659 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3660 restnm, egrptpALL_GENREST, bVerbose, wi);
3661 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3662 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3663 restnm, egrptpONE, bVerbose, wi);
3664 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3665 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3666 restnm, egrptpALL_GENREST, bVerbose, wi);
3668 /* QMMM input processing */
3669 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3670 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3671 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3672 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3674 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3675 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3677 /* group rest, if any, is always MM! */
3678 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3679 restnm, egrptpALL_GENREST, bVerbose, wi);
3680 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3681 ir->opts.ngQM = nQMg;
3682 snew(ir->opts.QMmethod, nr);
3683 snew(ir->opts.QMbasis, nr);
3684 for (i = 0; i < nr; i++)
3686 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3687 * converted to the corresponding enum in names.c
3689 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3691 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3695 str_nelem(is->QMmult, MAXPTR, ptr1);
3696 str_nelem(is->QMcharge, MAXPTR, ptr2);
3697 str_nelem(is->bSH, MAXPTR, ptr3);
3698 snew(ir->opts.QMmult, nr);
3699 snew(ir->opts.QMcharge, nr);
3700 snew(ir->opts.bSH, nr);
3702 for (i = 0; i < nr; i++)
3704 ir->opts.QMmult[i] = strtol(ptr1[i], &endptr, 10);
3707 warning_error(wi, "Invalid value for mdp option QMmult. QMmult should only consist of integers separated by spaces.");
3709 ir->opts.QMcharge[i] = strtol(ptr2[i], &endptr, 10);
3712 warning_error(wi, "Invalid value for mdp option QMcharge. QMcharge should only consist of integers separated by spaces.");
3714 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3717 str_nelem(is->CASelectrons, MAXPTR, ptr1);
3718 str_nelem(is->CASorbitals, MAXPTR, ptr2);
3719 snew(ir->opts.CASelectrons, nr);
3720 snew(ir->opts.CASorbitals, nr);
3721 for (i = 0; i < nr; i++)
3723 ir->opts.CASelectrons[i] = strtol(ptr1[i], &endptr, 10);
3726 warning_error(wi, "Invalid value for mdp option CASelectrons. CASelectrons should only consist of integers separated by spaces.");
3728 ir->opts.CASorbitals[i] = strtol(ptr2[i], &endptr, 10);
3731 warning_error(wi, "Invalid value for mdp option CASorbitals. CASorbitals should only consist of integers separated by spaces.");
3734 /* special optimization options */
3736 str_nelem(is->bOPT, MAXPTR, ptr1);
3737 str_nelem(is->bTS, MAXPTR, ptr2);
3738 snew(ir->opts.bOPT, nr);
3739 snew(ir->opts.bTS, nr);
3740 for (i = 0; i < nr; i++)
3742 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3743 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3745 str_nelem(is->SAon, MAXPTR, ptr1);
3746 str_nelem(is->SAoff, MAXPTR, ptr2);
3747 str_nelem(is->SAsteps, MAXPTR, ptr3);
3748 snew(ir->opts.SAon, nr);
3749 snew(ir->opts.SAoff, nr);
3750 snew(ir->opts.SAsteps, nr);
3752 for (i = 0; i < nr; i++)
3754 ir->opts.SAon[i] = strtod(ptr1[i], &endptr);
3757 warning_error(wi, "Invalid value for mdp option SAon. SAon should only consist of real numbers separated by spaces.");
3759 ir->opts.SAoff[i] = strtod(ptr2[i], &endptr);
3762 warning_error(wi, "Invalid value for mdp option SAoff. SAoff should only consist of real numbers separated by spaces.");
3764 ir->opts.SAsteps[i] = strtol(ptr3[i], &endptr, 10);
3767 warning_error(wi, "Invalid value for mdp option SAsteps. SAsteps should only consist of integers separated by spaces.");
3770 /* end of QMMM input */
3774 for (i = 0; (i < egcNR); i++)
3776 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3777 for (j = 0; (j < groups->grps[i].nr); j++)
3779 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3781 fprintf(stderr, "\n");
3785 nr = groups->grps[egcENER].nr;
3786 snew(ir->opts.egp_flags, nr*nr);
3788 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3789 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3791 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3793 if (bExcl && EEL_FULL(ir->coulombtype))
3795 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3798 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3799 if (bTable && !(ir->vdwtype == evdwUSER) &&
3800 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3801 !(ir->coulombtype == eelPMEUSERSWITCH))
3803 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3806 for (i = 0; (i < grps->nr); i++)
3818 static void check_disre(gmx_mtop_t *mtop)
3820 gmx_ffparams_t *ffparams;
3821 t_functype *functype;
3823 int i, ndouble, ftype;
3824 int label, old_label;
3826 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3828 ffparams = &mtop->ffparams;
3829 functype = ffparams->functype;
3830 ip = ffparams->iparams;
3833 for (i = 0; i < ffparams->ntypes; i++)
3835 ftype = functype[i];
3836 if (ftype == F_DISRES)
3838 label = ip[i].disres.label;
3839 if (label == old_label)
3841 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3849 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3850 "probably the parameters for multiple pairs in one restraint "
3851 "are not identical\n", ndouble);
3856 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3857 gmx_bool posres_only,
3861 gmx_mtop_ilistloop_t iloop;
3871 for (d = 0; d < DIM; d++)
3873 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3875 /* Check for freeze groups */
3876 for (g = 0; g < ir->opts.ngfrz; g++)
3878 for (d = 0; d < DIM; d++)
3880 if (ir->opts.nFreeze[g][d] != 0)
3888 /* Check for position restraints */
3889 iloop = gmx_mtop_ilistloop_init(sys);
3890 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3893 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3895 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3897 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3898 for (d = 0; d < DIM; d++)
3900 if (pr->posres.fcA[d] != 0)
3906 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3908 /* Check for flat-bottom posres */
3909 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3910 if (pr->fbposres.k != 0)
3912 switch (pr->fbposres.geom)
3914 case efbposresSPHERE:
3915 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3917 case efbposresCYLINDERX:
3918 AbsRef[YY] = AbsRef[ZZ] = 1;
3920 case efbposresCYLINDERY:
3921 AbsRef[XX] = AbsRef[ZZ] = 1;
3923 case efbposresCYLINDER:
3924 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3925 case efbposresCYLINDERZ:
3926 AbsRef[XX] = AbsRef[YY] = 1;
3928 case efbposresX: /* d=XX */
3929 case efbposresY: /* d=YY */
3930 case efbposresZ: /* d=ZZ */
3931 d = pr->fbposres.geom - efbposresX;
3935 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3936 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3944 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3948 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3949 gmx_bool *bC6ParametersWorkWithGeometricRules,
3950 gmx_bool *bC6ParametersWorkWithLBRules,
3951 gmx_bool *bLBRulesPossible)
3953 int ntypes, tpi, tpj;
3956 double c6i, c6j, c12i, c12j;
3957 double c6, c6_geometric, c6_LB;
3958 double sigmai, sigmaj, epsi, epsj;
3959 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3962 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3963 * force-field floating point parameters.
3966 ptr = getenv("GMX_LJCOMB_TOL");
3970 double gmx_unused canary;
3972 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3974 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
3979 *bC6ParametersWorkWithLBRules = TRUE;
3980 *bC6ParametersWorkWithGeometricRules = TRUE;
3981 bCanDoLBRules = TRUE;
3982 ntypes = mtop->ffparams.atnr;
3983 snew(typecount, ntypes);
3984 gmx_mtop_count_atomtypes(mtop, state, typecount);
3985 *bLBRulesPossible = TRUE;
3986 for (tpi = 0; tpi < ntypes; ++tpi)
3988 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3989 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3990 for (tpj = tpi; tpj < ntypes; ++tpj)
3992 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3993 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3994 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3995 c6_geometric = std::sqrt(c6i * c6j);
3996 if (!gmx_numzero(c6_geometric))
3998 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4000 sigmai = gmx::sixthroot(c12i / c6i);
4001 sigmaj = gmx::sixthroot(c12j / c6j);
4002 epsi = c6i * c6i /(4.0 * c12i);
4003 epsj = c6j * c6j /(4.0 * c12j);
4004 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4008 *bLBRulesPossible = FALSE;
4009 c6_LB = c6_geometric;
4011 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4014 if (FALSE == bCanDoLBRules)
4016 *bC6ParametersWorkWithLBRules = FALSE;
4019 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4021 if (FALSE == bCanDoGeometricRules)
4023 *bC6ParametersWorkWithGeometricRules = FALSE;
4031 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
4034 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4036 check_combination_rule_differences(mtop, 0,
4037 &bC6ParametersWorkWithGeometricRules,
4038 &bC6ParametersWorkWithLBRules,
4040 if (ir->ljpme_combination_rule == eljpmeLB)
4042 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
4044 warning(wi, "You are using arithmetic-geometric combination rules "
4045 "in LJ-PME, but your non-bonded C6 parameters do not "
4046 "follow these rules.");
4051 if (FALSE == bC6ParametersWorkWithGeometricRules)
4053 if (ir->eDispCorr != edispcNO)
4055 warning_note(wi, "You are using geometric combination rules in "
4056 "LJ-PME, but your non-bonded C6 parameters do "
4057 "not follow these rules. "
4058 "This will introduce very small errors in the forces and energies in "
4059 "your simulations. Dispersion correction will correct total energy "
4060 "and/or pressure for isotropic systems, but not forces or surface tensions.");
4064 warning_note(wi, "You are using geometric combination rules in "
4065 "LJ-PME, but your non-bonded C6 parameters do "
4066 "not follow these rules. "
4067 "This will introduce very small errors in the forces and energies in "
4068 "your simulations. If your system is homogeneous, consider using dispersion correction "
4069 "for the total energy and pressure.");
4075 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
4078 char err_buf[STRLEN];
4080 gmx_bool bCharge, bAcc;
4083 gmx_mtop_atomloop_block_t aloopb;
4084 gmx_mtop_atomloop_all_t aloop;
4086 char warn_buf[STRLEN];
4088 set_warning_line(wi, mdparin, -1);
4090 if (ir->cutoff_scheme == ecutsVERLET &&
4091 ir->verletbuf_tol > 0 &&
4093 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
4094 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4096 /* Check if a too small Verlet buffer might potentially
4097 * cause more drift than the thermostat can couple off.
4099 /* Temperature error fraction for warning and suggestion */
4100 const real T_error_warn = 0.002;
4101 const real T_error_suggest = 0.001;
4102 /* For safety: 2 DOF per atom (typical with constraints) */
4103 const real nrdf_at = 2;
4104 real T, tau, max_T_error;
4109 for (i = 0; i < ir->opts.ngtc; i++)
4111 T = std::max(T, ir->opts.ref_t[i]);
4112 tau = std::max(tau, ir->opts.tau_t[i]);
4116 /* This is a worst case estimate of the temperature error,
4117 * assuming perfect buffer estimation and no cancelation
4118 * of errors. The factor 0.5 is because energy distributes
4119 * equally over Ekin and Epot.
4121 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4122 if (max_T_error > T_error_warn)
4124 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.",
4125 ir->verletbuf_tol, T, tau,
4127 100*T_error_suggest,
4128 ir->verletbuf_tol*T_error_suggest/max_T_error);
4129 warning(wi, warn_buf);
4134 if (ETC_ANDERSEN(ir->etc))
4138 for (i = 0; i < ir->opts.ngtc; i++)
4140 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4141 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4142 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
4143 i, ir->opts.tau_t[i]);
4144 CHECK(ir->opts.tau_t[i] < 0);
4147 for (i = 0; i < ir->opts.ngtc; i++)
4149 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
4150 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);
4151 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
4155 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4156 ir->comm_mode == ecmNO &&
4157 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4158 !ETC_ANDERSEN(ir->etc))
4160 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");
4163 /* Check for pressure coupling with absolute position restraints */
4164 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4166 absolute_reference(ir, sys, TRUE, AbsRef);
4168 for (m = 0; m < DIM; m++)
4170 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4172 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4180 aloopb = gmx_mtop_atomloop_block_init(sys);
4182 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4184 if (atom->q != 0 || atom->qB != 0)
4192 if (EEL_FULL(ir->coulombtype))
4195 "You are using full electrostatics treatment %s for a system without charges.\n"
4196 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4197 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4198 warning(wi, err_buf);
4203 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4206 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4207 "You might want to consider using %s electrostatics.\n",
4209 warning_note(wi, err_buf);
4213 /* Check if combination rules used in LJ-PME are the same as in the force field */
4214 if (EVDW_PME(ir->vdwtype))
4216 check_combination_rules(ir, sys, wi);
4219 /* Generalized reaction field */
4220 if (ir->opts.ngtc == 0)
4222 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4224 CHECK(ir->coulombtype == eelGRF);
4228 sprintf(err_buf, "When using coulombtype = %s"
4229 " ref-t for temperature coupling should be > 0",
4231 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4235 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4237 for (m = 0; (m < DIM); m++)
4239 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4248 snew(mgrp, sys->groups.grps[egcACC].nr);
4249 aloop = gmx_mtop_atomloop_all_init(sys);
4251 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4253 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4256 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4258 for (m = 0; (m < DIM); m++)
4260 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4264 for (m = 0; (m < DIM); m++)
4266 if (fabs(acc[m]) > 1e-6)
4268 const char *dim[DIM] = { "X", "Y", "Z" };
4270 "Net Acceleration in %s direction, will %s be corrected\n",
4271 dim[m], ir->nstcomm != 0 ? "" : "not");
4272 if (ir->nstcomm != 0 && m < ndof_com(ir))
4275 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4277 ir->opts.acc[i][m] -= acc[m];
4285 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4286 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4288 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4296 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4298 if (ir->pull->coord[i].group[0] == 0 ||
4299 ir->pull->coord[i].group[1] == 0)
4301 absolute_reference(ir, sys, FALSE, AbsRef);
4302 for (m = 0; m < DIM; m++)
4304 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4306 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.");
4314 for (i = 0; i < 3; i++)
4316 for (m = 0; m <= i; m++)
4318 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4319 ir->deform[i][m] != 0)
4321 for (c = 0; c < ir->pull->ncoord; c++)
4323 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4324 ir->pull->coord[c].vec[m] != 0)
4326 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4337 void double_check(t_inputrec *ir, matrix box,
4338 gmx_bool bHasNormalConstraints,
4339 gmx_bool bHasAnyConstraints,
4343 char warn_buf[STRLEN];
4346 ptr = check_box(ir->ePBC, box);
4349 warning_error(wi, ptr);
4352 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4354 if (ir->shake_tol <= 0.0)
4356 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4358 warning_error(wi, warn_buf);
4362 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4364 /* If we have Lincs constraints: */
4365 if (ir->eI == eiMD && ir->etc == etcNO &&
4366 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4368 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4369 warning_note(wi, warn_buf);
4372 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4374 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4375 warning_note(wi, warn_buf);
4377 if (ir->epc == epcMTTK)
4379 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4383 if (bHasAnyConstraints && ir->epc == epcMTTK)
4385 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4388 if (ir->LincsWarnAngle > 90.0)
4390 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4391 warning(wi, warn_buf);
4392 ir->LincsWarnAngle = 90.0;
4395 if (ir->ePBC != epbcNONE)
4397 if (ir->nstlist == 0)
4399 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4401 if (ir->ns_type == ensGRID)
4403 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4405 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");
4406 warning_error(wi, warn_buf);
4411 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4412 if (2*ir->rlist >= min_size)
4414 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4415 warning_error(wi, warn_buf);
4418 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4425 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4429 real rvdw1, rvdw2, rcoul1, rcoul2;
4430 char warn_buf[STRLEN];
4432 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4436 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4441 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4447 if (rvdw1 + rvdw2 > ir->rlist ||
4448 rcoul1 + rcoul2 > ir->rlist)
4451 "The sum of the two largest charge group radii (%f) "
4452 "is larger than rlist (%f)\n",
4453 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4454 warning(wi, warn_buf);
4458 /* Here we do not use the zero at cut-off macro,
4459 * since user defined interactions might purposely
4460 * not be zero at the cut-off.
4462 if (ir_vdw_is_zero_at_cutoff(ir) &&
4463 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4465 sprintf(warn_buf, "The sum of the two largest charge group "
4466 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4467 "With exact cut-offs, better performance can be "
4468 "obtained with cutoff-scheme = %s, because it "
4469 "does not use charge groups at all.",
4471 ir->rlist, ir->rvdw,
4472 ecutscheme_names[ecutsVERLET]);
4475 warning(wi, warn_buf);
4479 warning_note(wi, warn_buf);
4482 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4483 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4485 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4486 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4488 ir->rlist, ir->rcoulomb,
4489 ecutscheme_names[ecutsVERLET]);
4492 warning(wi, warn_buf);
4496 warning_note(wi, warn_buf);