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/awh/read-params.h"
51 #include "gromacs/fileio/readinp.h"
52 #include "gromacs/fileio/warninp.h"
53 #include "gromacs/gmxlib/chargegroup.h"
54 #include "gromacs/gmxlib/network.h"
55 #include "gromacs/gmxpreprocess/keyvaluetreemdpwriter.h"
56 #include "gromacs/gmxpreprocess/toputil.h"
57 #include "gromacs/math/functions.h"
58 #include "gromacs/math/units.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/mdlib/calc_verletbuf.h"
61 #include "gromacs/mdrunutility/mdmodules.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/mdtypes/pull-params.h"
65 #include "gromacs/options/options.h"
66 #include "gromacs/options/treesupport.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/topology/block.h"
69 #include "gromacs/topology/ifunc.h"
70 #include "gromacs/topology/index.h"
71 #include "gromacs/topology/mtop_util.h"
72 #include "gromacs/topology/symtab.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/exceptions.h"
76 #include "gromacs/utility/fatalerror.h"
77 #include "gromacs/utility/filestream.h"
78 #include "gromacs/utility/gmxassert.h"
79 #include "gromacs/utility/ikeyvaluetreeerror.h"
80 #include "gromacs/utility/keyvaluetree.h"
81 #include "gromacs/utility/keyvaluetreebuilder.h"
82 #include "gromacs/utility/keyvaluetreetransform.h"
83 #include "gromacs/utility/smalloc.h"
84 #include "gromacs/utility/stringcompare.h"
85 #include "gromacs/utility/stringutil.h"
86 #include "gromacs/utility/textwriter.h"
91 /* Resource parameters
92 * Do not change any of these until you read the instruction
93 * in readinp.h. Some cpp's do not take spaces after the backslash
94 * (like the c-shell), which will give you a very weird compiler
98 typedef struct t_inputrec_strings
100 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
101 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
102 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
103 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
104 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
106 char fep_lambda[efptNR][STRLEN];
107 char lambda_weights[STRLEN];
110 char anneal[STRLEN], anneal_npoints[STRLEN],
111 anneal_time[STRLEN], anneal_temp[STRLEN];
112 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
113 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
114 SAoff[STRLEN], SAsteps[STRLEN];
116 } gmx_inputrec_strings;
118 static gmx_inputrec_strings *is = nullptr;
120 void init_inputrec_strings()
124 gmx_incons("Attempted to call init_inputrec_strings before calling done_inputrec_strings. Only one inputrec (i.e. .mdp file) can be parsed at a time.");
129 void done_inputrec_strings()
137 egrptpALL, /* All particles have to be a member of a group. */
138 egrptpALL_GENREST, /* A rest group with name is generated for particles *
139 * that are not part of any group. */
140 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
141 * for the rest group. */
142 egrptpONE /* Merge all selected groups into one group, *
143 * make a rest group for the remaining particles. */
146 static const char *constraints[eshNR+1] = {
147 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
150 static const char *couple_lam[ecouplamNR+1] = {
151 "vdw-q", "vdw", "q", "none", nullptr
154 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
159 for (i = 0; i < ntemps; i++)
161 /* simple linear scaling -- allows more control */
162 if (simtemp->eSimTempScale == esimtempLINEAR)
164 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
166 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
168 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
170 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
172 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
177 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
178 gmx_fatal(FARGS, errorstr);
185 static void _low_check(gmx_bool b, const char *s, warninp_t wi)
189 warning_error(wi, s);
193 static void check_nst(const char *desc_nst, int nst,
194 const char *desc_p, int *p,
199 if (*p > 0 && *p % nst != 0)
201 /* Round up to the next multiple of nst */
202 *p = ((*p)/nst + 1)*nst;
203 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
204 desc_p, desc_nst, desc_p, *p);
209 static gmx_bool ir_NVE(const t_inputrec *ir)
211 return (EI_MD(ir->eI) && ir->etc == etcNO);
214 static int lcd(int n1, int n2)
219 for (i = 2; (i <= n1 && i <= n2); i++)
221 if (n1 % i == 0 && n2 % i == 0)
230 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
232 if (*eintmod == eintmodPOTSHIFT_VERLET)
234 if (ir->cutoff_scheme == ecutsVERLET)
236 *eintmod = eintmodPOTSHIFT;
240 *eintmod = eintmodNONE;
245 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
247 /* Check internal consistency.
248 * NOTE: index groups are not set here yet, don't check things
249 * like temperature coupling group options here, but in triple_check
252 /* Strange macro: first one fills the err_buf, and then one can check
253 * the condition, which will print the message and increase the error
256 #define CHECK(b) _low_check(b, err_buf, wi)
257 char err_buf[256], warn_buf[STRLEN];
260 t_lambda *fep = ir->fepvals;
261 t_expanded *expand = ir->expandedvals;
263 set_warning_line(wi, mdparin, -1);
265 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
267 sprintf(warn_buf, "%s electrostatics is no longer supported",
268 eel_names[eelRF_NEC_UNSUPPORTED]);
269 warning_error(wi, warn_buf);
272 /* BASIC CUT-OFF STUFF */
273 if (ir->rcoulomb < 0)
275 warning_error(wi, "rcoulomb should be >= 0");
279 warning_error(wi, "rvdw should be >= 0");
282 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
284 warning_error(wi, "rlist should be >= 0");
286 sprintf(err_buf, "nstlist can not be smaller than 0. (If you were trying to use the heuristic neighbour-list update scheme for efficient buffering for improved energy conservation, please use the Verlet cut-off scheme instead.)");
287 CHECK(ir->nstlist < 0);
289 process_interaction_modifier(ir, &ir->coulomb_modifier);
290 process_interaction_modifier(ir, &ir->vdw_modifier);
292 if (ir->cutoff_scheme == ecutsGROUP)
295 "The group cutoff scheme is deprecated since GROMACS 5.0 and will be removed in a future "
296 "release when all interaction forms are supported for the verlet scheme. The verlet "
297 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
299 if (ir->rlist > 0 && ir->rlist < ir->rcoulomb)
301 gmx_fatal(FARGS, "rcoulomb must not be greater than rlist (twin-range schemes are not supported)");
303 if (ir->rlist > 0 && ir->rlist < ir->rvdw)
305 gmx_fatal(FARGS, "rvdw must not be greater than rlist (twin-range schemes are not supported)");
308 if (ir->rlist == 0 && ir->ePBC != epbcNONE)
310 warning_error(wi, "Can not have an infinite cut-off with PBC");
314 if (ir->cutoff_scheme == ecutsVERLET)
318 /* Normal Verlet type neighbor-list, currently only limited feature support */
319 if (inputrec2nboundeddim(ir) < 3)
321 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
324 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
325 if (ir->rcoulomb != ir->rvdw)
327 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
328 // for PME load balancing, we can support this exception.
329 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
330 ir->vdwtype == evdwCUT &&
331 ir->rcoulomb > ir->rvdw);
332 if (!bUsesPmeTwinRangeKernel)
334 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
338 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
340 if (ir->vdw_modifier == eintmodNONE ||
341 ir->vdw_modifier == eintmodPOTSHIFT)
343 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
345 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]);
346 warning_note(wi, warn_buf);
348 ir->vdwtype = evdwCUT;
352 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
353 warning_error(wi, warn_buf);
357 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
359 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
361 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
362 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
364 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
366 if (!(ir->coulomb_modifier == eintmodNONE ||
367 ir->coulomb_modifier == eintmodPOTSHIFT))
369 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
370 warning_error(wi, warn_buf);
373 if (ir->implicit_solvent != eisNO)
375 warning_error(wi, "Implicit solvent is not (yet) supported with the with Verlet lists.");
378 if (EEL_USER(ir->coulombtype))
380 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme", eel_names[ir->coulombtype]);
381 warning_error(wi, warn_buf);
384 if (ir->nstlist <= 0)
386 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
389 if (ir->nstlist < 10)
391 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.");
394 rc_max = std::max(ir->rvdw, ir->rcoulomb);
396 if (ir->verletbuf_tol <= 0)
398 if (ir->verletbuf_tol == 0)
400 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
403 if (ir->rlist < rc_max)
405 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
408 if (ir->rlist == rc_max && ir->nstlist > 1)
410 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.");
415 if (ir->rlist > rc_max)
417 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.");
420 if (ir->nstlist == 1)
422 /* No buffer required */
427 if (EI_DYNAMICS(ir->eI))
429 if (inputrec2nboundeddim(ir) < 3)
431 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.");
433 /* Set rlist temporarily so we can continue processing */
438 /* Set the buffer to 5% of the cut-off */
439 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
445 /* GENERAL INTEGRATOR STUFF */
448 if (ir->etc != etcNO)
450 if (EI_RANDOM(ir->eI))
452 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]);
456 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
458 warning_note(wi, warn_buf);
462 if (ir->eI == eiVVAK)
464 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]);
465 warning_note(wi, warn_buf);
467 if (!EI_DYNAMICS(ir->eI))
469 if (ir->epc != epcNO)
471 sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
472 warning_note(wi, warn_buf);
476 if (EI_DYNAMICS(ir->eI))
478 if (ir->nstcalcenergy < 0)
480 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
481 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
483 /* nstcalcenergy larger than nstener does not make sense.
484 * We ideally want nstcalcenergy=nstener.
488 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
492 ir->nstcalcenergy = ir->nstenergy;
496 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
497 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
498 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
501 const char *nsten = "nstenergy";
502 const char *nstdh = "nstdhdl";
503 const char *min_name = nsten;
504 int min_nst = ir->nstenergy;
506 /* find the smallest of ( nstenergy, nstdhdl ) */
507 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
508 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
510 min_nst = ir->fepvals->nstdhdl;
513 /* If the user sets nstenergy small, we should respect that */
515 "Setting nstcalcenergy (%d) equal to %s (%d)",
516 ir->nstcalcenergy, min_name, min_nst);
517 warning_note(wi, warn_buf);
518 ir->nstcalcenergy = min_nst;
521 if (ir->epc != epcNO)
523 if (ir->nstpcouple < 0)
525 ir->nstpcouple = ir_optimal_nstpcouple(ir);
529 if (ir->nstcalcenergy > 0)
531 if (ir->efep != efepNO)
533 /* nstdhdl should be a multiple of nstcalcenergy */
534 check_nst("nstcalcenergy", ir->nstcalcenergy,
535 "nstdhdl", &ir->fepvals->nstdhdl, wi);
536 /* nstexpanded should be a multiple of nstcalcenergy */
537 check_nst("nstcalcenergy", ir->nstcalcenergy,
538 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
540 /* for storing exact averages nstenergy should be
541 * a multiple of nstcalcenergy
543 check_nst("nstcalcenergy", ir->nstcalcenergy,
544 "nstenergy", &ir->nstenergy, wi);
548 if (ir->nsteps == 0 && !ir->bContinuation)
550 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
554 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
555 ir->bContinuation && ir->ld_seed != -1)
557 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)");
563 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
564 CHECK(ir->ePBC != epbcXYZ);
565 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
566 CHECK(ir->ns_type != ensGRID);
567 sprintf(err_buf, "with TPI nstlist should be larger than zero");
568 CHECK(ir->nstlist <= 0);
569 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
570 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
571 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
572 CHECK(ir->cutoff_scheme == ecutsVERLET);
576 if ( (opts->nshake > 0) && (opts->bMorse) )
579 "Using morse bond-potentials while constraining bonds is useless");
580 warning(wi, warn_buf);
583 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
584 ir->bContinuation && ir->ld_seed != -1)
586 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)");
588 /* verify simulated tempering options */
592 gmx_bool bAllTempZero = TRUE;
593 for (i = 0; i < fep->n_lambda; i++)
595 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]);
596 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
597 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
599 bAllTempZero = FALSE;
602 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
603 CHECK(bAllTempZero == TRUE);
605 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
606 CHECK(ir->eI != eiVV);
608 /* check compatability of the temperature coupling with simulated tempering */
610 if (ir->etc == etcNOSEHOOVER)
612 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
613 warning_note(wi, warn_buf);
616 /* check that the temperatures make sense */
618 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);
619 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
621 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
622 CHECK(ir->simtempvals->simtemp_high <= 0);
624 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
625 CHECK(ir->simtempvals->simtemp_low <= 0);
628 /* verify free energy options */
630 if (ir->efep != efepNO)
633 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
635 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
637 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
638 (int)fep->sc_r_power);
639 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
641 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
642 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
644 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
645 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
647 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
648 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
650 sprintf(err_buf, "Free-energy not implemented for Ewald");
651 CHECK(ir->coulombtype == eelEWALD);
653 /* check validty of lambda inputs */
654 if (fep->n_lambda == 0)
656 /* Clear output in case of no states:*/
657 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
658 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
662 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
663 CHECK((fep->init_fep_state >= fep->n_lambda));
666 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
667 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
669 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",
670 fep->init_lambda, fep->init_fep_state);
671 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
675 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
679 for (i = 0; i < efptNR; i++)
681 if (fep->separate_dvdl[i])
686 if (n_lambda_terms > 1)
688 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.");
689 warning(wi, warn_buf);
692 if (n_lambda_terms < 2 && fep->n_lambda > 0)
695 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
699 for (j = 0; j < efptNR; j++)
701 for (i = 0; i < fep->n_lambda; i++)
703 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]);
704 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
708 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
710 for (i = 0; i < fep->n_lambda; i++)
712 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],
713 fep->all_lambda[efptCOUL][i]);
714 CHECK((fep->sc_alpha > 0) &&
715 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
716 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
717 ((fep->all_lambda[efptVDW][i] > 0.0) &&
718 (fep->all_lambda[efptVDW][i] < 1.0))));
722 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
724 real sigma, lambda, r_sc;
727 /* Maximum estimate for A and B charges equal with lambda power 1 */
729 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);
730 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.",
732 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
733 warning_note(wi, warn_buf);
736 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
737 be treated differently, but that's the next step */
739 for (i = 0; i < efptNR; i++)
741 for (j = 0; j < fep->n_lambda; j++)
743 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
744 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
749 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
753 /* checking equilibration of weights inputs for validity */
755 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
756 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
757 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
759 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
760 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
761 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
763 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
764 expand->equil_steps, elmceq_names[elmceqSTEPS]);
765 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
767 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
768 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
769 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
771 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
772 expand->equil_ratio, elmceq_names[elmceqRATIO]);
773 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
775 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
776 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
777 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
779 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
780 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
781 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
783 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
784 expand->equil_steps, elmceq_names[elmceqSTEPS]);
785 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
787 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
788 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
789 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
791 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
792 expand->equil_ratio, elmceq_names[elmceqRATIO]);
793 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
795 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
796 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
797 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
799 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
800 CHECK((expand->lmc_repeats <= 0));
801 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
802 CHECK((expand->minvarmin <= 0));
803 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
804 CHECK((expand->c_range < 0));
805 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
806 fep->init_fep_state, expand->lmc_forced_nstart);
807 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
808 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
809 CHECK((expand->lmc_forced_nstart < 0));
810 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
811 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
813 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
814 CHECK((expand->init_wl_delta < 0));
815 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
816 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
817 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
818 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
820 /* if there is no temperature control, we need to specify an MC temperature */
821 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
822 if (expand->nstTij > 0)
824 sprintf(err_buf, "nstlog must be non-zero");
825 CHECK(ir->nstlog == 0);
826 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
827 expand->nstTij, ir->nstlog);
828 CHECK((expand->nstTij % ir->nstlog) != 0);
833 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
834 CHECK(ir->nwall && ir->ePBC != epbcXY);
837 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
839 if (ir->ePBC == epbcNONE)
841 if (ir->epc != epcNO)
843 warning(wi, "Turning off pressure coupling for vacuum system");
849 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
850 epbc_names[ir->ePBC]);
851 CHECK(ir->epc != epcNO);
853 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
854 CHECK(EEL_FULL(ir->coulombtype));
856 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
857 epbc_names[ir->ePBC]);
858 CHECK(ir->eDispCorr != edispcNO);
861 if (ir->rlist == 0.0)
863 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
864 "with coulombtype = %s or coulombtype = %s\n"
865 "without periodic boundary conditions (pbc = %s) and\n"
866 "rcoulomb and rvdw set to zero",
867 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
868 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
869 (ir->ePBC != epbcNONE) ||
870 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
874 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
879 if (ir->nstcomm == 0)
881 ir->comm_mode = ecmNO;
883 if (ir->comm_mode != ecmNO)
887 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");
888 ir->nstcomm = abs(ir->nstcomm);
891 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
893 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
894 ir->nstcomm = ir->nstcalcenergy;
897 if (ir->comm_mode == ecmANGULAR)
899 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
900 CHECK(ir->bPeriodicMols);
901 if (ir->ePBC != epbcNONE)
903 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.");
908 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
910 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.");
913 /* TEMPERATURE COUPLING */
914 if (ir->etc == etcYES)
916 ir->etc = etcBERENDSEN;
917 warning_note(wi, "Old option for temperature coupling given: "
918 "changing \"yes\" to \"Berendsen\"\n");
921 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
923 if (ir->opts.nhchainlength < 1)
925 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
926 ir->opts.nhchainlength = 1;
927 warning(wi, warn_buf);
930 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
932 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
933 ir->opts.nhchainlength = 1;
938 ir->opts.nhchainlength = 0;
941 if (ir->eI == eiVVAK)
943 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
945 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
948 if (ETC_ANDERSEN(ir->etc))
950 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
951 CHECK(!(EI_VV(ir->eI)));
953 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
955 sprintf(warn_buf, "Center of mass removal not necessary for %s. All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.", etcoupl_names[ir->etc]);
956 warning_note(wi, warn_buf);
959 sprintf(err_buf, "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step", ir->nstcomm, etcoupl_names[ir->etc]);
960 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
963 if (ir->etc == etcBERENDSEN)
965 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
966 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
967 warning_note(wi, warn_buf);
970 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
971 && ir->epc == epcBERENDSEN)
973 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
974 "true ensemble for the thermostat");
975 warning(wi, warn_buf);
978 /* PRESSURE COUPLING */
979 if (ir->epc == epcISOTROPIC)
981 ir->epc = epcBERENDSEN;
982 warning_note(wi, "Old option for pressure coupling given: "
983 "changing \"Isotropic\" to \"Berendsen\"\n");
986 if (ir->epc != epcNO)
988 dt_pcoupl = ir->nstpcouple*ir->delta_t;
990 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
991 CHECK(ir->tau_p <= 0);
993 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
995 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
996 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
997 warning(wi, warn_buf);
1000 sprintf(err_buf, "compressibility must be > 0 when using pressure"
1001 " coupling %s\n", EPCOUPLTYPE(ir->epc));
1002 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
1003 ir->compress[ZZ][ZZ] < 0 ||
1004 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
1005 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1007 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1010 "You are generating velocities so I am assuming you "
1011 "are equilibrating a system. You are using "
1012 "%s pressure coupling, but this can be "
1013 "unstable for equilibration. If your system crashes, try "
1014 "equilibrating first with Berendsen pressure coupling. If "
1015 "you are not equilibrating the system, you can probably "
1016 "ignore this warning.",
1017 epcoupl_names[ir->epc]);
1018 warning(wi, warn_buf);
1024 if (ir->epc > epcNO)
1026 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1028 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.");
1034 if (ir->epc == epcMTTK)
1036 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1040 /* ELECTROSTATICS */
1041 /* More checks are in triple check (grompp.c) */
1043 if (ir->coulombtype == eelSWITCH)
1045 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1046 "artifacts, advice: use coulombtype = %s",
1047 eel_names[ir->coulombtype],
1048 eel_names[eelRF_ZERO]);
1049 warning(wi, warn_buf);
1052 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1054 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1055 warning_note(wi, warn_buf);
1058 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1060 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);
1061 warning(wi, warn_buf);
1062 ir->epsilon_rf = ir->epsilon_r;
1063 ir->epsilon_r = 1.0;
1066 if (ir->epsilon_r == 0)
1069 "It is pointless to use long-range or Generalized Born electrostatics with infinite relative permittivity."
1070 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1071 CHECK(EEL_FULL(ir->coulombtype) || ir->implicit_solvent == eisGBSA);
1074 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1076 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1077 CHECK(ir->epsilon_r < 0);
1080 if (EEL_RF(ir->coulombtype))
1082 /* reaction field (at the cut-off) */
1084 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1086 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1087 eel_names[ir->coulombtype]);
1088 warning(wi, warn_buf);
1089 ir->epsilon_rf = 0.0;
1092 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1093 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1094 (ir->epsilon_r == 0));
1095 if (ir->epsilon_rf == ir->epsilon_r)
1097 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1098 eel_names[ir->coulombtype]);
1099 warning(wi, warn_buf);
1102 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1103 * means the interaction is zero outside rcoulomb, but it helps to
1104 * provide accurate energy conservation.
1106 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1108 if (ir_coulomb_switched(ir))
1111 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1112 eel_names[ir->coulombtype]);
1113 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1116 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1118 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1120 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1121 eel_names[ir->coulombtype]);
1122 CHECK(ir->rlist > ir->rcoulomb);
1126 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1129 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1130 CHECK( ir->coulomb_modifier != eintmodNONE);
1132 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1135 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1136 CHECK( ir->vdw_modifier != eintmodNONE);
1139 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1140 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1143 "The switch/shift interaction settings are just for compatibility; you will get better "
1144 "performance from applying potential modifiers to your interactions!\n");
1145 warning_note(wi, warn_buf);
1148 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1150 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1152 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1153 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.",
1154 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1155 warning(wi, warn_buf);
1159 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1161 if (ir->rvdw_switch == 0)
1163 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.");
1164 warning(wi, warn_buf);
1168 if (EEL_FULL(ir->coulombtype))
1170 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1171 ir->coulombtype == eelPMEUSERSWITCH)
1173 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1174 eel_names[ir->coulombtype]);
1175 CHECK(ir->rcoulomb > ir->rlist);
1177 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1179 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1182 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist.\n"
1183 "For optimal energy conservation,consider using\n"
1184 "a potential modifier.", eel_names[ir->coulombtype]);
1185 CHECK(ir->rcoulomb != ir->rlist);
1190 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1192 // TODO: Move these checks into the ewald module with the options class
1194 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1196 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1198 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
1199 warning_error(wi, warn_buf);
1203 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1205 if (ir->ewald_geometry == eewg3D)
1207 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1208 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1209 warning(wi, warn_buf);
1211 /* This check avoids extra pbc coding for exclusion corrections */
1212 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1213 CHECK(ir->wall_ewald_zfac < 2);
1215 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1216 EEL_FULL(ir->coulombtype))
1218 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1219 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1220 warning(wi, warn_buf);
1222 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1224 if (ir->cutoff_scheme == ecutsVERLET)
1226 sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
1227 eel_names[ir->coulombtype]);
1228 warning(wi, warn_buf);
1232 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",
1233 eel_names[ir->coulombtype]);
1234 warning_note(wi, warn_buf);
1238 if (ir_vdw_switched(ir))
1240 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1241 CHECK(ir->rvdw_switch >= ir->rvdw);
1243 if (ir->rvdw_switch < 0.5*ir->rvdw)
1245 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.",
1246 ir->rvdw_switch, ir->rvdw);
1247 warning_note(wi, warn_buf);
1250 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1252 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1254 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1255 CHECK(ir->rlist > ir->rvdw);
1259 if (ir->vdwtype == evdwPME)
1261 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1263 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1264 evdw_names[ir->vdwtype],
1265 eintmod_names[eintmodPOTSHIFT],
1266 eintmod_names[eintmodNONE]);
1267 warning_error(wi, err_buf);
1271 if (ir->cutoff_scheme == ecutsGROUP)
1273 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1274 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)))
1276 warning_note(wi, "With exact cut-offs, rlist should be "
1277 "larger than rcoulomb and rvdw, so that there "
1278 "is a buffer region for particle motion "
1279 "between neighborsearch steps");
1282 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlist <= ir->rcoulomb)
1284 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rcoulomb.");
1285 warning_note(wi, warn_buf);
1287 if (ir_vdw_switched(ir) && (ir->rlist <= ir->rvdw))
1289 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rvdw.");
1290 warning_note(wi, warn_buf);
1294 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1296 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.");
1299 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1302 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1305 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1307 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1310 /* ENERGY CONSERVATION */
1311 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1313 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1315 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1316 evdw_names[evdwSHIFT]);
1317 warning_note(wi, warn_buf);
1319 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1321 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1322 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1323 warning_note(wi, warn_buf);
1327 /* IMPLICIT SOLVENT */
1328 if (ir->coulombtype == eelGB_NOTUSED)
1330 sprintf(warn_buf, "Invalid option %s for coulombtype",
1331 eel_names[ir->coulombtype]);
1332 warning_error(wi, warn_buf);
1335 if (ir->sa_algorithm == esaSTILL)
1337 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1338 CHECK(ir->sa_algorithm == esaSTILL);
1341 if (ir->implicit_solvent == eisGBSA)
1343 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1344 CHECK(ir->rgbradii != ir->rlist);
1346 if (ir->coulombtype != eelCUT)
1348 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1349 CHECK(ir->coulombtype != eelCUT);
1351 if (ir->vdwtype != evdwCUT)
1353 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1354 CHECK(ir->vdwtype != evdwCUT);
1356 if (ir->nstgbradii < 1)
1358 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1359 warning_note(wi, warn_buf);
1362 if (ir->sa_algorithm == esaNO)
1364 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1365 warning_note(wi, warn_buf);
1367 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1369 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");
1370 warning_note(wi, warn_buf);
1372 if (ir->gb_algorithm == egbSTILL)
1374 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1378 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1381 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1383 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1384 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1391 if (ir->cutoff_scheme != ecutsGROUP)
1393 warning_error(wi, "QMMM is currently only supported with cutoff-scheme=group");
1395 if (!EI_DYNAMICS(ir->eI))
1398 sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
1399 warning_error(wi, buf);
1405 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1409 /* count the number of text elemets separated by whitespace in a string.
1410 str = the input string
1411 maxptr = the maximum number of allowed elements
1412 ptr = the output array of pointers to the first character of each element
1413 returns: the number of elements. */
1414 int str_nelem(const char *str, int maxptr, char *ptr[])
1419 copy0 = gmx_strdup(str);
1422 while (*copy != '\0')
1426 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1434 while ((*copy != '\0') && !isspace(*copy))
1453 /* interpret a number of doubles from a string and put them in an array,
1454 after allocating space for them.
1455 str = the input string
1456 n = the (pre-allocated) number of doubles read
1457 r = the output array of doubles. */
1458 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1463 char warn_buf[STRLEN];
1465 *n = str_nelem(str, MAXPTR, ptr);
1468 for (i = 0; i < *n; i++)
1470 (*r)[i] = strtod(ptr[i], &endptr);
1473 sprintf(warn_buf, "Invalid value %s in string in mdp file. Expected a real number.", ptr[i]);
1474 warning_error(wi, warn_buf);
1479 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1482 int i, j, max_n_lambda, nweights, nfep[efptNR];
1483 t_lambda *fep = ir->fepvals;
1484 t_expanded *expand = ir->expandedvals;
1485 real **count_fep_lambdas;
1486 gmx_bool bOneLambda = TRUE;
1488 snew(count_fep_lambdas, efptNR);
1490 /* FEP input processing */
1491 /* first, identify the number of lambda values for each type.
1492 All that are nonzero must have the same number */
1494 for (i = 0; i < efptNR; i++)
1496 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1499 /* now, determine the number of components. All must be either zero, or equal. */
1502 for (i = 0; i < efptNR; i++)
1504 if (nfep[i] > max_n_lambda)
1506 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1507 must have the same number if its not zero.*/
1512 for (i = 0; i < efptNR; i++)
1516 ir->fepvals->separate_dvdl[i] = FALSE;
1518 else if (nfep[i] == max_n_lambda)
1520 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1521 respect to the temperature currently */
1523 ir->fepvals->separate_dvdl[i] = TRUE;
1528 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1529 nfep[i], efpt_names[i], max_n_lambda);
1532 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1533 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1535 /* the number of lambdas is the number we've read in, which is either zero
1536 or the same for all */
1537 fep->n_lambda = max_n_lambda;
1539 /* allocate space for the array of lambda values */
1540 snew(fep->all_lambda, efptNR);
1541 /* if init_lambda is defined, we need to set lambda */
1542 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1544 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1546 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1547 for (i = 0; i < efptNR; i++)
1549 snew(fep->all_lambda[i], fep->n_lambda);
1550 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1553 for (j = 0; j < fep->n_lambda; j++)
1555 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1557 sfree(count_fep_lambdas[i]);
1560 sfree(count_fep_lambdas);
1562 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1563 bookkeeping -- for now, init_lambda */
1565 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1567 for (i = 0; i < fep->n_lambda; i++)
1569 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1573 /* check to see if only a single component lambda is defined, and soft core is defined.
1574 In this case, turn on coulomb soft core */
1576 if (max_n_lambda == 0)
1582 for (i = 0; i < efptNR; i++)
1584 if ((nfep[i] != 0) && (i != efptFEP))
1590 if ((bOneLambda) && (fep->sc_alpha > 0))
1592 fep->bScCoul = TRUE;
1595 /* Fill in the others with the efptFEP if they are not explicitly
1596 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1597 they are all zero. */
1599 for (i = 0; i < efptNR; i++)
1601 if ((nfep[i] == 0) && (i != efptFEP))
1603 for (j = 0; j < fep->n_lambda; j++)
1605 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1611 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1612 if (fep->sc_r_power == 48)
1614 if (fep->sc_alpha > 0.1)
1616 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1620 /* now read in the weights */
1621 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1624 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1626 else if (nweights != fep->n_lambda)
1628 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1629 nweights, fep->n_lambda);
1631 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1633 expand->nstexpanded = fep->nstdhdl;
1634 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1636 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1638 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1639 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1640 2*tau_t just to be careful so it's not to frequent */
1645 static void do_simtemp_params(t_inputrec *ir)
1648 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1649 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1654 static void do_wall_params(t_inputrec *ir,
1655 char *wall_atomtype, char *wall_density,
1659 char *names[MAXPTR];
1662 opts->wall_atomtype[0] = nullptr;
1663 opts->wall_atomtype[1] = nullptr;
1665 ir->wall_atomtype[0] = -1;
1666 ir->wall_atomtype[1] = -1;
1667 ir->wall_density[0] = 0;
1668 ir->wall_density[1] = 0;
1672 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1673 if (nstr != ir->nwall)
1675 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1678 for (i = 0; i < ir->nwall; i++)
1680 opts->wall_atomtype[i] = gmx_strdup(names[i]);
1683 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1685 nstr = str_nelem(wall_density, MAXPTR, names);
1686 if (nstr != ir->nwall)
1688 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1690 for (i = 0; i < ir->nwall; i++)
1692 if (sscanf(names[i], "%lf", &dbl) != 1)
1694 gmx_fatal(FARGS, "Could not parse wall-density value from string '%s'", names[i]);
1698 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1700 ir->wall_density[i] = dbl;
1706 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1714 srenew(groups->grpname, groups->ngrpname+nwall);
1715 grps = &(groups->grps[egcENER]);
1716 srenew(grps->nm_ind, grps->nr+nwall);
1717 for (i = 0; i < nwall; i++)
1719 sprintf(str, "wall%d", i);
1720 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1721 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1726 static void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1727 t_expanded *expand, warninp_t wi)
1735 /* read expanded ensemble parameters */
1736 CCTYPE ("expanded ensemble variables");
1737 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1738 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1739 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1740 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1741 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1742 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1743 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1744 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1745 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1746 CCTYPE("Seed for Monte Carlo in lambda space");
1747 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1748 RTYPE ("mc-temperature", expand->mc_temp, -1);
1749 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1750 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1751 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1752 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1753 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1754 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1755 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1756 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1757 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1758 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1759 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1767 /*! \brief Return whether an end state with the given coupling-lambda
1768 * value describes fully-interacting VDW.
1770 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1771 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1773 static gmx_bool couple_lambda_has_vdw_on(int couple_lambda_value)
1775 return (couple_lambda_value == ecouplamVDW ||
1776 couple_lambda_value == ecouplamVDWQ);
1782 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1785 explicit MdpErrorHandler(warninp_t wi)
1786 : wi_(wi), mapping_(nullptr)
1790 void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
1792 mapping_ = &mapping;
1795 virtual bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context)
1797 ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
1798 getOptionName(context).c_str()));
1799 std::string message = gmx::formatExceptionMessageToString(*ex);
1800 warning_error(wi_, message.c_str());
1805 std::string getOptionName(const gmx::KeyValueTreePath &context)
1807 if (mapping_ != nullptr)
1809 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1810 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1813 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1818 const gmx::IKeyValueTreeBackMapping *mapping_;
1823 void get_ir(const char *mdparin, const char *mdparout,
1824 gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
1825 WriteMdpHeader writeMdpHeader, warninp_t wi)
1828 double dumdub[2][6];
1832 char warn_buf[STRLEN];
1833 t_lambda *fep = ir->fepvals;
1834 t_expanded *expand = ir->expandedvals;
1836 init_inputrec_strings();
1837 gmx::TextInputFile stream(mdparin);
1838 inp = read_inpfile(&stream, mdparin, &ninp, wi);
1840 snew(dumstr[0], STRLEN);
1841 snew(dumstr[1], STRLEN);
1843 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1846 "%s did not specify a value for the .mdp option "
1847 "\"cutoff-scheme\". Probably it was first intended for use "
1848 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1849 "introduced, but the group scheme was still the default. "
1850 "The default is now the Verlet scheme, so you will observe "
1851 "different behaviour.", mdparin);
1852 warning_note(wi, warn_buf);
1855 /* ignore the following deprecated commands */
1858 REM_TYPE("domain-decomposition");
1859 REM_TYPE("andersen-seed");
1861 REM_TYPE("dihre-fc");
1862 REM_TYPE("dihre-tau");
1863 REM_TYPE("nstdihreout");
1864 REM_TYPE("nstcheckpoint");
1865 REM_TYPE("optimize-fft");
1866 REM_TYPE("adress_type");
1867 REM_TYPE("adress_const_wf");
1868 REM_TYPE("adress_ex_width");
1869 REM_TYPE("adress_hy_width");
1870 REM_TYPE("adress_ex_forcecap");
1871 REM_TYPE("adress_interface_correction");
1872 REM_TYPE("adress_site");
1873 REM_TYPE("adress_reference_coords");
1874 REM_TYPE("adress_tf_grp_names");
1875 REM_TYPE("adress_cg_grp_names");
1876 REM_TYPE("adress_do_hybridpairs");
1877 REM_TYPE("rlistlong");
1878 REM_TYPE("nstcalclr");
1879 REM_TYPE("pull-print-com2");
1881 /* replace the following commands with the clearer new versions*/
1882 REPL_TYPE("unconstrained-start", "continuation");
1883 REPL_TYPE("foreign-lambda", "fep-lambdas");
1884 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1885 REPL_TYPE("nstxtcout", "nstxout-compressed");
1886 REPL_TYPE("xtc-grps", "compressed-x-grps");
1887 REPL_TYPE("xtc-precision", "compressed-x-precision");
1888 REPL_TYPE("pull-print-com1", "pull-print-com");
1890 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1891 CTYPE ("Preprocessor information: use cpp syntax.");
1892 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1893 STYPE ("include", opts->include, nullptr);
1894 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1895 STYPE ("define", opts->define, nullptr);
1897 CCTYPE ("RUN CONTROL PARAMETERS");
1898 EETYPE("integrator", ir->eI, ei_names);
1899 CTYPE ("Start time and timestep in ps");
1900 RTYPE ("tinit", ir->init_t, 0.0);
1901 RTYPE ("dt", ir->delta_t, 0.001);
1902 STEPTYPE ("nsteps", ir->nsteps, 0);
1903 CTYPE ("For exact run continuation or redoing part of a run");
1904 STEPTYPE ("init-step", ir->init_step, 0);
1905 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1906 ITYPE ("simulation-part", ir->simulation_part, 1);
1907 CTYPE ("mode for center of mass motion removal");
1908 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1909 CTYPE ("number of steps for center of mass motion removal");
1910 ITYPE ("nstcomm", ir->nstcomm, 100);
1911 CTYPE ("group(s) for center of mass motion removal");
1912 STYPE ("comm-grps", is->vcm, nullptr);
1914 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1915 CTYPE ("Friction coefficient (amu/ps) and random seed");
1916 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1917 STEPTYPE ("ld-seed", ir->ld_seed, -1);
1920 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1921 CTYPE ("Force tolerance and initial step-size");
1922 RTYPE ("emtol", ir->em_tol, 10.0);
1923 RTYPE ("emstep", ir->em_stepsize, 0.01);
1924 CTYPE ("Max number of iterations in relax-shells");
1925 ITYPE ("niter", ir->niter, 20);
1926 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1927 RTYPE ("fcstep", ir->fc_stepsize, 0);
1928 CTYPE ("Frequency of steepest descents steps when doing CG");
1929 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1930 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1932 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1933 RTYPE ("rtpi", ir->rtpi, 0.05);
1935 /* Output options */
1936 CCTYPE ("OUTPUT CONTROL OPTIONS");
1937 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1938 ITYPE ("nstxout", ir->nstxout, 0);
1939 ITYPE ("nstvout", ir->nstvout, 0);
1940 ITYPE ("nstfout", ir->nstfout, 0);
1941 CTYPE ("Output frequency for energies to log file and energy file");
1942 ITYPE ("nstlog", ir->nstlog, 1000);
1943 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1944 ITYPE ("nstenergy", ir->nstenergy, 1000);
1945 CTYPE ("Output frequency and precision for .xtc file");
1946 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
1947 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1948 CTYPE ("This selects the subset of atoms for the compressed");
1949 CTYPE ("trajectory file. You can select multiple groups. By");
1950 CTYPE ("default, all atoms will be written.");
1951 STYPE ("compressed-x-grps", is->x_compressed_groups, nullptr);
1952 CTYPE ("Selection of energy groups");
1953 STYPE ("energygrps", is->energy, nullptr);
1955 /* Neighbor searching */
1956 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1957 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1958 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1959 CTYPE ("nblist update frequency");
1960 ITYPE ("nstlist", ir->nstlist, 10);
1961 CTYPE ("ns algorithm (simple or grid)");
1962 EETYPE("ns-type", ir->ns_type, ens_names);
1963 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1964 EETYPE("pbc", ir->ePBC, epbc_names);
1965 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1966 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1967 CTYPE ("a value of -1 means: use rlist");
1968 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1969 CTYPE ("nblist cut-off");
1970 RTYPE ("rlist", ir->rlist, 1.0);
1971 CTYPE ("long-range cut-off for switched potentials");
1973 /* Electrostatics */
1974 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1975 CTYPE ("Method for doing electrostatics");
1976 EETYPE("coulombtype", ir->coulombtype, eel_names);
1977 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1978 CTYPE ("cut-off lengths");
1979 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1980 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1981 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1982 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1983 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1984 CTYPE ("Method for doing Van der Waals");
1985 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1986 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1987 CTYPE ("cut-off lengths");
1988 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1989 RTYPE ("rvdw", ir->rvdw, 1.0);
1990 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1991 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1992 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1993 RTYPE ("table-extension", ir->tabext, 1.0);
1994 CTYPE ("Separate tables between energy group pairs");
1995 STYPE ("energygrp-table", is->egptable, nullptr);
1996 CTYPE ("Spacing for the PME/PPPM FFT grid");
1997 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1998 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1999 ITYPE ("fourier-nx", ir->nkx, 0);
2000 ITYPE ("fourier-ny", ir->nky, 0);
2001 ITYPE ("fourier-nz", ir->nkz, 0);
2002 CTYPE ("EWALD/PME/PPPM parameters");
2003 ITYPE ("pme-order", ir->pme_order, 4);
2004 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
2005 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
2006 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
2007 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
2008 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
2010 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
2011 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
2013 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
2014 CTYPE ("Algorithm for calculating Born radii");
2015 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
2016 CTYPE ("Frequency of calculating the Born radii inside rlist");
2017 ITYPE ("nstgbradii", ir->nstgbradii, 1);
2018 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
2019 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
2020 RTYPE ("rgbradii", ir->rgbradii, 1.0);
2021 CTYPE ("Dielectric coefficient of the implicit solvent");
2022 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
2023 CTYPE ("Salt concentration in M for Generalized Born models");
2024 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
2025 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
2026 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
2027 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
2028 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
2029 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
2030 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
2031 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
2032 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
2033 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
2035 /* Coupling stuff */
2036 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
2037 CTYPE ("Temperature coupling");
2038 EETYPE("tcoupl", ir->etc, etcoupl_names);
2039 ITYPE ("nsttcouple", ir->nsttcouple, -1);
2040 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
2041 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
2042 CTYPE ("Groups to couple separately");
2043 STYPE ("tc-grps", is->tcgrps, nullptr);
2044 CTYPE ("Time constant (ps) and reference temperature (K)");
2045 STYPE ("tau-t", is->tau_t, nullptr);
2046 STYPE ("ref-t", is->ref_t, nullptr);
2047 CTYPE ("pressure coupling");
2048 EETYPE("pcoupl", ir->epc, epcoupl_names);
2049 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
2050 ITYPE ("nstpcouple", ir->nstpcouple, -1);
2051 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
2052 RTYPE ("tau-p", ir->tau_p, 1.0);
2053 STYPE ("compressibility", dumstr[0], nullptr);
2054 STYPE ("ref-p", dumstr[1], nullptr);
2055 CTYPE ("Scaling of reference coordinates, No, All or COM");
2056 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
2059 CCTYPE ("OPTIONS FOR QMMM calculations");
2060 EETYPE("QMMM", ir->bQMMM, yesno_names);
2061 CTYPE ("Groups treated Quantum Mechanically");
2062 STYPE ("QMMM-grps", is->QMMM, nullptr);
2063 CTYPE ("QM method");
2064 STYPE("QMmethod", is->QMmethod, nullptr);
2065 CTYPE ("QMMM scheme");
2066 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
2067 CTYPE ("QM basisset");
2068 STYPE("QMbasis", is->QMbasis, nullptr);
2069 CTYPE ("QM charge");
2070 STYPE ("QMcharge", is->QMcharge, nullptr);
2071 CTYPE ("QM multiplicity");
2072 STYPE ("QMmult", is->QMmult, nullptr);
2073 CTYPE ("Surface Hopping");
2074 STYPE ("SH", is->bSH, nullptr);
2075 CTYPE ("CAS space options");
2076 STYPE ("CASorbitals", is->CASorbitals, nullptr);
2077 STYPE ("CASelectrons", is->CASelectrons, nullptr);
2078 STYPE ("SAon", is->SAon, nullptr);
2079 STYPE ("SAoff", is->SAoff, nullptr);
2080 STYPE ("SAsteps", is->SAsteps, nullptr);
2081 CTYPE ("Scale factor for MM charges");
2082 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
2084 /* Simulated annealing */
2085 CCTYPE("SIMULATED ANNEALING");
2086 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
2087 STYPE ("annealing", is->anneal, nullptr);
2088 CTYPE ("Number of time points to use for specifying annealing in each group");
2089 STYPE ("annealing-npoints", is->anneal_npoints, nullptr);
2090 CTYPE ("List of times at the annealing points for each group");
2091 STYPE ("annealing-time", is->anneal_time, nullptr);
2092 CTYPE ("Temp. at each annealing point, for each group.");
2093 STYPE ("annealing-temp", is->anneal_temp, nullptr);
2096 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
2097 EETYPE("gen-vel", opts->bGenVel, yesno_names);
2098 RTYPE ("gen-temp", opts->tempi, 300.0);
2099 ITYPE ("gen-seed", opts->seed, -1);
2102 CCTYPE ("OPTIONS FOR BONDS");
2103 EETYPE("constraints", opts->nshake, constraints);
2104 CTYPE ("Type of constraint algorithm");
2105 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
2106 CTYPE ("Do not constrain the start configuration");
2107 EETYPE("continuation", ir->bContinuation, yesno_names);
2108 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2109 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2110 CTYPE ("Relative tolerance of shake");
2111 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2112 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2113 ITYPE ("lincs-order", ir->nProjOrder, 4);
2114 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2115 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2116 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2117 ITYPE ("lincs-iter", ir->nLincsIter, 1);
2118 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2119 CTYPE ("rotates over more degrees than");
2120 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2121 CTYPE ("Convert harmonic bonds to morse potentials");
2122 EETYPE("morse", opts->bMorse, yesno_names);
2124 /* Energy group exclusions */
2125 CCTYPE ("ENERGY GROUP EXCLUSIONS");
2126 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2127 STYPE ("energygrp-excl", is->egpexcl, nullptr);
2131 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2132 ITYPE ("nwall", ir->nwall, 0);
2133 EETYPE("wall-type", ir->wall_type, ewt_names);
2134 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2135 STYPE ("wall-atomtype", is->wall_atomtype, nullptr);
2136 STYPE ("wall-density", is->wall_density, nullptr);
2137 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2140 CCTYPE("COM PULLING");
2141 EETYPE("pull", ir->bPull, yesno_names);
2145 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, wi);
2149 NOTE: needs COM pulling input */
2150 CCTYPE("AWH biasing");
2151 EETYPE("awh", ir->bDoAwh, yesno_names);
2156 ir->awhParams = gmx::readAndCheckAwhParams(&ninp, &inp, ir, wi);
2160 gmx_fatal(FARGS, "AWH biasing is only compatible with COM pulling turned on");
2164 /* Enforced rotation */
2165 CCTYPE("ENFORCED ROTATION");
2166 CTYPE("Enforced rotation: No or Yes");
2167 EETYPE("rotation", ir->bRot, yesno_names);
2171 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2174 /* Interactive MD */
2176 CCTYPE("Group to display and/or manipulate in interactive MD session");
2177 STYPE ("IMD-group", is->imd_grp, nullptr);
2178 if (is->imd_grp[0] != '\0')
2185 CCTYPE("NMR refinement stuff");
2186 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2187 EETYPE("disre", ir->eDisre, edisre_names);
2188 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2189 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2190 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2191 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2192 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2193 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2194 CTYPE ("Output frequency for pair distances to energy file");
2195 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2196 CTYPE ("Orientation restraints: No or Yes");
2197 EETYPE("orire", opts->bOrire, yesno_names);
2198 CTYPE ("Orientation restraints force constant and tau for time averaging");
2199 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2200 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2201 STYPE ("orire-fitgrp", is->orirefitgrp, nullptr);
2202 CTYPE ("Output frequency for trace(SD) and S to energy file");
2203 ITYPE ("nstorireout", ir->nstorireout, 100);
2205 /* free energy variables */
2206 CCTYPE ("Free energy variables");
2207 EETYPE("free-energy", ir->efep, efep_names);
2208 STYPE ("couple-moltype", is->couple_moltype, nullptr);
2209 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2210 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2211 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2213 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2215 it was not entered */
2216 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2217 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2218 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2219 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], nullptr);
2220 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], nullptr);
2221 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
2222 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
2223 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
2224 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
2225 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
2226 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2227 STYPE ("init-lambda-weights", is->lambda_weights, nullptr);
2228 EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
2229 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2230 ITYPE ("sc-power", fep->sc_power, 1);
2231 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2232 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2233 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2234 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2235 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2236 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2237 separate_dhdl_file_names);
2238 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2239 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2240 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2242 /* Non-equilibrium MD stuff */
2243 CCTYPE("Non-equilibrium MD stuff");
2244 STYPE ("acc-grps", is->accgrps, nullptr);
2245 STYPE ("accelerate", is->acc, nullptr);
2246 STYPE ("freezegrps", is->freeze, nullptr);
2247 STYPE ("freezedim", is->frdim, nullptr);
2248 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2249 STYPE ("deform", is->deform, nullptr);
2251 /* simulated tempering variables */
2252 CCTYPE("simulated tempering variables");
2253 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2254 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2255 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2256 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2258 /* expanded ensemble variables */
2259 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2261 read_expandedparams(&ninp, &inp, expand, wi);
2264 /* Electric fields */
2266 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(ninp, inp);
2267 gmx::KeyValueTreeTransformer transform;
2268 transform.rules()->addRule()
2269 .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2270 mdModules->initMdpTransform(transform.rules());
2271 for (const auto &path : transform.mappedPaths())
2273 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2274 mark_einp_set(ninp, inp, path[0].c_str());
2276 MdpErrorHandler errorHandler(wi);
2278 = transform.transform(convertedValues, &errorHandler);
2279 ir->params = new gmx::KeyValueTreeObject(result.object());
2280 mdModules->adjustInputrecBasedOnModules(ir);
2281 errorHandler.setBackMapping(result.backMapping());
2282 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2285 /* Ion/water position swapping ("computational electrophysiology") */
2286 CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2287 CTYPE("Swap positions along direction: no, X, Y, Z");
2288 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2289 if (ir->eSwapCoords != eswapNO)
2296 CTYPE("Swap attempt frequency");
2297 ITYPE("swap-frequency", ir->swap->nstswap, 1);
2298 CTYPE("Number of ion types to be controlled");
2299 ITYPE("iontypes", nIonTypes, 1);
2302 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2304 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2305 snew(ir->swap->grp, ir->swap->ngrp);
2306 for (i = 0; i < ir->swap->ngrp; i++)
2308 snew(ir->swap->grp[i].molname, STRLEN);
2310 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2311 STYPE("split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2312 STYPE("split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2313 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2314 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2315 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2317 CTYPE("Name of solvent molecules");
2318 STYPE("solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2320 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2321 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2322 CTYPE("however, if correctly defined, the permeation events are recorded per channel");
2323 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2324 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2325 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2326 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2327 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2328 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2330 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2331 ITYPE("coupl-steps", ir->swap->nAverage, 10);
2333 CTYPE("Names of the ion types that can be exchanged with solvent molecules,");
2334 CTYPE("and the requested number of ions of this type in compartments A and B");
2335 CTYPE("-1 means fix the numbers as found in step 0");
2336 for (i = 0; i < nIonTypes; i++)
2338 int ig = eSwapFixedGrpNR + i;
2340 sprintf(buf, "iontype%d-name", i);
2341 STYPE(buf, ir->swap->grp[ig].molname, nullptr);
2342 sprintf(buf, "iontype%d-in-A", i);
2343 ITYPE(buf, ir->swap->grp[ig].nmolReq[0], -1);
2344 sprintf(buf, "iontype%d-in-B", i);
2345 ITYPE(buf, ir->swap->grp[ig].nmolReq[1], -1);
2348 CTYPE("By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2349 CTYPE("at maximum distance (= bulk concentration) to the split group layers. However,");
2350 CTYPE("an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2351 CTYPE("towards one of the compartment-partitioning layers (at +/- 1.0).");
2352 RTYPE("bulk-offsetA", ir->swap->bulkOffset[0], 0.0);
2353 RTYPE("bulk-offsetB", ir->swap->bulkOffset[1], 0.0);
2354 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2355 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2357 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2360 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2361 RTYPE("threshold", ir->swap->threshold, 1.0);
2364 /* AdResS is no longer supported, but we need mdrun to be able to refuse to run old AdResS .tpr files */
2365 EETYPE("adress", ir->bAdress, yesno_names);
2367 /* User defined thingies */
2368 CCTYPE ("User defined thingies");
2369 STYPE ("user1-grps", is->user1, nullptr);
2370 STYPE ("user2-grps", is->user2, nullptr);
2371 ITYPE ("userint1", ir->userint1, 0);
2372 ITYPE ("userint2", ir->userint2, 0);
2373 ITYPE ("userint3", ir->userint3, 0);
2374 ITYPE ("userint4", ir->userint4, 0);
2375 RTYPE ("userreal1", ir->userreal1, 0);
2376 RTYPE ("userreal2", ir->userreal2, 0);
2377 RTYPE ("userreal3", ir->userreal3, 0);
2378 RTYPE ("userreal4", ir->userreal4, 0);
2382 gmx::TextOutputFile stream(mdparout);
2383 write_inpfile(&stream, mdparout, ninp, inp, FALSE, writeMdpHeader, wi);
2385 // Transform module data into a flat key-value tree for output.
2386 gmx::KeyValueTreeBuilder builder;
2387 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2388 mdModules->buildMdpOutput(&builderObject);
2390 gmx::TextWriter writer(&stream);
2391 writeKeyValueTreeAsMdp(&writer, builder.build());
2396 for (i = 0; (i < ninp); i++)
2399 sfree(inp[i].value);
2403 /* Process options if necessary */
2404 for (m = 0; m < 2; m++)
2406 for (i = 0; i < 2*DIM; i++)
2415 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2417 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
2419 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2421 case epctSEMIISOTROPIC:
2422 case epctSURFACETENSION:
2423 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2425 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
2427 dumdub[m][YY] = dumdub[m][XX];
2429 case epctANISOTROPIC:
2430 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2431 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2432 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2434 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
2438 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2439 epcoupltype_names[ir->epct]);
2443 clear_mat(ir->ref_p);
2444 clear_mat(ir->compress);
2445 for (i = 0; i < DIM; i++)
2447 ir->ref_p[i][i] = dumdub[1][i];
2448 ir->compress[i][i] = dumdub[0][i];
2450 if (ir->epct == epctANISOTROPIC)
2452 ir->ref_p[XX][YY] = dumdub[1][3];
2453 ir->ref_p[XX][ZZ] = dumdub[1][4];
2454 ir->ref_p[YY][ZZ] = dumdub[1][5];
2455 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2457 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2459 ir->compress[XX][YY] = dumdub[0][3];
2460 ir->compress[XX][ZZ] = dumdub[0][4];
2461 ir->compress[YY][ZZ] = dumdub[0][5];
2462 for (i = 0; i < DIM; i++)
2464 for (m = 0; m < i; m++)
2466 ir->ref_p[i][m] = ir->ref_p[m][i];
2467 ir->compress[i][m] = ir->compress[m][i];
2472 if (ir->comm_mode == ecmNO)
2477 opts->couple_moltype = nullptr;
2478 if (strlen(is->couple_moltype) > 0)
2480 if (ir->efep != efepNO)
2482 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2483 if (opts->couple_lam0 == opts->couple_lam1)
2485 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2487 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2488 opts->couple_lam1 == ecouplamNONE))
2490 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2495 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2498 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2499 if (ir->efep != efepNO)
2501 if (fep->delta_lambda > 0)
2503 ir->efep = efepSLOWGROWTH;
2507 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2509 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2510 warning_note(wi, "Old option for dhdl-print-energy given: "
2511 "changing \"yes\" to \"total\"\n");
2514 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2516 /* always print out the energy to dhdl if we are doing
2517 expanded ensemble, since we need the total energy for
2518 analysis if the temperature is changing. In some
2519 conditions one may only want the potential energy, so
2520 we will allow that if the appropriate mdp setting has
2521 been enabled. Otherwise, total it is:
2523 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2526 if ((ir->efep != efepNO) || ir->bSimTemp)
2528 ir->bExpanded = FALSE;
2529 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2531 ir->bExpanded = TRUE;
2533 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2534 if (ir->bSimTemp) /* done after fep params */
2536 do_simtemp_params(ir);
2539 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2540 * setup and not on the old way of specifying the free-energy setup,
2541 * we should check for using soft-core when not needed, since that
2542 * can complicate the sampling significantly.
2543 * Note that we only check for the automated coupling setup.
2544 * If the (advanced) user does FEP through manual topology changes,
2545 * this check will not be triggered.
2547 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2548 ir->fepvals->sc_alpha != 0 &&
2549 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2550 couple_lambda_has_vdw_on(opts->couple_lam1)))
2552 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.");
2557 ir->fepvals->n_lambda = 0;
2560 /* WALL PARAMETERS */
2562 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2564 /* ORIENTATION RESTRAINT PARAMETERS */
2566 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, nullptr) != 1)
2568 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2571 /* DEFORMATION PARAMETERS */
2573 clear_mat(ir->deform);
2574 for (i = 0; i < 6; i++)
2579 double gmx_unused canary;
2580 int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
2581 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2582 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
2584 if (strlen(is->deform) > 0 && ndeform != 6)
2586 warning_error(wi, gmx::formatString("Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform).c_str());
2588 for (i = 0; i < 3; i++)
2590 ir->deform[i][i] = dumdub[0][i];
2592 ir->deform[YY][XX] = dumdub[0][3];
2593 ir->deform[ZZ][XX] = dumdub[0][4];
2594 ir->deform[ZZ][YY] = dumdub[0][5];
2595 if (ir->epc != epcNO)
2597 for (i = 0; i < 3; i++)
2599 for (j = 0; j <= i; j++)
2601 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2603 warning_error(wi, "A box element has deform set and compressibility > 0");
2607 for (i = 0; i < 3; i++)
2609 for (j = 0; j < i; j++)
2611 if (ir->deform[i][j] != 0)
2613 for (m = j; m < DIM; m++)
2615 if (ir->compress[m][j] != 0)
2617 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.");
2618 warning(wi, warn_buf);
2626 /* Ion/water position swapping checks */
2627 if (ir->eSwapCoords != eswapNO)
2629 if (ir->swap->nstswap < 1)
2631 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2633 if (ir->swap->nAverage < 1)
2635 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2637 if (ir->swap->threshold < 1.0)
2639 warning_error(wi, "Ion count threshold must be at least 1.\n");
2647 static int search_QMstring(const char *s, int ng, const char *gn[])
2649 /* same as normal search_string, but this one searches QM strings */
2652 for (i = 0; (i < ng); i++)
2654 if (gmx_strcasecmp(s, gn[i]) == 0)
2660 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2664 } /* search_QMstring */
2666 /* We would like gn to be const as well, but C doesn't allow this */
2667 /* TODO this is utility functionality (search for the index of a
2668 string in a collection), so should be refactored and located more
2670 int search_string(const char *s, int ng, char *gn[])
2674 for (i = 0; (i < ng); i++)
2676 if (gmx_strcasecmp(s, gn[i]) == 0)
2683 "Group %s referenced in the .mdp file was not found in the index file.\n"
2684 "Group names must match either [moleculetype] names or custom index group\n"
2685 "names, in which case you must supply an index file to the '-n' option\n"
2692 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2693 t_blocka *block, char *gnames[],
2694 int gtype, int restnm,
2695 int grptp, gmx_bool bVerbose,
2698 unsigned short *cbuf;
2699 t_grps *grps = &(groups->grps[gtype]);
2700 int i, j, gid, aj, ognr, ntot = 0;
2703 char warn_buf[STRLEN];
2707 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2710 title = gtypes[gtype];
2713 /* Mark all id's as not set */
2714 for (i = 0; (i < natoms); i++)
2719 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2720 for (i = 0; (i < ng); i++)
2722 /* Lookup the group name in the block structure */
2723 gid = search_string(ptrs[i], block->nr, gnames);
2724 if ((grptp != egrptpONE) || (i == 0))
2726 grps->nm_ind[grps->nr++] = gid;
2730 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2733 /* Now go over the atoms in the group */
2734 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2739 /* Range checking */
2740 if ((aj < 0) || (aj >= natoms))
2742 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2744 /* Lookup up the old group number */
2748 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2749 aj+1, title, ognr+1, i+1);
2753 /* Store the group number in buffer */
2754 if (grptp == egrptpONE)
2767 /* Now check whether we have done all atoms */
2771 if (grptp == egrptpALL)
2773 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2774 natoms-ntot, title);
2776 else if (grptp == egrptpPART)
2778 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2779 natoms-ntot, title);
2780 warning_note(wi, warn_buf);
2782 /* Assign all atoms currently unassigned to a rest group */
2783 for (j = 0; (j < natoms); j++)
2785 if (cbuf[j] == NOGID)
2791 if (grptp != egrptpPART)
2796 "Making dummy/rest group for %s containing %d elements\n",
2797 title, natoms-ntot);
2799 /* Add group name "rest" */
2800 grps->nm_ind[grps->nr] = restnm;
2802 /* Assign the rest name to all atoms not currently assigned to a group */
2803 for (j = 0; (j < natoms); j++)
2805 if (cbuf[j] == NOGID)
2814 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2816 /* All atoms are part of one (or no) group, no index required */
2817 groups->ngrpnr[gtype] = 0;
2818 groups->grpnr[gtype] = nullptr;
2822 groups->ngrpnr[gtype] = natoms;
2823 snew(groups->grpnr[gtype], natoms);
2824 for (j = 0; (j < natoms); j++)
2826 groups->grpnr[gtype][j] = cbuf[j];
2832 return (bRest && grptp == egrptpPART);
2835 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2838 gmx_groups_t *groups;
2839 pull_params_t *pull;
2840 int natoms, ai, aj, i, j, d, g, imin, jmin;
2842 int *nrdf2, *na_vcm, na_tot;
2843 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2845 gmx_mtop_atomloop_all_t aloop;
2846 int mb, mol, ftype, as;
2847 gmx_molblock_t *molb;
2848 gmx_moltype_t *molt;
2851 * First calc 3xnr-atoms for each group
2852 * then subtract half a degree of freedom for each constraint
2854 * Only atoms and nuclei contribute to the degrees of freedom...
2859 groups = &mtop->groups;
2860 natoms = mtop->natoms;
2862 /* Allocate one more for a possible rest group */
2863 /* We need to sum degrees of freedom into doubles,
2864 * since floats give too low nrdf's above 3 million atoms.
2866 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2867 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2868 snew(dof_vcm, groups->grps[egcVCM].nr+1);
2869 snew(na_vcm, groups->grps[egcVCM].nr+1);
2870 snew(nrdf_vcm_sub, groups->grps[egcVCM].nr+1);
2872 for (i = 0; i < groups->grps[egcTC].nr; i++)
2876 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2879 clear_ivec(dof_vcm[i]);
2881 nrdf_vcm_sub[i] = 0;
2884 snew(nrdf2, natoms);
2885 aloop = gmx_mtop_atomloop_all_init(mtop);
2887 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2890 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2892 g = ggrpnr(groups, egcFREEZE, i);
2893 for (d = 0; d < DIM; d++)
2895 if (opts->nFreeze[g][d] == 0)
2897 /* Add one DOF for particle i (counted as 2*1) */
2899 /* VCM group i has dim d as a DOF */
2900 dof_vcm[ggrpnr(groups, egcVCM, i)][d] = 1;
2903 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2904 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2909 for (mb = 0; mb < mtop->nmolblock; mb++)
2911 molb = &mtop->molblock[mb];
2912 molt = &mtop->moltype[molb->type];
2913 atom = molt->atoms.atom;
2914 for (mol = 0; mol < molb->nmol; mol++)
2916 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2918 ia = molt->ilist[ftype].iatoms;
2919 for (i = 0; i < molt->ilist[ftype].nr; )
2921 /* Subtract degrees of freedom for the constraints,
2922 * if the particles still have degrees of freedom left.
2923 * If one of the particles is a vsite or a shell, then all
2924 * constraint motion will go there, but since they do not
2925 * contribute to the constraints the degrees of freedom do not
2930 if (((atom[ia[1]].ptype == eptNucleus) ||
2931 (atom[ia[1]].ptype == eptAtom)) &&
2932 ((atom[ia[2]].ptype == eptNucleus) ||
2933 (atom[ia[2]].ptype == eptAtom)))
2951 imin = std::min(imin, nrdf2[ai]);
2952 jmin = std::min(jmin, nrdf2[aj]);
2955 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2956 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2957 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2958 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2960 ia += interaction_function[ftype].nratoms+1;
2961 i += interaction_function[ftype].nratoms+1;
2964 ia = molt->ilist[F_SETTLE].iatoms;
2965 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2967 /* Subtract 1 dof from every atom in the SETTLE */
2968 for (j = 0; j < 3; j++)
2971 imin = std::min(2, nrdf2[ai]);
2973 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2974 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2979 as += molt->atoms.nr;
2985 /* Correct nrdf for the COM constraints.
2986 * We correct using the TC and VCM group of the first atom
2987 * in the reference and pull group. If atoms in one pull group
2988 * belong to different TC or VCM groups it is anyhow difficult
2989 * to determine the optimal nrdf assignment.
2993 for (i = 0; i < pull->ncoord; i++)
2995 if (pull->coord[i].eType != epullCONSTRAINT)
3002 for (j = 0; j < 2; j++)
3004 const t_pull_group *pgrp;
3006 pgrp = &pull->group[pull->coord[i].group[j]];
3010 /* Subtract 1/2 dof from each group */
3012 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
3013 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
3014 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
3016 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)]]);
3021 /* We need to subtract the whole DOF from group j=1 */
3028 if (ir->nstcomm != 0)
3032 /* We remove COM motion up to dim ndof_com() */
3033 ndim_rm_vcm = ndof_com(ir);
3035 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3036 * the number of degrees of freedom in each vcm group when COM
3037 * translation is removed and 6 when rotation is removed as well.
3039 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3041 switch (ir->comm_mode)
3044 case ecmLINEAR_ACCELERATION_CORRECTION:
3045 nrdf_vcm_sub[j] = 0;
3046 for (d = 0; d < ndim_rm_vcm; d++)
3055 nrdf_vcm_sub[j] = 6;
3058 gmx_incons("Checking comm_mode");
3062 for (i = 0; i < groups->grps[egcTC].nr; i++)
3064 /* Count the number of atoms of TC group i for every VCM group */
3065 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3070 for (ai = 0; ai < natoms; ai++)
3072 if (ggrpnr(groups, egcTC, ai) == i)
3074 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
3078 /* Correct for VCM removal according to the fraction of each VCM
3079 * group present in this TC group.
3081 nrdf_uc = nrdf_tc[i];
3084 fprintf(debug, "T-group[%d] nrdf_uc = %g\n", i, nrdf_uc);
3087 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3089 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3091 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
3092 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
3096 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
3097 j, nrdf_vcm[j], nrdf_tc[i]);
3102 for (i = 0; (i < groups->grps[egcTC].nr); i++)
3104 opts->nrdf[i] = nrdf_tc[i];
3105 if (opts->nrdf[i] < 0)
3110 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3111 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
3119 sfree(nrdf_vcm_sub);
3122 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
3123 const char *option, const char *val, int flag)
3125 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3126 * But since this is much larger than STRLEN, such a line can not be parsed.
3127 * The real maximum is the number of names that fit in a string: STRLEN/2.
3129 #define EGP_MAX (STRLEN/2)
3130 int nelem, i, j, k, nr;
3131 char *names[EGP_MAX];
3135 gnames = groups->grpname;
3137 nelem = str_nelem(val, EGP_MAX, names);
3140 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3142 nr = groups->grps[egcENER].nr;
3144 for (i = 0; i < nelem/2; i++)
3148 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
3154 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3155 names[2*i], option);
3159 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3165 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3166 names[2*i+1], option);
3168 if ((j < nr) && (k < nr))
3170 ir->opts.egp_flags[nr*j+k] |= flag;
3171 ir->opts.egp_flags[nr*k+j] |= flag;
3180 static void make_swap_groups(
3185 int ig = -1, i = 0, gind;
3189 /* Just a quick check here, more thorough checks are in mdrun */
3190 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3192 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3195 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3196 for (ig = 0; ig < swap->ngrp; ig++)
3198 swapg = &swap->grp[ig];
3199 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3200 swapg->nat = grps->index[gind+1] - grps->index[gind];
3204 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3205 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3206 swap->grp[ig].molname, swapg->nat);
3207 snew(swapg->ind, swapg->nat);
3208 for (i = 0; i < swapg->nat; i++)
3210 swapg->ind[i] = grps->a[grps->index[gind]+i];
3215 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3221 static void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3226 ig = search_string(IMDgname, grps->nr, gnames);
3227 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3229 if (IMDgroup->nat > 0)
3231 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3232 IMDgname, IMDgroup->nat);
3233 snew(IMDgroup->ind, IMDgroup->nat);
3234 for (i = 0; i < IMDgroup->nat; i++)
3236 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3242 void do_index(const char* mdparin, const char *ndx,
3249 gmx_groups_t *groups;
3253 char warnbuf[STRLEN], **gnames;
3254 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3257 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3258 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3259 int i, j, k, restnm;
3260 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3261 int nQMmethod, nQMbasis, nQMg;
3262 char warn_buf[STRLEN];
3267 fprintf(stderr, "processing index file...\n");
3272 snew(grps->index, 1);
3274 atoms_all = gmx_mtop_global_atoms(mtop);
3275 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3276 done_atom(&atoms_all);
3280 grps = init_index(ndx, &gnames);
3283 groups = &mtop->groups;
3284 natoms = mtop->natoms;
3285 symtab = &mtop->symtab;
3287 snew(groups->grpname, grps->nr+1);
3289 for (i = 0; (i < grps->nr); i++)
3291 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3293 groups->grpname[i] = put_symtab(symtab, "rest");
3295 srenew(gnames, grps->nr+1);
3296 gnames[restnm] = *(groups->grpname[i]);
3297 groups->ngrpname = grps->nr+1;
3299 set_warning_line(wi, mdparin, -1);
3301 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3302 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3303 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3304 if ((ntau_t != ntcg) || (nref_t != ntcg))
3306 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3307 "%d tau-t values", ntcg, nref_t, ntau_t);
3310 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3311 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3312 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3313 nr = groups->grps[egcTC].nr;
3315 snew(ir->opts.nrdf, nr);
3316 snew(ir->opts.tau_t, nr);
3317 snew(ir->opts.ref_t, nr);
3318 if (ir->eI == eiBD && ir->bd_fric == 0)
3320 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3327 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3331 for (i = 0; (i < nr); i++)
3333 ir->opts.tau_t[i] = strtod(ptr1[i], &endptr);
3336 warning_error(wi, "Invalid value for mdp option tau-t. tau-t should only consist of real numbers separated by spaces.");
3338 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3340 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3341 warning_error(wi, warn_buf);
3344 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3346 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.");
3349 if (ir->opts.tau_t[i] >= 0)
3351 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3354 if (ir->etc != etcNO && ir->nsttcouple == -1)
3356 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3361 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3363 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");
3365 if (ir->epc == epcMTTK)
3367 if (ir->etc != etcNOSEHOOVER)
3369 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3373 if (ir->nstpcouple != ir->nsttcouple)
3375 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3376 ir->nstpcouple = ir->nsttcouple = mincouple;
3377 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);
3378 warning_note(wi, warn_buf);
3383 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3384 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3386 if (ETC_ANDERSEN(ir->etc))
3388 if (ir->nsttcouple != 1)
3391 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");
3392 warning_note(wi, warn_buf);
3395 nstcmin = tcouple_min_integration_steps(ir->etc);
3398 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3400 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3401 ETCOUPLTYPE(ir->etc),
3403 ir->nsttcouple*ir->delta_t);
3404 warning(wi, warn_buf);
3407 for (i = 0; (i < nr); i++)
3409 ir->opts.ref_t[i] = strtod(ptr2[i], &endptr);
3412 warning_error(wi, "Invalid value for mdp option ref-t. ref-t should only consist of real numbers separated by spaces.");
3414 if (ir->opts.ref_t[i] < 0)
3416 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3419 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3420 if we are in this conditional) if mc_temp is negative */
3421 if (ir->expandedvals->mc_temp < 0)
3423 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3427 /* Simulated annealing for each group. There are nr groups */
3428 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3429 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3433 if (nSA > 0 && nSA != nr)
3435 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3439 snew(ir->opts.annealing, nr);
3440 snew(ir->opts.anneal_npoints, nr);
3441 snew(ir->opts.anneal_time, nr);
3442 snew(ir->opts.anneal_temp, nr);
3443 for (i = 0; i < nr; i++)
3445 ir->opts.annealing[i] = eannNO;
3446 ir->opts.anneal_npoints[i] = 0;
3447 ir->opts.anneal_time[i] = nullptr;
3448 ir->opts.anneal_temp[i] = nullptr;
3453 for (i = 0; i < nr; i++)
3455 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3457 ir->opts.annealing[i] = eannNO;
3459 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3461 ir->opts.annealing[i] = eannSINGLE;
3464 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3466 ir->opts.annealing[i] = eannPERIODIC;
3472 /* Read the other fields too */
3473 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3474 if (nSA_points != nSA)
3476 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3478 for (k = 0, i = 0; i < nr; i++)
3480 ir->opts.anneal_npoints[i] = strtol(ptr1[i], &endptr, 10);
3483 warning_error(wi, "Invalid value for mdp option annealing-npoints. annealing should only consist of integers separated by spaces.");
3485 if (ir->opts.anneal_npoints[i] == 1)
3487 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3489 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3490 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3491 k += ir->opts.anneal_npoints[i];
3494 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3497 gmx_fatal(FARGS, "Found %d annealing-time values, wanted %d\n", nSA_time, k);
3499 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3502 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3505 for (i = 0, k = 0; i < nr; i++)
3508 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3510 ir->opts.anneal_time[i][j] = strtod(ptr1[k], &endptr);
3513 warning_error(wi, "Invalid value for mdp option anneal-time. anneal-time should only consist of real numbers separated by spaces.");
3515 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], &endptr);
3518 warning_error(wi, "Invalid value for anneal-temp. anneal-temp should only consist of real numbers separated by spaces.");
3522 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3524 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3530 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3532 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3533 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3536 if (ir->opts.anneal_temp[i][j] < 0)
3538 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3543 /* Print out some summary information, to make sure we got it right */
3544 for (i = 0, k = 0; i < nr; i++)
3546 if (ir->opts.annealing[i] != eannNO)
3548 j = groups->grps[egcTC].nm_ind[i];
3549 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3550 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3551 ir->opts.anneal_npoints[i]);
3552 fprintf(stderr, "Time (ps) Temperature (K)\n");
3553 /* All terms except the last one */
3554 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3556 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3559 /* Finally the last one */
3560 j = ir->opts.anneal_npoints[i]-1;
3561 if (ir->opts.annealing[i] == eannSINGLE)
3563 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3567 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3568 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3570 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3581 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3583 make_pull_coords(ir->pull);
3588 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3591 if (ir->eSwapCoords != eswapNO)
3593 make_swap_groups(ir->swap, grps, gnames);
3596 /* Make indices for IMD session */
3599 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3602 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3603 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3604 if (nacg*DIM != nacc)
3606 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3609 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3610 restnm, egrptpALL_GENREST, bVerbose, wi);
3611 nr = groups->grps[egcACC].nr;
3612 snew(ir->opts.acc, nr);
3613 ir->opts.ngacc = nr;
3615 for (i = k = 0; (i < nacg); i++)
3617 for (j = 0; (j < DIM); j++, k++)
3619 ir->opts.acc[i][j] = strtod(ptr1[k], &endptr);
3622 warning_error(wi, "Invalid value for mdp option accelerate. accelerate should only consist of real numbers separated by spaces.");
3626 for (; (i < nr); i++)
3628 for (j = 0; (j < DIM); j++)
3630 ir->opts.acc[i][j] = 0;
3634 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3635 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3636 if (nfrdim != DIM*nfreeze)
3638 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3641 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3642 restnm, egrptpALL_GENREST, bVerbose, wi);
3643 nr = groups->grps[egcFREEZE].nr;
3644 ir->opts.ngfrz = nr;
3645 snew(ir->opts.nFreeze, nr);
3646 for (i = k = 0; (i < nfreeze); i++)
3648 for (j = 0; (j < DIM); j++, k++)
3650 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3651 if (!ir->opts.nFreeze[i][j])
3653 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3655 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3656 "(not %s)", ptr1[k]);
3657 warning(wi, warn_buf);
3662 for (; (i < nr); i++)
3664 for (j = 0; (j < DIM); j++)
3666 ir->opts.nFreeze[i][j] = 0;
3670 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3671 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3672 restnm, egrptpALL_GENREST, bVerbose, wi);
3673 add_wall_energrps(groups, ir->nwall, symtab);
3674 ir->opts.ngener = groups->grps[egcENER].nr;
3675 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3677 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3678 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3681 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3682 "This may lead to artifacts.\n"
3683 "In most cases one should use one group for the whole system.");
3686 /* Now we have filled the freeze struct, so we can calculate NRDF */
3687 calc_nrdf(mtop, ir, gnames);
3689 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3690 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3691 restnm, egrptpALL_GENREST, bVerbose, wi);
3692 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3693 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3694 restnm, egrptpALL_GENREST, bVerbose, wi);
3695 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3696 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3697 restnm, egrptpONE, bVerbose, wi);
3698 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3699 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3700 restnm, egrptpALL_GENREST, bVerbose, wi);
3702 /* QMMM input processing */
3703 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3704 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3705 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3706 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3708 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3709 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3711 /* group rest, if any, is always MM! */
3712 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3713 restnm, egrptpALL_GENREST, bVerbose, wi);
3714 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3715 ir->opts.ngQM = nQMg;
3716 snew(ir->opts.QMmethod, nr);
3717 snew(ir->opts.QMbasis, nr);
3718 for (i = 0; i < nr; i++)
3720 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3721 * converted to the corresponding enum in names.c
3723 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3725 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3729 str_nelem(is->QMmult, MAXPTR, ptr1);
3730 str_nelem(is->QMcharge, MAXPTR, ptr2);
3731 str_nelem(is->bSH, MAXPTR, ptr3);
3732 snew(ir->opts.QMmult, nr);
3733 snew(ir->opts.QMcharge, nr);
3734 snew(ir->opts.bSH, nr);
3736 for (i = 0; i < nr; i++)
3738 ir->opts.QMmult[i] = strtol(ptr1[i], &endptr, 10);
3741 warning_error(wi, "Invalid value for mdp option QMmult. QMmult should only consist of integers separated by spaces.");
3743 ir->opts.QMcharge[i] = strtol(ptr2[i], &endptr, 10);
3746 warning_error(wi, "Invalid value for mdp option QMcharge. QMcharge should only consist of integers separated by spaces.");
3748 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3751 str_nelem(is->CASelectrons, MAXPTR, ptr1);
3752 str_nelem(is->CASorbitals, MAXPTR, ptr2);
3753 snew(ir->opts.CASelectrons, nr);
3754 snew(ir->opts.CASorbitals, nr);
3755 for (i = 0; i < nr; i++)
3757 ir->opts.CASelectrons[i] = strtol(ptr1[i], &endptr, 10);
3760 warning_error(wi, "Invalid value for mdp option CASelectrons. CASelectrons should only consist of integers separated by spaces.");
3762 ir->opts.CASorbitals[i] = strtol(ptr2[i], &endptr, 10);
3765 warning_error(wi, "Invalid value for mdp option CASorbitals. CASorbitals should only consist of integers separated by spaces.");
3769 str_nelem(is->SAon, MAXPTR, ptr1);
3770 str_nelem(is->SAoff, MAXPTR, ptr2);
3771 str_nelem(is->SAsteps, MAXPTR, ptr3);
3772 snew(ir->opts.SAon, nr);
3773 snew(ir->opts.SAoff, nr);
3774 snew(ir->opts.SAsteps, nr);
3776 for (i = 0; i < nr; i++)
3778 ir->opts.SAon[i] = strtod(ptr1[i], &endptr);
3781 warning_error(wi, "Invalid value for mdp option SAon. SAon should only consist of real numbers separated by spaces.");
3783 ir->opts.SAoff[i] = strtod(ptr2[i], &endptr);
3786 warning_error(wi, "Invalid value for mdp option SAoff. SAoff should only consist of real numbers separated by spaces.");
3788 ir->opts.SAsteps[i] = strtol(ptr3[i], &endptr, 10);
3791 warning_error(wi, "Invalid value for mdp option SAsteps. SAsteps should only consist of integers separated by spaces.");
3794 /* end of QMMM input */
3798 for (i = 0; (i < egcNR); i++)
3800 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3801 for (j = 0; (j < groups->grps[i].nr); j++)
3803 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3805 fprintf(stderr, "\n");
3809 nr = groups->grps[egcENER].nr;
3810 snew(ir->opts.egp_flags, nr*nr);
3812 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3813 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3815 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3817 if (bExcl && EEL_FULL(ir->coulombtype))
3819 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3822 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3823 if (bTable && !(ir->vdwtype == evdwUSER) &&
3824 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3825 !(ir->coulombtype == eelPMEUSERSWITCH))
3827 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3830 for (i = 0; (i < grps->nr); i++)
3842 static void check_disre(gmx_mtop_t *mtop)
3844 gmx_ffparams_t *ffparams;
3845 t_functype *functype;
3847 int i, ndouble, ftype;
3848 int label, old_label;
3850 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3852 ffparams = &mtop->ffparams;
3853 functype = ffparams->functype;
3854 ip = ffparams->iparams;
3857 for (i = 0; i < ffparams->ntypes; i++)
3859 ftype = functype[i];
3860 if (ftype == F_DISRES)
3862 label = ip[i].disres.label;
3863 if (label == old_label)
3865 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3873 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3874 "probably the parameters for multiple pairs in one restraint "
3875 "are not identical\n", ndouble);
3880 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3881 gmx_bool posres_only,
3885 gmx_mtop_ilistloop_t iloop;
3895 for (d = 0; d < DIM; d++)
3897 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3899 /* Check for freeze groups */
3900 for (g = 0; g < ir->opts.ngfrz; g++)
3902 for (d = 0; d < DIM; d++)
3904 if (ir->opts.nFreeze[g][d] != 0)
3912 /* Check for position restraints */
3913 iloop = gmx_mtop_ilistloop_init(sys);
3914 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3917 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3919 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3921 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3922 for (d = 0; d < DIM; d++)
3924 if (pr->posres.fcA[d] != 0)
3930 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3932 /* Check for flat-bottom posres */
3933 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3934 if (pr->fbposres.k != 0)
3936 switch (pr->fbposres.geom)
3938 case efbposresSPHERE:
3939 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3941 case efbposresCYLINDERX:
3942 AbsRef[YY] = AbsRef[ZZ] = 1;
3944 case efbposresCYLINDERY:
3945 AbsRef[XX] = AbsRef[ZZ] = 1;
3947 case efbposresCYLINDER:
3948 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3949 case efbposresCYLINDERZ:
3950 AbsRef[XX] = AbsRef[YY] = 1;
3952 case efbposresX: /* d=XX */
3953 case efbposresY: /* d=YY */
3954 case efbposresZ: /* d=ZZ */
3955 d = pr->fbposres.geom - efbposresX;
3959 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3960 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3968 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3972 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3973 gmx_bool *bC6ParametersWorkWithGeometricRules,
3974 gmx_bool *bC6ParametersWorkWithLBRules,
3975 gmx_bool *bLBRulesPossible)
3977 int ntypes, tpi, tpj;
3980 double c6i, c6j, c12i, c12j;
3981 double c6, c6_geometric, c6_LB;
3982 double sigmai, sigmaj, epsi, epsj;
3983 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3986 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3987 * force-field floating point parameters.
3990 ptr = getenv("GMX_LJCOMB_TOL");
3994 double gmx_unused canary;
3996 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3998 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4003 *bC6ParametersWorkWithLBRules = TRUE;
4004 *bC6ParametersWorkWithGeometricRules = TRUE;
4005 bCanDoLBRules = TRUE;
4006 ntypes = mtop->ffparams.atnr;
4007 snew(typecount, ntypes);
4008 gmx_mtop_count_atomtypes(mtop, state, typecount);
4009 *bLBRulesPossible = TRUE;
4010 for (tpi = 0; tpi < ntypes; ++tpi)
4012 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4013 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4014 for (tpj = tpi; tpj < ntypes; ++tpj)
4016 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4017 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4018 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4019 c6_geometric = std::sqrt(c6i * c6j);
4020 if (!gmx_numzero(c6_geometric))
4022 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4024 sigmai = gmx::sixthroot(c12i / c6i);
4025 sigmaj = gmx::sixthroot(c12j / c6j);
4026 epsi = c6i * c6i /(4.0 * c12i);
4027 epsj = c6j * c6j /(4.0 * c12j);
4028 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4032 *bLBRulesPossible = FALSE;
4033 c6_LB = c6_geometric;
4035 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4038 if (FALSE == bCanDoLBRules)
4040 *bC6ParametersWorkWithLBRules = FALSE;
4043 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4045 if (FALSE == bCanDoGeometricRules)
4047 *bC6ParametersWorkWithGeometricRules = FALSE;
4055 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
4058 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4060 check_combination_rule_differences(mtop, 0,
4061 &bC6ParametersWorkWithGeometricRules,
4062 &bC6ParametersWorkWithLBRules,
4064 if (ir->ljpme_combination_rule == eljpmeLB)
4066 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
4068 warning(wi, "You are using arithmetic-geometric combination rules "
4069 "in LJ-PME, but your non-bonded C6 parameters do not "
4070 "follow these rules.");
4075 if (FALSE == bC6ParametersWorkWithGeometricRules)
4077 if (ir->eDispCorr != edispcNO)
4079 warning_note(wi, "You are using geometric combination rules in "
4080 "LJ-PME, but your non-bonded C6 parameters do "
4081 "not follow these rules. "
4082 "This will introduce very small errors in the forces and energies in "
4083 "your simulations. Dispersion correction will correct total energy "
4084 "and/or pressure for isotropic systems, but not forces or surface tensions.");
4088 warning_note(wi, "You are using geometric combination rules in "
4089 "LJ-PME, but your non-bonded C6 parameters do "
4090 "not follow these rules. "
4091 "This will introduce very small errors in the forces and energies in "
4092 "your simulations. If your system is homogeneous, consider using dispersion correction "
4093 "for the total energy and pressure.");
4099 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
4102 char err_buf[STRLEN];
4104 gmx_bool bCharge, bAcc;
4107 gmx_mtop_atomloop_block_t aloopb;
4108 gmx_mtop_atomloop_all_t aloop;
4110 char warn_buf[STRLEN];
4112 set_warning_line(wi, mdparin, -1);
4114 if (ir->cutoff_scheme == ecutsVERLET &&
4115 ir->verletbuf_tol > 0 &&
4117 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
4118 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4120 /* Check if a too small Verlet buffer might potentially
4121 * cause more drift than the thermostat can couple off.
4123 /* Temperature error fraction for warning and suggestion */
4124 const real T_error_warn = 0.002;
4125 const real T_error_suggest = 0.001;
4126 /* For safety: 2 DOF per atom (typical with constraints) */
4127 const real nrdf_at = 2;
4128 real T, tau, max_T_error;
4133 for (i = 0; i < ir->opts.ngtc; i++)
4135 T = std::max(T, ir->opts.ref_t[i]);
4136 tau = std::max(tau, ir->opts.tau_t[i]);
4140 /* This is a worst case estimate of the temperature error,
4141 * assuming perfect buffer estimation and no cancelation
4142 * of errors. The factor 0.5 is because energy distributes
4143 * equally over Ekin and Epot.
4145 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4146 if (max_T_error > T_error_warn)
4148 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.",
4149 ir->verletbuf_tol, T, tau,
4151 100*T_error_suggest,
4152 ir->verletbuf_tol*T_error_suggest/max_T_error);
4153 warning(wi, warn_buf);
4158 if (ETC_ANDERSEN(ir->etc))
4162 for (i = 0; i < ir->opts.ngtc; i++)
4164 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4165 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4166 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
4167 i, ir->opts.tau_t[i]);
4168 CHECK(ir->opts.tau_t[i] < 0);
4171 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
4173 for (i = 0; i < ir->opts.ngtc; i++)
4175 int nsteps = static_cast<int>(ir->opts.tau_t[i]/ir->delta_t + 0.5);
4176 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);
4177 CHECK(nsteps % ir->nstcomm != 0);
4182 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4183 ir->comm_mode == ecmNO &&
4184 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4185 !ETC_ANDERSEN(ir->etc))
4187 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");
4190 /* Check for pressure coupling with absolute position restraints */
4191 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4193 absolute_reference(ir, sys, TRUE, AbsRef);
4195 for (m = 0; m < DIM; m++)
4197 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4199 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4207 aloopb = gmx_mtop_atomloop_block_init(sys);
4209 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4211 if (atom->q != 0 || atom->qB != 0)
4219 if (EEL_FULL(ir->coulombtype))
4222 "You are using full electrostatics treatment %s for a system without charges.\n"
4223 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4224 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4225 warning(wi, err_buf);
4230 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4233 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4234 "You might want to consider using %s electrostatics.\n",
4236 warning_note(wi, err_buf);
4240 /* Check if combination rules used in LJ-PME are the same as in the force field */
4241 if (EVDW_PME(ir->vdwtype))
4243 check_combination_rules(ir, sys, wi);
4246 /* Generalized reaction field */
4247 if (ir->opts.ngtc == 0)
4249 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4251 CHECK(ir->coulombtype == eelGRF);
4255 sprintf(err_buf, "When using coulombtype = %s"
4256 " ref-t for temperature coupling should be > 0",
4258 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4262 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4264 for (m = 0; (m < DIM); m++)
4266 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4275 snew(mgrp, sys->groups.grps[egcACC].nr);
4276 aloop = gmx_mtop_atomloop_all_init(sys);
4278 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4280 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4283 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4285 for (m = 0; (m < DIM); m++)
4287 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4291 for (m = 0; (m < DIM); m++)
4293 if (fabs(acc[m]) > 1e-6)
4295 const char *dim[DIM] = { "X", "Y", "Z" };
4297 "Net Acceleration in %s direction, will %s be corrected\n",
4298 dim[m], ir->nstcomm != 0 ? "" : "not");
4299 if (ir->nstcomm != 0 && m < ndof_com(ir))
4302 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4304 ir->opts.acc[i][m] -= acc[m];
4312 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4313 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4315 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4323 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4325 if (ir->pull->coord[i].group[0] == 0 ||
4326 ir->pull->coord[i].group[1] == 0)
4328 absolute_reference(ir, sys, FALSE, AbsRef);
4329 for (m = 0; m < DIM; m++)
4331 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4333 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.");
4341 for (i = 0; i < 3; i++)
4343 for (m = 0; m <= i; m++)
4345 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4346 ir->deform[i][m] != 0)
4348 for (c = 0; c < ir->pull->ncoord; c++)
4350 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4351 ir->pull->coord[c].vec[m] != 0)
4353 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4364 void double_check(t_inputrec *ir, matrix box,
4365 gmx_bool bHasNormalConstraints,
4366 gmx_bool bHasAnyConstraints,
4370 char warn_buf[STRLEN];
4373 ptr = check_box(ir->ePBC, box);
4376 warning_error(wi, ptr);
4379 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4381 if (ir->shake_tol <= 0.0)
4383 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4385 warning_error(wi, warn_buf);
4389 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4391 /* If we have Lincs constraints: */
4392 if (ir->eI == eiMD && ir->etc == etcNO &&
4393 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4395 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4396 warning_note(wi, warn_buf);
4399 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4401 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4402 warning_note(wi, warn_buf);
4404 if (ir->epc == epcMTTK)
4406 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4410 if (bHasAnyConstraints && ir->epc == epcMTTK)
4412 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4415 if (ir->LincsWarnAngle > 90.0)
4417 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4418 warning(wi, warn_buf);
4419 ir->LincsWarnAngle = 90.0;
4422 if (ir->ePBC != epbcNONE)
4424 if (ir->nstlist == 0)
4426 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4428 if (ir->ns_type == ensGRID)
4430 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4432 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");
4433 warning_error(wi, warn_buf);
4438 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4439 if (2*ir->rlist >= min_size)
4441 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4442 warning_error(wi, warn_buf);
4445 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4452 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4456 real rvdw1, rvdw2, rcoul1, rcoul2;
4457 char warn_buf[STRLEN];
4459 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4463 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4468 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4474 if (rvdw1 + rvdw2 > ir->rlist ||
4475 rcoul1 + rcoul2 > ir->rlist)
4478 "The sum of the two largest charge group radii (%f) "
4479 "is larger than rlist (%f)\n",
4480 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4481 warning(wi, warn_buf);
4485 /* Here we do not use the zero at cut-off macro,
4486 * since user defined interactions might purposely
4487 * not be zero at the cut-off.
4489 if (ir_vdw_is_zero_at_cutoff(ir) &&
4490 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4492 sprintf(warn_buf, "The sum of the two largest charge group "
4493 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4494 "With exact cut-offs, better performance can be "
4495 "obtained with cutoff-scheme = %s, because it "
4496 "does not use charge groups at all.",
4498 ir->rlist, ir->rvdw,
4499 ecutscheme_names[ecutsVERLET]);
4502 warning(wi, warn_buf);
4506 warning_note(wi, warn_buf);
4509 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4510 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4512 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4513 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4515 ir->rlist, ir->rcoulomb,
4516 ecutscheme_names[ecutsVERLET]);
4519 warning(wi, warn_buf);
4523 warning_note(wi, warn_buf);