2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018, 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) && !EI_SD(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
910 sprintf(warn_buf, "Tumbling and flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR or use integrator = %s.", ei_names[eiSD1]);
911 warning_note(wi, warn_buf);
914 /* TEMPERATURE COUPLING */
915 if (ir->etc == etcYES)
917 ir->etc = etcBERENDSEN;
918 warning_note(wi, "Old option for temperature coupling given: "
919 "changing \"yes\" to \"Berendsen\"\n");
922 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
924 if (ir->opts.nhchainlength < 1)
926 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
927 ir->opts.nhchainlength = 1;
928 warning(wi, warn_buf);
931 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
933 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
934 ir->opts.nhchainlength = 1;
939 ir->opts.nhchainlength = 0;
942 if (ir->eI == eiVVAK)
944 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
946 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
949 if (ETC_ANDERSEN(ir->etc))
951 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
952 CHECK(!(EI_VV(ir->eI)));
954 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
956 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]);
957 warning_note(wi, warn_buf);
960 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]);
961 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
964 if (ir->etc == etcBERENDSEN)
966 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
967 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
968 warning_note(wi, warn_buf);
971 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
972 && ir->epc == epcBERENDSEN)
974 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
975 "true ensemble for the thermostat");
976 warning(wi, warn_buf);
979 /* PRESSURE COUPLING */
980 if (ir->epc == epcISOTROPIC)
982 ir->epc = epcBERENDSEN;
983 warning_note(wi, "Old option for pressure coupling given: "
984 "changing \"Isotropic\" to \"Berendsen\"\n");
987 if (ir->epc != epcNO)
989 dt_pcoupl = ir->nstpcouple*ir->delta_t;
991 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
992 CHECK(ir->tau_p <= 0);
994 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
996 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
997 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
998 warning(wi, warn_buf);
1001 sprintf(err_buf, "compressibility must be > 0 when using pressure"
1002 " coupling %s\n", EPCOUPLTYPE(ir->epc));
1003 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
1004 ir->compress[ZZ][ZZ] < 0 ||
1005 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
1006 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1008 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1011 "You are generating velocities so I am assuming you "
1012 "are equilibrating a system. You are using "
1013 "%s pressure coupling, but this can be "
1014 "unstable for equilibration. If your system crashes, try "
1015 "equilibrating first with Berendsen pressure coupling. If "
1016 "you are not equilibrating the system, you can probably "
1017 "ignore this warning.",
1018 epcoupl_names[ir->epc]);
1019 warning(wi, warn_buf);
1025 if (ir->epc > epcNO)
1027 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1029 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.");
1035 if (ir->epc == epcMTTK)
1037 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1041 /* ELECTROSTATICS */
1042 /* More checks are in triple check (grompp.c) */
1044 if (ir->coulombtype == eelSWITCH)
1046 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1047 "artifacts, advice: use coulombtype = %s",
1048 eel_names[ir->coulombtype],
1049 eel_names[eelRF_ZERO]);
1050 warning(wi, warn_buf);
1053 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1055 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1056 warning_note(wi, warn_buf);
1059 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1061 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);
1062 warning(wi, warn_buf);
1063 ir->epsilon_rf = ir->epsilon_r;
1064 ir->epsilon_r = 1.0;
1067 if (ir->epsilon_r == 0)
1070 "It is pointless to use long-range or Generalized Born electrostatics with infinite relative permittivity."
1071 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1072 CHECK(EEL_FULL(ir->coulombtype) || ir->implicit_solvent == eisGBSA);
1075 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1077 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1078 CHECK(ir->epsilon_r < 0);
1081 if (EEL_RF(ir->coulombtype))
1083 /* reaction field (at the cut-off) */
1085 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1087 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1088 eel_names[ir->coulombtype]);
1089 warning(wi, warn_buf);
1090 ir->epsilon_rf = 0.0;
1093 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1094 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1095 (ir->epsilon_r == 0));
1096 if (ir->epsilon_rf == ir->epsilon_r)
1098 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1099 eel_names[ir->coulombtype]);
1100 warning(wi, warn_buf);
1103 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1104 * means the interaction is zero outside rcoulomb, but it helps to
1105 * provide accurate energy conservation.
1107 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1109 if (ir_coulomb_switched(ir))
1112 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1113 eel_names[ir->coulombtype]);
1114 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1117 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1119 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1121 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1122 eel_names[ir->coulombtype]);
1123 CHECK(ir->rlist > ir->rcoulomb);
1127 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1130 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1131 CHECK( ir->coulomb_modifier != eintmodNONE);
1133 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1136 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1137 CHECK( ir->vdw_modifier != eintmodNONE);
1140 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1141 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1144 "The switch/shift interaction settings are just for compatibility; you will get better "
1145 "performance from applying potential modifiers to your interactions!\n");
1146 warning_note(wi, warn_buf);
1149 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1151 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1153 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1154 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.",
1155 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1156 warning(wi, warn_buf);
1160 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1162 if (ir->rvdw_switch == 0)
1164 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.");
1165 warning(wi, warn_buf);
1169 if (EEL_FULL(ir->coulombtype))
1171 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1172 ir->coulombtype == eelPMEUSERSWITCH)
1174 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1175 eel_names[ir->coulombtype]);
1176 CHECK(ir->rcoulomb > ir->rlist);
1178 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1180 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1183 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist.\n"
1184 "For optimal energy conservation,consider using\n"
1185 "a potential modifier.", eel_names[ir->coulombtype]);
1186 CHECK(ir->rcoulomb != ir->rlist);
1191 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1193 // TODO: Move these checks into the ewald module with the options class
1195 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1197 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1199 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
1200 warning_error(wi, warn_buf);
1204 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1206 if (ir->ewald_geometry == eewg3D)
1208 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1209 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1210 warning(wi, warn_buf);
1212 /* This check avoids extra pbc coding for exclusion corrections */
1213 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1214 CHECK(ir->wall_ewald_zfac < 2);
1216 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1217 EEL_FULL(ir->coulombtype))
1219 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1220 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1221 warning(wi, warn_buf);
1223 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1225 if (ir->cutoff_scheme == ecutsVERLET)
1227 sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
1228 eel_names[ir->coulombtype]);
1229 warning(wi, warn_buf);
1233 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",
1234 eel_names[ir->coulombtype]);
1235 warning_note(wi, warn_buf);
1239 if (ir_vdw_switched(ir))
1241 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1242 CHECK(ir->rvdw_switch >= ir->rvdw);
1244 if (ir->rvdw_switch < 0.5*ir->rvdw)
1246 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.",
1247 ir->rvdw_switch, ir->rvdw);
1248 warning_note(wi, warn_buf);
1251 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1253 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1255 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1256 CHECK(ir->rlist > ir->rvdw);
1260 if (ir->vdwtype == evdwPME)
1262 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1264 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1265 evdw_names[ir->vdwtype],
1266 eintmod_names[eintmodPOTSHIFT],
1267 eintmod_names[eintmodNONE]);
1268 warning_error(wi, err_buf);
1272 if (ir->cutoff_scheme == ecutsGROUP)
1274 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1275 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)))
1277 warning_note(wi, "With exact cut-offs, rlist should be "
1278 "larger than rcoulomb and rvdw, so that there "
1279 "is a buffer region for particle motion "
1280 "between neighborsearch steps");
1283 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlist <= ir->rcoulomb)
1285 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rcoulomb.");
1286 warning_note(wi, warn_buf);
1288 if (ir_vdw_switched(ir) && (ir->rlist <= ir->rvdw))
1290 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rvdw.");
1291 warning_note(wi, warn_buf);
1295 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1297 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.");
1300 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1303 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1306 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1308 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1311 /* ENERGY CONSERVATION */
1312 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1314 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1316 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1317 evdw_names[evdwSHIFT]);
1318 warning_note(wi, warn_buf);
1320 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1322 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1323 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1324 warning_note(wi, warn_buf);
1328 /* IMPLICIT SOLVENT */
1329 if (ir->coulombtype == eelGB_NOTUSED)
1331 sprintf(warn_buf, "Invalid option %s for coulombtype",
1332 eel_names[ir->coulombtype]);
1333 warning_error(wi, warn_buf);
1336 if (ir->sa_algorithm == esaSTILL)
1338 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1339 CHECK(ir->sa_algorithm == esaSTILL);
1342 if (ir->implicit_solvent == eisGBSA)
1344 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1345 CHECK(ir->rgbradii != ir->rlist);
1347 if (ir->coulombtype != eelCUT)
1349 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1350 CHECK(ir->coulombtype != eelCUT);
1352 if (ir->vdwtype != evdwCUT)
1354 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1355 CHECK(ir->vdwtype != evdwCUT);
1357 if (ir->nstgbradii < 1)
1359 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1360 warning_note(wi, warn_buf);
1363 if (ir->sa_algorithm == esaNO)
1365 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1366 warning_note(wi, warn_buf);
1368 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1370 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");
1371 warning_note(wi, warn_buf);
1373 if (ir->gb_algorithm == egbSTILL)
1375 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1379 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1382 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1384 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1385 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1392 if (ir->cutoff_scheme != ecutsGROUP)
1394 warning_error(wi, "QMMM is currently only supported with cutoff-scheme=group");
1396 if (!EI_DYNAMICS(ir->eI))
1399 sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
1400 warning_error(wi, buf);
1406 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1410 /* count the number of text elemets separated by whitespace in a string.
1411 str = the input string
1412 maxptr = the maximum number of allowed elements
1413 ptr = the output array of pointers to the first character of each element
1414 returns: the number of elements. */
1415 int str_nelem(const char *str, int maxptr, char *ptr[])
1420 copy0 = gmx_strdup(str);
1423 while (*copy != '\0')
1427 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1435 while ((*copy != '\0') && !isspace(*copy))
1454 /* interpret a number of doubles from a string and put them in an array,
1455 after allocating space for them.
1456 str = the input string
1457 n = the (pre-allocated) number of doubles read
1458 r = the output array of doubles. */
1459 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1464 char warn_buf[STRLEN];
1466 *n = str_nelem(str, MAXPTR, ptr);
1469 for (i = 0; i < *n; i++)
1471 (*r)[i] = strtod(ptr[i], &endptr);
1474 sprintf(warn_buf, "Invalid value %s in string in mdp file. Expected a real number.", ptr[i]);
1475 warning_error(wi, warn_buf);
1480 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1483 int i, j, max_n_lambda, nweights, nfep[efptNR];
1484 t_lambda *fep = ir->fepvals;
1485 t_expanded *expand = ir->expandedvals;
1486 real **count_fep_lambdas;
1487 gmx_bool bOneLambda = TRUE;
1489 snew(count_fep_lambdas, efptNR);
1491 /* FEP input processing */
1492 /* first, identify the number of lambda values for each type.
1493 All that are nonzero must have the same number */
1495 for (i = 0; i < efptNR; i++)
1497 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1500 /* now, determine the number of components. All must be either zero, or equal. */
1503 for (i = 0; i < efptNR; i++)
1505 if (nfep[i] > max_n_lambda)
1507 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1508 must have the same number if its not zero.*/
1513 for (i = 0; i < efptNR; i++)
1517 ir->fepvals->separate_dvdl[i] = FALSE;
1519 else if (nfep[i] == max_n_lambda)
1521 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1522 respect to the temperature currently */
1524 ir->fepvals->separate_dvdl[i] = TRUE;
1529 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1530 nfep[i], efpt_names[i], max_n_lambda);
1533 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1534 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1536 /* the number of lambdas is the number we've read in, which is either zero
1537 or the same for all */
1538 fep->n_lambda = max_n_lambda;
1540 /* allocate space for the array of lambda values */
1541 snew(fep->all_lambda, efptNR);
1542 /* if init_lambda is defined, we need to set lambda */
1543 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1545 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1547 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1548 for (i = 0; i < efptNR; i++)
1550 snew(fep->all_lambda[i], fep->n_lambda);
1551 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1554 for (j = 0; j < fep->n_lambda; j++)
1556 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1558 sfree(count_fep_lambdas[i]);
1561 sfree(count_fep_lambdas);
1563 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1564 bookkeeping -- for now, init_lambda */
1566 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1568 for (i = 0; i < fep->n_lambda; i++)
1570 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1574 /* check to see if only a single component lambda is defined, and soft core is defined.
1575 In this case, turn on coulomb soft core */
1577 if (max_n_lambda == 0)
1583 for (i = 0; i < efptNR; i++)
1585 if ((nfep[i] != 0) && (i != efptFEP))
1591 if ((bOneLambda) && (fep->sc_alpha > 0))
1593 fep->bScCoul = TRUE;
1596 /* Fill in the others with the efptFEP if they are not explicitly
1597 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1598 they are all zero. */
1600 for (i = 0; i < efptNR; i++)
1602 if ((nfep[i] == 0) && (i != efptFEP))
1604 for (j = 0; j < fep->n_lambda; j++)
1606 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1612 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1613 if (fep->sc_r_power == 48)
1615 if (fep->sc_alpha > 0.1)
1617 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1621 /* now read in the weights */
1622 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1625 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1627 else if (nweights != fep->n_lambda)
1629 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1630 nweights, fep->n_lambda);
1632 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1634 expand->nstexpanded = fep->nstdhdl;
1635 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1637 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1639 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1640 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1641 2*tau_t just to be careful so it's not to frequent */
1646 static void do_simtemp_params(t_inputrec *ir)
1649 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1650 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1655 static void do_wall_params(t_inputrec *ir,
1656 char *wall_atomtype, char *wall_density,
1660 char *names[MAXPTR];
1663 opts->wall_atomtype[0] = nullptr;
1664 opts->wall_atomtype[1] = nullptr;
1666 ir->wall_atomtype[0] = -1;
1667 ir->wall_atomtype[1] = -1;
1668 ir->wall_density[0] = 0;
1669 ir->wall_density[1] = 0;
1673 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1674 if (nstr != ir->nwall)
1676 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1679 for (i = 0; i < ir->nwall; i++)
1681 opts->wall_atomtype[i] = gmx_strdup(names[i]);
1684 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1686 nstr = str_nelem(wall_density, MAXPTR, names);
1687 if (nstr != ir->nwall)
1689 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1691 for (i = 0; i < ir->nwall; i++)
1693 if (sscanf(names[i], "%lf", &dbl) != 1)
1695 gmx_fatal(FARGS, "Could not parse wall-density value from string '%s'", names[i]);
1699 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1701 ir->wall_density[i] = dbl;
1707 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1715 srenew(groups->grpname, groups->ngrpname+nwall);
1716 grps = &(groups->grps[egcENER]);
1717 srenew(grps->nm_ind, grps->nr+nwall);
1718 for (i = 0; i < nwall; i++)
1720 sprintf(str, "wall%d", i);
1721 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1722 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1727 static void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1728 t_expanded *expand, warninp_t wi)
1736 /* read expanded ensemble parameters */
1737 CCTYPE ("expanded ensemble variables");
1738 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1739 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1740 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1741 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1742 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1743 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1744 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1745 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1746 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1747 CCTYPE("Seed for Monte Carlo in lambda space");
1748 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1749 RTYPE ("mc-temperature", expand->mc_temp, -1);
1750 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1751 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1752 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1753 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1754 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1755 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1756 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1757 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1758 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1759 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1760 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1768 /*! \brief Return whether an end state with the given coupling-lambda
1769 * value describes fully-interacting VDW.
1771 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1772 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1774 static gmx_bool couple_lambda_has_vdw_on(int couple_lambda_value)
1776 return (couple_lambda_value == ecouplamVDW ||
1777 couple_lambda_value == ecouplamVDWQ);
1783 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1786 explicit MdpErrorHandler(warninp_t wi)
1787 : wi_(wi), mapping_(nullptr)
1791 void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
1793 mapping_ = &mapping;
1796 virtual bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context)
1798 ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
1799 getOptionName(context).c_str()));
1800 std::string message = gmx::formatExceptionMessageToString(*ex);
1801 warning_error(wi_, message.c_str());
1806 std::string getOptionName(const gmx::KeyValueTreePath &context)
1808 if (mapping_ != nullptr)
1810 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1811 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1814 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1819 const gmx::IKeyValueTreeBackMapping *mapping_;
1824 void get_ir(const char *mdparin, const char *mdparout,
1825 gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
1826 WriteMdpHeader writeMdpHeader, warninp_t wi)
1829 double dumdub[2][6];
1833 char warn_buf[STRLEN];
1834 t_lambda *fep = ir->fepvals;
1835 t_expanded *expand = ir->expandedvals;
1837 init_inputrec_strings();
1838 gmx::TextInputFile stream(mdparin);
1839 inp = read_inpfile(&stream, mdparin, &ninp, wi);
1841 snew(dumstr[0], STRLEN);
1842 snew(dumstr[1], STRLEN);
1844 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1847 "%s did not specify a value for the .mdp option "
1848 "\"cutoff-scheme\". Probably it was first intended for use "
1849 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1850 "introduced, but the group scheme was still the default. "
1851 "The default is now the Verlet scheme, so you will observe "
1852 "different behaviour.", mdparin);
1853 warning_note(wi, warn_buf);
1856 /* ignore the following deprecated commands */
1859 REM_TYPE("domain-decomposition");
1860 REM_TYPE("andersen-seed");
1862 REM_TYPE("dihre-fc");
1863 REM_TYPE("dihre-tau");
1864 REM_TYPE("nstdihreout");
1865 REM_TYPE("nstcheckpoint");
1866 REM_TYPE("optimize-fft");
1867 REM_TYPE("adress_type");
1868 REM_TYPE("adress_const_wf");
1869 REM_TYPE("adress_ex_width");
1870 REM_TYPE("adress_hy_width");
1871 REM_TYPE("adress_ex_forcecap");
1872 REM_TYPE("adress_interface_correction");
1873 REM_TYPE("adress_site");
1874 REM_TYPE("adress_reference_coords");
1875 REM_TYPE("adress_tf_grp_names");
1876 REM_TYPE("adress_cg_grp_names");
1877 REM_TYPE("adress_do_hybridpairs");
1878 REM_TYPE("rlistlong");
1879 REM_TYPE("nstcalclr");
1880 REM_TYPE("pull-print-com2");
1882 /* replace the following commands with the clearer new versions*/
1883 REPL_TYPE("unconstrained-start", "continuation");
1884 REPL_TYPE("foreign-lambda", "fep-lambdas");
1885 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1886 REPL_TYPE("nstxtcout", "nstxout-compressed");
1887 REPL_TYPE("xtc-grps", "compressed-x-grps");
1888 REPL_TYPE("xtc-precision", "compressed-x-precision");
1889 REPL_TYPE("pull-print-com1", "pull-print-com");
1891 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1892 CTYPE ("Preprocessor information: use cpp syntax.");
1893 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1894 STYPE ("include", opts->include, nullptr);
1895 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1896 STYPE ("define", opts->define, nullptr);
1898 CCTYPE ("RUN CONTROL PARAMETERS");
1899 EETYPE("integrator", ir->eI, ei_names);
1900 CTYPE ("Start time and timestep in ps");
1901 RTYPE ("tinit", ir->init_t, 0.0);
1902 RTYPE ("dt", ir->delta_t, 0.001);
1903 STEPTYPE ("nsteps", ir->nsteps, 0);
1904 CTYPE ("For exact run continuation or redoing part of a run");
1905 STEPTYPE ("init-step", ir->init_step, 0);
1906 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1907 ITYPE ("simulation-part", ir->simulation_part, 1);
1908 CTYPE ("mode for center of mass motion removal");
1909 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1910 CTYPE ("number of steps for center of mass motion removal");
1911 ITYPE ("nstcomm", ir->nstcomm, 100);
1912 CTYPE ("group(s) for center of mass motion removal");
1913 STYPE ("comm-grps", is->vcm, nullptr);
1915 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1916 CTYPE ("Friction coefficient (amu/ps) and random seed");
1917 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1918 STEPTYPE ("ld-seed", ir->ld_seed, -1);
1921 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1922 CTYPE ("Force tolerance and initial step-size");
1923 RTYPE ("emtol", ir->em_tol, 10.0);
1924 RTYPE ("emstep", ir->em_stepsize, 0.01);
1925 CTYPE ("Max number of iterations in relax-shells");
1926 ITYPE ("niter", ir->niter, 20);
1927 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1928 RTYPE ("fcstep", ir->fc_stepsize, 0);
1929 CTYPE ("Frequency of steepest descents steps when doing CG");
1930 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1931 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1933 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1934 RTYPE ("rtpi", ir->rtpi, 0.05);
1936 /* Output options */
1937 CCTYPE ("OUTPUT CONTROL OPTIONS");
1938 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1939 ITYPE ("nstxout", ir->nstxout, 0);
1940 ITYPE ("nstvout", ir->nstvout, 0);
1941 ITYPE ("nstfout", ir->nstfout, 0);
1942 CTYPE ("Output frequency for energies to log file and energy file");
1943 ITYPE ("nstlog", ir->nstlog, 1000);
1944 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1945 ITYPE ("nstenergy", ir->nstenergy, 1000);
1946 CTYPE ("Output frequency and precision for .xtc file");
1947 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
1948 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1949 CTYPE ("This selects the subset of atoms for the compressed");
1950 CTYPE ("trajectory file. You can select multiple groups. By");
1951 CTYPE ("default, all atoms will be written.");
1952 STYPE ("compressed-x-grps", is->x_compressed_groups, nullptr);
1953 CTYPE ("Selection of energy groups");
1954 STYPE ("energygrps", is->energy, nullptr);
1956 /* Neighbor searching */
1957 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1958 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1959 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1960 CTYPE ("nblist update frequency");
1961 ITYPE ("nstlist", ir->nstlist, 10);
1962 CTYPE ("ns algorithm (simple or grid)");
1963 EETYPE("ns-type", ir->ns_type, ens_names);
1964 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1965 EETYPE("pbc", ir->ePBC, epbc_names);
1966 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1967 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1968 CTYPE ("a value of -1 means: use rlist");
1969 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1970 CTYPE ("nblist cut-off");
1971 RTYPE ("rlist", ir->rlist, 1.0);
1972 CTYPE ("long-range cut-off for switched potentials");
1974 /* Electrostatics */
1975 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1976 CTYPE ("Method for doing electrostatics");
1977 EETYPE("coulombtype", ir->coulombtype, eel_names);
1978 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1979 CTYPE ("cut-off lengths");
1980 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1981 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1982 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1983 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1984 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1985 CTYPE ("Method for doing Van der Waals");
1986 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1987 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1988 CTYPE ("cut-off lengths");
1989 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1990 RTYPE ("rvdw", ir->rvdw, 1.0);
1991 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1992 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1993 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1994 RTYPE ("table-extension", ir->tabext, 1.0);
1995 CTYPE ("Separate tables between energy group pairs");
1996 STYPE ("energygrp-table", is->egptable, nullptr);
1997 CTYPE ("Spacing for the PME/PPPM FFT grid");
1998 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1999 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
2000 ITYPE ("fourier-nx", ir->nkx, 0);
2001 ITYPE ("fourier-ny", ir->nky, 0);
2002 ITYPE ("fourier-nz", ir->nkz, 0);
2003 CTYPE ("EWALD/PME/PPPM parameters");
2004 ITYPE ("pme-order", ir->pme_order, 4);
2005 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
2006 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
2007 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
2008 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
2009 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
2011 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
2012 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
2014 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
2015 CTYPE ("Algorithm for calculating Born radii");
2016 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
2017 CTYPE ("Frequency of calculating the Born radii inside rlist");
2018 ITYPE ("nstgbradii", ir->nstgbradii, 1);
2019 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
2020 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
2021 RTYPE ("rgbradii", ir->rgbradii, 1.0);
2022 CTYPE ("Dielectric coefficient of the implicit solvent");
2023 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
2024 CTYPE ("Salt concentration in M for Generalized Born models");
2025 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
2026 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
2027 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
2028 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
2029 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
2030 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
2031 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
2032 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
2033 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
2034 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
2036 /* Coupling stuff */
2037 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
2038 CTYPE ("Temperature coupling");
2039 EETYPE("tcoupl", ir->etc, etcoupl_names);
2040 ITYPE ("nsttcouple", ir->nsttcouple, -1);
2041 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
2042 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
2043 CTYPE ("Groups to couple separately");
2044 STYPE ("tc-grps", is->tcgrps, nullptr);
2045 CTYPE ("Time constant (ps) and reference temperature (K)");
2046 STYPE ("tau-t", is->tau_t, nullptr);
2047 STYPE ("ref-t", is->ref_t, nullptr);
2048 CTYPE ("pressure coupling");
2049 EETYPE("pcoupl", ir->epc, epcoupl_names);
2050 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
2051 ITYPE ("nstpcouple", ir->nstpcouple, -1);
2052 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
2053 RTYPE ("tau-p", ir->tau_p, 1.0);
2054 STYPE ("compressibility", dumstr[0], nullptr);
2055 STYPE ("ref-p", dumstr[1], nullptr);
2056 CTYPE ("Scaling of reference coordinates, No, All or COM");
2057 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
2060 CCTYPE ("OPTIONS FOR QMMM calculations");
2061 EETYPE("QMMM", ir->bQMMM, yesno_names);
2062 CTYPE ("Groups treated Quantum Mechanically");
2063 STYPE ("QMMM-grps", is->QMMM, nullptr);
2064 CTYPE ("QM method");
2065 STYPE("QMmethod", is->QMmethod, nullptr);
2066 CTYPE ("QMMM scheme");
2067 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
2068 CTYPE ("QM basisset");
2069 STYPE("QMbasis", is->QMbasis, nullptr);
2070 CTYPE ("QM charge");
2071 STYPE ("QMcharge", is->QMcharge, nullptr);
2072 CTYPE ("QM multiplicity");
2073 STYPE ("QMmult", is->QMmult, nullptr);
2074 CTYPE ("Surface Hopping");
2075 STYPE ("SH", is->bSH, nullptr);
2076 CTYPE ("CAS space options");
2077 STYPE ("CASorbitals", is->CASorbitals, nullptr);
2078 STYPE ("CASelectrons", is->CASelectrons, nullptr);
2079 STYPE ("SAon", is->SAon, nullptr);
2080 STYPE ("SAoff", is->SAoff, nullptr);
2081 STYPE ("SAsteps", is->SAsteps, nullptr);
2082 CTYPE ("Scale factor for MM charges");
2083 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
2085 /* Simulated annealing */
2086 CCTYPE("SIMULATED ANNEALING");
2087 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
2088 STYPE ("annealing", is->anneal, nullptr);
2089 CTYPE ("Number of time points to use for specifying annealing in each group");
2090 STYPE ("annealing-npoints", is->anneal_npoints, nullptr);
2091 CTYPE ("List of times at the annealing points for each group");
2092 STYPE ("annealing-time", is->anneal_time, nullptr);
2093 CTYPE ("Temp. at each annealing point, for each group.");
2094 STYPE ("annealing-temp", is->anneal_temp, nullptr);
2097 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
2098 EETYPE("gen-vel", opts->bGenVel, yesno_names);
2099 RTYPE ("gen-temp", opts->tempi, 300.0);
2100 ITYPE ("gen-seed", opts->seed, -1);
2103 CCTYPE ("OPTIONS FOR BONDS");
2104 EETYPE("constraints", opts->nshake, constraints);
2105 CTYPE ("Type of constraint algorithm");
2106 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
2107 CTYPE ("Do not constrain the start configuration");
2108 EETYPE("continuation", ir->bContinuation, yesno_names);
2109 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2110 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2111 CTYPE ("Relative tolerance of shake");
2112 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2113 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2114 ITYPE ("lincs-order", ir->nProjOrder, 4);
2115 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2116 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2117 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2118 ITYPE ("lincs-iter", ir->nLincsIter, 1);
2119 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2120 CTYPE ("rotates over more degrees than");
2121 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2122 CTYPE ("Convert harmonic bonds to morse potentials");
2123 EETYPE("morse", opts->bMorse, yesno_names);
2125 /* Energy group exclusions */
2126 CCTYPE ("ENERGY GROUP EXCLUSIONS");
2127 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2128 STYPE ("energygrp-excl", is->egpexcl, nullptr);
2132 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2133 ITYPE ("nwall", ir->nwall, 0);
2134 EETYPE("wall-type", ir->wall_type, ewt_names);
2135 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2136 STYPE ("wall-atomtype", is->wall_atomtype, nullptr);
2137 STYPE ("wall-density", is->wall_density, nullptr);
2138 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2141 CCTYPE("COM PULLING");
2142 EETYPE("pull", ir->bPull, yesno_names);
2146 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, wi);
2150 NOTE: needs COM pulling input */
2151 CCTYPE("AWH biasing");
2152 EETYPE("awh", ir->bDoAwh, yesno_names);
2157 ir->awhParams = gmx::readAndCheckAwhParams(&ninp, &inp, ir, wi);
2161 gmx_fatal(FARGS, "AWH biasing is only compatible with COM pulling turned on");
2165 /* Enforced rotation */
2166 CCTYPE("ENFORCED ROTATION");
2167 CTYPE("Enforced rotation: No or Yes");
2168 EETYPE("rotation", ir->bRot, yesno_names);
2172 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2175 /* Interactive MD */
2177 CCTYPE("Group to display and/or manipulate in interactive MD session");
2178 STYPE ("IMD-group", is->imd_grp, nullptr);
2179 if (is->imd_grp[0] != '\0')
2186 CCTYPE("NMR refinement stuff");
2187 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2188 EETYPE("disre", ir->eDisre, edisre_names);
2189 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2190 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2191 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2192 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2193 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2194 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2195 CTYPE ("Output frequency for pair distances to energy file");
2196 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2197 CTYPE ("Orientation restraints: No or Yes");
2198 EETYPE("orire", opts->bOrire, yesno_names);
2199 CTYPE ("Orientation restraints force constant and tau for time averaging");
2200 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2201 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2202 STYPE ("orire-fitgrp", is->orirefitgrp, nullptr);
2203 CTYPE ("Output frequency for trace(SD) and S to energy file");
2204 ITYPE ("nstorireout", ir->nstorireout, 100);
2206 /* free energy variables */
2207 CCTYPE ("Free energy variables");
2208 EETYPE("free-energy", ir->efep, efep_names);
2209 STYPE ("couple-moltype", is->couple_moltype, nullptr);
2210 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2211 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2212 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2214 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2216 it was not entered */
2217 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2218 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2219 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2220 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], nullptr);
2221 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], nullptr);
2222 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
2223 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
2224 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
2225 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
2226 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
2227 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2228 STYPE ("init-lambda-weights", is->lambda_weights, nullptr);
2229 EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
2230 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2231 ITYPE ("sc-power", fep->sc_power, 1);
2232 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2233 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2234 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2235 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2236 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2237 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2238 separate_dhdl_file_names);
2239 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2240 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2241 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2243 /* Non-equilibrium MD stuff */
2244 CCTYPE("Non-equilibrium MD stuff");
2245 STYPE ("acc-grps", is->accgrps, nullptr);
2246 STYPE ("accelerate", is->acc, nullptr);
2247 STYPE ("freezegrps", is->freeze, nullptr);
2248 STYPE ("freezedim", is->frdim, nullptr);
2249 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2250 STYPE ("deform", is->deform, nullptr);
2252 /* simulated tempering variables */
2253 CCTYPE("simulated tempering variables");
2254 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2255 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2256 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2257 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2259 /* expanded ensemble variables */
2260 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2262 read_expandedparams(&ninp, &inp, expand, wi);
2265 /* Electric fields */
2267 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(ninp, inp);
2268 gmx::KeyValueTreeTransformer transform;
2269 transform.rules()->addRule()
2270 .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2271 mdModules->initMdpTransform(transform.rules());
2272 for (const auto &path : transform.mappedPaths())
2274 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2275 mark_einp_set(ninp, inp, path[0].c_str());
2277 MdpErrorHandler errorHandler(wi);
2279 = transform.transform(convertedValues, &errorHandler);
2280 ir->params = new gmx::KeyValueTreeObject(result.object());
2281 mdModules->adjustInputrecBasedOnModules(ir);
2282 errorHandler.setBackMapping(result.backMapping());
2283 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2286 /* Ion/water position swapping ("computational electrophysiology") */
2287 CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2288 CTYPE("Swap positions along direction: no, X, Y, Z");
2289 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2290 if (ir->eSwapCoords != eswapNO)
2297 CTYPE("Swap attempt frequency");
2298 ITYPE("swap-frequency", ir->swap->nstswap, 1);
2299 CTYPE("Number of ion types to be controlled");
2300 ITYPE("iontypes", nIonTypes, 1);
2303 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2305 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2306 snew(ir->swap->grp, ir->swap->ngrp);
2307 for (i = 0; i < ir->swap->ngrp; i++)
2309 snew(ir->swap->grp[i].molname, STRLEN);
2311 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2312 STYPE("split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2313 STYPE("split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2314 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2315 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2316 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2318 CTYPE("Name of solvent molecules");
2319 STYPE("solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2321 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2322 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2323 CTYPE("however, if correctly defined, the permeation events are recorded per channel");
2324 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2325 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2326 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2327 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2328 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2329 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2331 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2332 ITYPE("coupl-steps", ir->swap->nAverage, 10);
2334 CTYPE("Names of the ion types that can be exchanged with solvent molecules,");
2335 CTYPE("and the requested number of ions of this type in compartments A and B");
2336 CTYPE("-1 means fix the numbers as found in step 0");
2337 for (i = 0; i < nIonTypes; i++)
2339 int ig = eSwapFixedGrpNR + i;
2341 sprintf(buf, "iontype%d-name", i);
2342 STYPE(buf, ir->swap->grp[ig].molname, nullptr);
2343 sprintf(buf, "iontype%d-in-A", i);
2344 ITYPE(buf, ir->swap->grp[ig].nmolReq[0], -1);
2345 sprintf(buf, "iontype%d-in-B", i);
2346 ITYPE(buf, ir->swap->grp[ig].nmolReq[1], -1);
2349 CTYPE("By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2350 CTYPE("at maximum distance (= bulk concentration) to the split group layers. However,");
2351 CTYPE("an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2352 CTYPE("towards one of the compartment-partitioning layers (at +/- 1.0).");
2353 RTYPE("bulk-offsetA", ir->swap->bulkOffset[0], 0.0);
2354 RTYPE("bulk-offsetB", ir->swap->bulkOffset[1], 0.0);
2355 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2356 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2358 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2361 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2362 RTYPE("threshold", ir->swap->threshold, 1.0);
2365 /* AdResS is no longer supported, but we need mdrun to be able to refuse to run old AdResS .tpr files */
2366 EETYPE("adress", ir->bAdress, yesno_names);
2368 /* User defined thingies */
2369 CCTYPE ("User defined thingies");
2370 STYPE ("user1-grps", is->user1, nullptr);
2371 STYPE ("user2-grps", is->user2, nullptr);
2372 ITYPE ("userint1", ir->userint1, 0);
2373 ITYPE ("userint2", ir->userint2, 0);
2374 ITYPE ("userint3", ir->userint3, 0);
2375 ITYPE ("userint4", ir->userint4, 0);
2376 RTYPE ("userreal1", ir->userreal1, 0);
2377 RTYPE ("userreal2", ir->userreal2, 0);
2378 RTYPE ("userreal3", ir->userreal3, 0);
2379 RTYPE ("userreal4", ir->userreal4, 0);
2383 gmx::TextOutputFile stream(mdparout);
2384 write_inpfile(&stream, mdparout, ninp, inp, FALSE, writeMdpHeader, wi);
2386 // Transform module data into a flat key-value tree for output.
2387 gmx::KeyValueTreeBuilder builder;
2388 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2389 mdModules->buildMdpOutput(&builderObject);
2391 gmx::TextWriter writer(&stream);
2392 writeKeyValueTreeAsMdp(&writer, builder.build());
2397 for (i = 0; (i < ninp); i++)
2400 sfree(inp[i].value);
2404 /* Process options if necessary */
2405 for (m = 0; m < 2; m++)
2407 for (i = 0; i < 2*DIM; i++)
2416 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2418 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
2420 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2422 case epctSEMIISOTROPIC:
2423 case epctSURFACETENSION:
2424 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2426 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
2428 dumdub[m][YY] = dumdub[m][XX];
2430 case epctANISOTROPIC:
2431 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2432 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2433 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2435 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
2439 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2440 epcoupltype_names[ir->epct]);
2444 clear_mat(ir->ref_p);
2445 clear_mat(ir->compress);
2446 for (i = 0; i < DIM; i++)
2448 ir->ref_p[i][i] = dumdub[1][i];
2449 ir->compress[i][i] = dumdub[0][i];
2451 if (ir->epct == epctANISOTROPIC)
2453 ir->ref_p[XX][YY] = dumdub[1][3];
2454 ir->ref_p[XX][ZZ] = dumdub[1][4];
2455 ir->ref_p[YY][ZZ] = dumdub[1][5];
2456 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2458 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2460 ir->compress[XX][YY] = dumdub[0][3];
2461 ir->compress[XX][ZZ] = dumdub[0][4];
2462 ir->compress[YY][ZZ] = dumdub[0][5];
2463 for (i = 0; i < DIM; i++)
2465 for (m = 0; m < i; m++)
2467 ir->ref_p[i][m] = ir->ref_p[m][i];
2468 ir->compress[i][m] = ir->compress[m][i];
2473 if (ir->comm_mode == ecmNO)
2478 opts->couple_moltype = nullptr;
2479 if (strlen(is->couple_moltype) > 0)
2481 if (ir->efep != efepNO)
2483 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2484 if (opts->couple_lam0 == opts->couple_lam1)
2486 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2488 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2489 opts->couple_lam1 == ecouplamNONE))
2491 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2496 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2499 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2500 if (ir->efep != efepNO)
2502 if (fep->delta_lambda > 0)
2504 ir->efep = efepSLOWGROWTH;
2508 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2510 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2511 warning_note(wi, "Old option for dhdl-print-energy given: "
2512 "changing \"yes\" to \"total\"\n");
2515 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2517 /* always print out the energy to dhdl if we are doing
2518 expanded ensemble, since we need the total energy for
2519 analysis if the temperature is changing. In some
2520 conditions one may only want the potential energy, so
2521 we will allow that if the appropriate mdp setting has
2522 been enabled. Otherwise, total it is:
2524 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2527 if ((ir->efep != efepNO) || ir->bSimTemp)
2529 ir->bExpanded = FALSE;
2530 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2532 ir->bExpanded = TRUE;
2534 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2535 if (ir->bSimTemp) /* done after fep params */
2537 do_simtemp_params(ir);
2540 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2541 * setup and not on the old way of specifying the free-energy setup,
2542 * we should check for using soft-core when not needed, since that
2543 * can complicate the sampling significantly.
2544 * Note that we only check for the automated coupling setup.
2545 * If the (advanced) user does FEP through manual topology changes,
2546 * this check will not be triggered.
2548 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2549 ir->fepvals->sc_alpha != 0 &&
2550 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2551 couple_lambda_has_vdw_on(opts->couple_lam1)))
2553 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.");
2558 ir->fepvals->n_lambda = 0;
2561 /* WALL PARAMETERS */
2563 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2565 /* ORIENTATION RESTRAINT PARAMETERS */
2567 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, nullptr) != 1)
2569 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2572 /* DEFORMATION PARAMETERS */
2574 clear_mat(ir->deform);
2575 for (i = 0; i < 6; i++)
2580 double gmx_unused canary;
2581 int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
2582 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2583 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
2585 if (strlen(is->deform) > 0 && ndeform != 6)
2587 warning_error(wi, gmx::formatString("Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform).c_str());
2589 for (i = 0; i < 3; i++)
2591 ir->deform[i][i] = dumdub[0][i];
2593 ir->deform[YY][XX] = dumdub[0][3];
2594 ir->deform[ZZ][XX] = dumdub[0][4];
2595 ir->deform[ZZ][YY] = dumdub[0][5];
2596 if (ir->epc != epcNO)
2598 for (i = 0; i < 3; i++)
2600 for (j = 0; j <= i; j++)
2602 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2604 warning_error(wi, "A box element has deform set and compressibility > 0");
2608 for (i = 0; i < 3; i++)
2610 for (j = 0; j < i; j++)
2612 if (ir->deform[i][j] != 0)
2614 for (m = j; m < DIM; m++)
2616 if (ir->compress[m][j] != 0)
2618 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.");
2619 warning(wi, warn_buf);
2627 /* Ion/water position swapping checks */
2628 if (ir->eSwapCoords != eswapNO)
2630 if (ir->swap->nstswap < 1)
2632 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2634 if (ir->swap->nAverage < 1)
2636 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2638 if (ir->swap->threshold < 1.0)
2640 warning_error(wi, "Ion count threshold must be at least 1.\n");
2648 static int search_QMstring(const char *s, int ng, const char *gn[])
2650 /* same as normal search_string, but this one searches QM strings */
2653 for (i = 0; (i < ng); i++)
2655 if (gmx_strcasecmp(s, gn[i]) == 0)
2661 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2665 } /* search_QMstring */
2667 /* We would like gn to be const as well, but C doesn't allow this */
2668 /* TODO this is utility functionality (search for the index of a
2669 string in a collection), so should be refactored and located more
2671 int search_string(const char *s, int ng, char *gn[])
2675 for (i = 0; (i < ng); i++)
2677 if (gmx_strcasecmp(s, gn[i]) == 0)
2684 "Group %s referenced in the .mdp file was not found in the index file.\n"
2685 "Group names must match either [moleculetype] names or custom index group\n"
2686 "names, in which case you must supply an index file to the '-n' option\n"
2693 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2694 t_blocka *block, char *gnames[],
2695 int gtype, int restnm,
2696 int grptp, gmx_bool bVerbose,
2699 unsigned short *cbuf;
2700 t_grps *grps = &(groups->grps[gtype]);
2701 int i, j, gid, aj, ognr, ntot = 0;
2704 char warn_buf[STRLEN];
2708 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2711 title = gtypes[gtype];
2714 /* Mark all id's as not set */
2715 for (i = 0; (i < natoms); i++)
2720 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2721 for (i = 0; (i < ng); i++)
2723 /* Lookup the group name in the block structure */
2724 gid = search_string(ptrs[i], block->nr, gnames);
2725 if ((grptp != egrptpONE) || (i == 0))
2727 grps->nm_ind[grps->nr++] = gid;
2731 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2734 /* Now go over the atoms in the group */
2735 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2740 /* Range checking */
2741 if ((aj < 0) || (aj >= natoms))
2743 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2745 /* Lookup up the old group number */
2749 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2750 aj+1, title, ognr+1, i+1);
2754 /* Store the group number in buffer */
2755 if (grptp == egrptpONE)
2768 /* Now check whether we have done all atoms */
2772 if (grptp == egrptpALL)
2774 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2775 natoms-ntot, title);
2777 else if (grptp == egrptpPART)
2779 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2780 natoms-ntot, title);
2781 warning_note(wi, warn_buf);
2783 /* Assign all atoms currently unassigned to a rest group */
2784 for (j = 0; (j < natoms); j++)
2786 if (cbuf[j] == NOGID)
2792 if (grptp != egrptpPART)
2797 "Making dummy/rest group for %s containing %d elements\n",
2798 title, natoms-ntot);
2800 /* Add group name "rest" */
2801 grps->nm_ind[grps->nr] = restnm;
2803 /* Assign the rest name to all atoms not currently assigned to a group */
2804 for (j = 0; (j < natoms); j++)
2806 if (cbuf[j] == NOGID)
2815 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2817 /* All atoms are part of one (or no) group, no index required */
2818 groups->ngrpnr[gtype] = 0;
2819 groups->grpnr[gtype] = nullptr;
2823 groups->ngrpnr[gtype] = natoms;
2824 snew(groups->grpnr[gtype], natoms);
2825 for (j = 0; (j < natoms); j++)
2827 groups->grpnr[gtype][j] = cbuf[j];
2833 return (bRest && grptp == egrptpPART);
2836 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2839 gmx_groups_t *groups;
2840 pull_params_t *pull;
2841 int natoms, ai, aj, i, j, d, g, imin, jmin;
2843 int *nrdf2, *na_vcm, na_tot;
2844 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2846 gmx_mtop_atomloop_all_t aloop;
2847 int mb, mol, ftype, as;
2848 gmx_molblock_t *molb;
2849 gmx_moltype_t *molt;
2852 * First calc 3xnr-atoms for each group
2853 * then subtract half a degree of freedom for each constraint
2855 * Only atoms and nuclei contribute to the degrees of freedom...
2860 groups = &mtop->groups;
2861 natoms = mtop->natoms;
2863 /* Allocate one more for a possible rest group */
2864 /* We need to sum degrees of freedom into doubles,
2865 * since floats give too low nrdf's above 3 million atoms.
2867 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2868 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2869 snew(dof_vcm, groups->grps[egcVCM].nr+1);
2870 snew(na_vcm, groups->grps[egcVCM].nr+1);
2871 snew(nrdf_vcm_sub, groups->grps[egcVCM].nr+1);
2873 for (i = 0; i < groups->grps[egcTC].nr; i++)
2877 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2880 clear_ivec(dof_vcm[i]);
2882 nrdf_vcm_sub[i] = 0;
2885 snew(nrdf2, natoms);
2886 aloop = gmx_mtop_atomloop_all_init(mtop);
2888 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2891 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2893 g = ggrpnr(groups, egcFREEZE, i);
2894 for (d = 0; d < DIM; d++)
2896 if (opts->nFreeze[g][d] == 0)
2898 /* Add one DOF for particle i (counted as 2*1) */
2900 /* VCM group i has dim d as a DOF */
2901 dof_vcm[ggrpnr(groups, egcVCM, i)][d] = 1;
2904 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2905 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2910 for (mb = 0; mb < mtop->nmolblock; mb++)
2912 molb = &mtop->molblock[mb];
2913 molt = &mtop->moltype[molb->type];
2914 atom = molt->atoms.atom;
2915 for (mol = 0; mol < molb->nmol; mol++)
2917 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2919 ia = molt->ilist[ftype].iatoms;
2920 for (i = 0; i < molt->ilist[ftype].nr; )
2922 /* Subtract degrees of freedom for the constraints,
2923 * if the particles still have degrees of freedom left.
2924 * If one of the particles is a vsite or a shell, then all
2925 * constraint motion will go there, but since they do not
2926 * contribute to the constraints the degrees of freedom do not
2931 if (((atom[ia[1]].ptype == eptNucleus) ||
2932 (atom[ia[1]].ptype == eptAtom)) &&
2933 ((atom[ia[2]].ptype == eptNucleus) ||
2934 (atom[ia[2]].ptype == eptAtom)))
2952 imin = std::min(imin, nrdf2[ai]);
2953 jmin = std::min(jmin, nrdf2[aj]);
2956 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2957 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2958 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2959 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2961 ia += interaction_function[ftype].nratoms+1;
2962 i += interaction_function[ftype].nratoms+1;
2965 ia = molt->ilist[F_SETTLE].iatoms;
2966 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2968 /* Subtract 1 dof from every atom in the SETTLE */
2969 for (j = 0; j < 3; j++)
2972 imin = std::min(2, nrdf2[ai]);
2974 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2975 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2980 as += molt->atoms.nr;
2986 /* Correct nrdf for the COM constraints.
2987 * We correct using the TC and VCM group of the first atom
2988 * in the reference and pull group. If atoms in one pull group
2989 * belong to different TC or VCM groups it is anyhow difficult
2990 * to determine the optimal nrdf assignment.
2994 for (i = 0; i < pull->ncoord; i++)
2996 if (pull->coord[i].eType != epullCONSTRAINT)
3003 for (j = 0; j < 2; j++)
3005 const t_pull_group *pgrp;
3007 pgrp = &pull->group[pull->coord[i].group[j]];
3011 /* Subtract 1/2 dof from each group */
3013 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
3014 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
3015 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
3017 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)]]);
3022 /* We need to subtract the whole DOF from group j=1 */
3029 if (ir->nstcomm != 0)
3033 /* We remove COM motion up to dim ndof_com() */
3034 ndim_rm_vcm = ndof_com(ir);
3036 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3037 * the number of degrees of freedom in each vcm group when COM
3038 * translation is removed and 6 when rotation is removed as well.
3040 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3042 switch (ir->comm_mode)
3045 case ecmLINEAR_ACCELERATION_CORRECTION:
3046 nrdf_vcm_sub[j] = 0;
3047 for (d = 0; d < ndim_rm_vcm; d++)
3056 nrdf_vcm_sub[j] = 6;
3059 gmx_incons("Checking comm_mode");
3063 for (i = 0; i < groups->grps[egcTC].nr; i++)
3065 /* Count the number of atoms of TC group i for every VCM group */
3066 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3071 for (ai = 0; ai < natoms; ai++)
3073 if (ggrpnr(groups, egcTC, ai) == i)
3075 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
3079 /* Correct for VCM removal according to the fraction of each VCM
3080 * group present in this TC group.
3082 nrdf_uc = nrdf_tc[i];
3085 fprintf(debug, "T-group[%d] nrdf_uc = %g\n", i, nrdf_uc);
3088 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3090 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3092 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
3093 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
3097 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
3098 j, nrdf_vcm[j], nrdf_tc[i]);
3103 for (i = 0; (i < groups->grps[egcTC].nr); i++)
3105 opts->nrdf[i] = nrdf_tc[i];
3106 if (opts->nrdf[i] < 0)
3111 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3112 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
3120 sfree(nrdf_vcm_sub);
3123 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
3124 const char *option, const char *val, int flag)
3126 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3127 * But since this is much larger than STRLEN, such a line can not be parsed.
3128 * The real maximum is the number of names that fit in a string: STRLEN/2.
3130 #define EGP_MAX (STRLEN/2)
3131 int nelem, i, j, k, nr;
3132 char *names[EGP_MAX];
3136 gnames = groups->grpname;
3138 nelem = str_nelem(val, EGP_MAX, names);
3141 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3143 nr = groups->grps[egcENER].nr;
3145 for (i = 0; i < nelem/2; i++)
3149 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
3155 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3156 names[2*i], option);
3160 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3166 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3167 names[2*i+1], option);
3169 if ((j < nr) && (k < nr))
3171 ir->opts.egp_flags[nr*j+k] |= flag;
3172 ir->opts.egp_flags[nr*k+j] |= flag;
3181 static void make_swap_groups(
3186 int ig = -1, i = 0, gind;
3190 /* Just a quick check here, more thorough checks are in mdrun */
3191 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3193 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3196 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3197 for (ig = 0; ig < swap->ngrp; ig++)
3199 swapg = &swap->grp[ig];
3200 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3201 swapg->nat = grps->index[gind+1] - grps->index[gind];
3205 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3206 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3207 swap->grp[ig].molname, swapg->nat);
3208 snew(swapg->ind, swapg->nat);
3209 for (i = 0; i < swapg->nat; i++)
3211 swapg->ind[i] = grps->a[grps->index[gind]+i];
3216 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3222 static void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3227 ig = search_string(IMDgname, grps->nr, gnames);
3228 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3230 if (IMDgroup->nat > 0)
3232 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3233 IMDgname, IMDgroup->nat);
3234 snew(IMDgroup->ind, IMDgroup->nat);
3235 for (i = 0; i < IMDgroup->nat; i++)
3237 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3243 void do_index(const char* mdparin, const char *ndx,
3250 gmx_groups_t *groups;
3254 char warnbuf[STRLEN], **gnames;
3255 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3258 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3259 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3260 int i, j, k, restnm;
3261 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3262 int nQMmethod, nQMbasis, nQMg;
3263 char warn_buf[STRLEN];
3268 fprintf(stderr, "processing index file...\n");
3273 snew(grps->index, 1);
3275 atoms_all = gmx_mtop_global_atoms(mtop);
3276 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3277 done_atom(&atoms_all);
3281 grps = init_index(ndx, &gnames);
3284 groups = &mtop->groups;
3285 natoms = mtop->natoms;
3286 symtab = &mtop->symtab;
3288 snew(groups->grpname, grps->nr+1);
3290 for (i = 0; (i < grps->nr); i++)
3292 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3294 groups->grpname[i] = put_symtab(symtab, "rest");
3296 srenew(gnames, grps->nr+1);
3297 gnames[restnm] = *(groups->grpname[i]);
3298 groups->ngrpname = grps->nr+1;
3300 set_warning_line(wi, mdparin, -1);
3302 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3303 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3304 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3305 if ((ntau_t != ntcg) || (nref_t != ntcg))
3307 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3308 "%d tau-t values", ntcg, nref_t, ntau_t);
3311 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3312 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3313 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3314 nr = groups->grps[egcTC].nr;
3316 snew(ir->opts.nrdf, nr);
3317 snew(ir->opts.tau_t, nr);
3318 snew(ir->opts.ref_t, nr);
3319 if (ir->eI == eiBD && ir->bd_fric == 0)
3321 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3328 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3332 for (i = 0; (i < nr); i++)
3334 ir->opts.tau_t[i] = strtod(ptr1[i], &endptr);
3337 warning_error(wi, "Invalid value for mdp option tau-t. tau-t should only consist of real numbers separated by spaces.");
3339 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3341 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3342 warning_error(wi, warn_buf);
3345 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3347 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.");
3350 if (ir->opts.tau_t[i] >= 0)
3352 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3355 if (ir->etc != etcNO && ir->nsttcouple == -1)
3357 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3362 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3364 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");
3366 if (ir->epc == epcMTTK)
3368 if (ir->etc != etcNOSEHOOVER)
3370 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3374 if (ir->nstpcouple != ir->nsttcouple)
3376 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3377 ir->nstpcouple = ir->nsttcouple = mincouple;
3378 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);
3379 warning_note(wi, warn_buf);
3384 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3385 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3387 if (ETC_ANDERSEN(ir->etc))
3389 if (ir->nsttcouple != 1)
3392 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");
3393 warning_note(wi, warn_buf);
3396 nstcmin = tcouple_min_integration_steps(ir->etc);
3399 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3401 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3402 ETCOUPLTYPE(ir->etc),
3404 ir->nsttcouple*ir->delta_t);
3405 warning(wi, warn_buf);
3408 for (i = 0; (i < nr); i++)
3410 ir->opts.ref_t[i] = strtod(ptr2[i], &endptr);
3413 warning_error(wi, "Invalid value for mdp option ref-t. ref-t should only consist of real numbers separated by spaces.");
3415 if (ir->opts.ref_t[i] < 0)
3417 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3420 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3421 if we are in this conditional) if mc_temp is negative */
3422 if (ir->expandedvals->mc_temp < 0)
3424 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3428 /* Simulated annealing for each group. There are nr groups */
3429 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3430 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3434 if (nSA > 0 && nSA != nr)
3436 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3440 snew(ir->opts.annealing, nr);
3441 snew(ir->opts.anneal_npoints, nr);
3442 snew(ir->opts.anneal_time, nr);
3443 snew(ir->opts.anneal_temp, nr);
3444 for (i = 0; i < nr; i++)
3446 ir->opts.annealing[i] = eannNO;
3447 ir->opts.anneal_npoints[i] = 0;
3448 ir->opts.anneal_time[i] = nullptr;
3449 ir->opts.anneal_temp[i] = nullptr;
3454 for (i = 0; i < nr; i++)
3456 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3458 ir->opts.annealing[i] = eannNO;
3460 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3462 ir->opts.annealing[i] = eannSINGLE;
3465 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3467 ir->opts.annealing[i] = eannPERIODIC;
3473 /* Read the other fields too */
3474 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3475 if (nSA_points != nSA)
3477 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3479 for (k = 0, i = 0; i < nr; i++)
3481 ir->opts.anneal_npoints[i] = strtol(ptr1[i], &endptr, 10);
3484 warning_error(wi, "Invalid value for mdp option annealing-npoints. annealing should only consist of integers separated by spaces.");
3486 if (ir->opts.anneal_npoints[i] == 1)
3488 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3490 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3491 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3492 k += ir->opts.anneal_npoints[i];
3495 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3498 gmx_fatal(FARGS, "Found %d annealing-time values, wanted %d\n", nSA_time, k);
3500 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3503 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3506 for (i = 0, k = 0; i < nr; i++)
3509 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3511 ir->opts.anneal_time[i][j] = strtod(ptr1[k], &endptr);
3514 warning_error(wi, "Invalid value for mdp option anneal-time. anneal-time should only consist of real numbers separated by spaces.");
3516 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], &endptr);
3519 warning_error(wi, "Invalid value for anneal-temp. anneal-temp should only consist of real numbers separated by spaces.");
3523 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3525 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3531 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3533 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3534 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3537 if (ir->opts.anneal_temp[i][j] < 0)
3539 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3544 /* Print out some summary information, to make sure we got it right */
3545 for (i = 0, k = 0; i < nr; i++)
3547 if (ir->opts.annealing[i] != eannNO)
3549 j = groups->grps[egcTC].nm_ind[i];
3550 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3551 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3552 ir->opts.anneal_npoints[i]);
3553 fprintf(stderr, "Time (ps) Temperature (K)\n");
3554 /* All terms except the last one */
3555 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3557 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3560 /* Finally the last one */
3561 j = ir->opts.anneal_npoints[i]-1;
3562 if (ir->opts.annealing[i] == eannSINGLE)
3564 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3568 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3569 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3571 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3582 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3584 make_pull_coords(ir->pull);
3589 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3592 if (ir->eSwapCoords != eswapNO)
3594 make_swap_groups(ir->swap, grps, gnames);
3597 /* Make indices for IMD session */
3600 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3603 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3604 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3605 if (nacg*DIM != nacc)
3607 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3610 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3611 restnm, egrptpALL_GENREST, bVerbose, wi);
3612 nr = groups->grps[egcACC].nr;
3613 snew(ir->opts.acc, nr);
3614 ir->opts.ngacc = nr;
3616 for (i = k = 0; (i < nacg); i++)
3618 for (j = 0; (j < DIM); j++, k++)
3620 ir->opts.acc[i][j] = strtod(ptr1[k], &endptr);
3623 warning_error(wi, "Invalid value for mdp option accelerate. accelerate should only consist of real numbers separated by spaces.");
3627 for (; (i < nr); i++)
3629 for (j = 0; (j < DIM); j++)
3631 ir->opts.acc[i][j] = 0;
3635 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3636 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3637 if (nfrdim != DIM*nfreeze)
3639 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3642 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3643 restnm, egrptpALL_GENREST, bVerbose, wi);
3644 nr = groups->grps[egcFREEZE].nr;
3645 ir->opts.ngfrz = nr;
3646 snew(ir->opts.nFreeze, nr);
3647 for (i = k = 0; (i < nfreeze); i++)
3649 for (j = 0; (j < DIM); j++, k++)
3651 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3652 if (!ir->opts.nFreeze[i][j])
3654 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3656 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3657 "(not %s)", ptr1[k]);
3658 warning(wi, warn_buf);
3663 for (; (i < nr); i++)
3665 for (j = 0; (j < DIM); j++)
3667 ir->opts.nFreeze[i][j] = 0;
3671 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3672 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3673 restnm, egrptpALL_GENREST, bVerbose, wi);
3674 add_wall_energrps(groups, ir->nwall, symtab);
3675 ir->opts.ngener = groups->grps[egcENER].nr;
3676 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3678 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3679 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3682 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3683 "This may lead to artifacts.\n"
3684 "In most cases one should use one group for the whole system.");
3687 /* Now we have filled the freeze struct, so we can calculate NRDF */
3688 calc_nrdf(mtop, ir, gnames);
3690 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3691 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3692 restnm, egrptpALL_GENREST, bVerbose, wi);
3693 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3694 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3695 restnm, egrptpALL_GENREST, bVerbose, wi);
3696 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3697 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3698 restnm, egrptpONE, bVerbose, wi);
3699 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3700 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3701 restnm, egrptpALL_GENREST, bVerbose, wi);
3703 /* QMMM input processing */
3704 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3705 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3706 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3707 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3709 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3710 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3712 /* group rest, if any, is always MM! */
3713 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3714 restnm, egrptpALL_GENREST, bVerbose, wi);
3715 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3716 ir->opts.ngQM = nQMg;
3717 snew(ir->opts.QMmethod, nr);
3718 snew(ir->opts.QMbasis, nr);
3719 for (i = 0; i < nr; i++)
3721 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3722 * converted to the corresponding enum in names.c
3724 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3726 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3730 str_nelem(is->QMmult, MAXPTR, ptr1);
3731 str_nelem(is->QMcharge, MAXPTR, ptr2);
3732 str_nelem(is->bSH, MAXPTR, ptr3);
3733 snew(ir->opts.QMmult, nr);
3734 snew(ir->opts.QMcharge, nr);
3735 snew(ir->opts.bSH, nr);
3737 for (i = 0; i < nr; i++)
3739 ir->opts.QMmult[i] = strtol(ptr1[i], &endptr, 10);
3742 warning_error(wi, "Invalid value for mdp option QMmult. QMmult should only consist of integers separated by spaces.");
3744 ir->opts.QMcharge[i] = strtol(ptr2[i], &endptr, 10);
3747 warning_error(wi, "Invalid value for mdp option QMcharge. QMcharge should only consist of integers separated by spaces.");
3749 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3752 str_nelem(is->CASelectrons, MAXPTR, ptr1);
3753 str_nelem(is->CASorbitals, MAXPTR, ptr2);
3754 snew(ir->opts.CASelectrons, nr);
3755 snew(ir->opts.CASorbitals, nr);
3756 for (i = 0; i < nr; i++)
3758 ir->opts.CASelectrons[i] = strtol(ptr1[i], &endptr, 10);
3761 warning_error(wi, "Invalid value for mdp option CASelectrons. CASelectrons should only consist of integers separated by spaces.");
3763 ir->opts.CASorbitals[i] = strtol(ptr2[i], &endptr, 10);
3766 warning_error(wi, "Invalid value for mdp option CASorbitals. CASorbitals should only consist of integers separated by spaces.");
3770 str_nelem(is->SAon, MAXPTR, ptr1);
3771 str_nelem(is->SAoff, MAXPTR, ptr2);
3772 str_nelem(is->SAsteps, MAXPTR, ptr3);
3773 snew(ir->opts.SAon, nr);
3774 snew(ir->opts.SAoff, nr);
3775 snew(ir->opts.SAsteps, nr);
3777 for (i = 0; i < nr; i++)
3779 ir->opts.SAon[i] = strtod(ptr1[i], &endptr);
3782 warning_error(wi, "Invalid value for mdp option SAon. SAon should only consist of real numbers separated by spaces.");
3784 ir->opts.SAoff[i] = strtod(ptr2[i], &endptr);
3787 warning_error(wi, "Invalid value for mdp option SAoff. SAoff should only consist of real numbers separated by spaces.");
3789 ir->opts.SAsteps[i] = strtol(ptr3[i], &endptr, 10);
3792 warning_error(wi, "Invalid value for mdp option SAsteps. SAsteps should only consist of integers separated by spaces.");
3795 /* end of QMMM input */
3799 for (i = 0; (i < egcNR); i++)
3801 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3802 for (j = 0; (j < groups->grps[i].nr); j++)
3804 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3806 fprintf(stderr, "\n");
3810 nr = groups->grps[egcENER].nr;
3811 snew(ir->opts.egp_flags, nr*nr);
3813 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3814 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3816 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3818 if (bExcl && EEL_FULL(ir->coulombtype))
3820 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3823 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3824 if (bTable && !(ir->vdwtype == evdwUSER) &&
3825 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3826 !(ir->coulombtype == eelPMEUSERSWITCH))
3828 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3831 for (i = 0; (i < grps->nr); i++)
3843 static void check_disre(gmx_mtop_t *mtop)
3845 gmx_ffparams_t *ffparams;
3846 t_functype *functype;
3848 int i, ndouble, ftype;
3849 int label, old_label;
3851 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3853 ffparams = &mtop->ffparams;
3854 functype = ffparams->functype;
3855 ip = ffparams->iparams;
3858 for (i = 0; i < ffparams->ntypes; i++)
3860 ftype = functype[i];
3861 if (ftype == F_DISRES)
3863 label = ip[i].disres.label;
3864 if (label == old_label)
3866 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3874 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3875 "probably the parameters for multiple pairs in one restraint "
3876 "are not identical\n", ndouble);
3881 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3882 gmx_bool posres_only,
3886 gmx_mtop_ilistloop_t iloop;
3896 for (d = 0; d < DIM; d++)
3898 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3900 /* Check for freeze groups */
3901 for (g = 0; g < ir->opts.ngfrz; g++)
3903 for (d = 0; d < DIM; d++)
3905 if (ir->opts.nFreeze[g][d] != 0)
3913 /* Check for position restraints */
3914 iloop = gmx_mtop_ilistloop_init(sys);
3915 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3918 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3920 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3922 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3923 for (d = 0; d < DIM; d++)
3925 if (pr->posres.fcA[d] != 0)
3931 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3933 /* Check for flat-bottom posres */
3934 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3935 if (pr->fbposres.k != 0)
3937 switch (pr->fbposres.geom)
3939 case efbposresSPHERE:
3940 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3942 case efbposresCYLINDERX:
3943 AbsRef[YY] = AbsRef[ZZ] = 1;
3945 case efbposresCYLINDERY:
3946 AbsRef[XX] = AbsRef[ZZ] = 1;
3948 case efbposresCYLINDER:
3949 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3950 case efbposresCYLINDERZ:
3951 AbsRef[XX] = AbsRef[YY] = 1;
3953 case efbposresX: /* d=XX */
3954 case efbposresY: /* d=YY */
3955 case efbposresZ: /* d=ZZ */
3956 d = pr->fbposres.geom - efbposresX;
3960 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3961 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3969 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3973 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3974 gmx_bool *bC6ParametersWorkWithGeometricRules,
3975 gmx_bool *bC6ParametersWorkWithLBRules,
3976 gmx_bool *bLBRulesPossible)
3978 int ntypes, tpi, tpj;
3981 double c6i, c6j, c12i, c12j;
3982 double c6, c6_geometric, c6_LB;
3983 double sigmai, sigmaj, epsi, epsj;
3984 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3987 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3988 * force-field floating point parameters.
3991 ptr = getenv("GMX_LJCOMB_TOL");
3995 double gmx_unused canary;
3997 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3999 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
4004 *bC6ParametersWorkWithLBRules = TRUE;
4005 *bC6ParametersWorkWithGeometricRules = TRUE;
4006 bCanDoLBRules = TRUE;
4007 ntypes = mtop->ffparams.atnr;
4008 snew(typecount, ntypes);
4009 gmx_mtop_count_atomtypes(mtop, state, typecount);
4010 *bLBRulesPossible = TRUE;
4011 for (tpi = 0; tpi < ntypes; ++tpi)
4013 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
4014 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
4015 for (tpj = tpi; tpj < ntypes; ++tpj)
4017 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
4018 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
4019 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
4020 c6_geometric = std::sqrt(c6i * c6j);
4021 if (!gmx_numzero(c6_geometric))
4023 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4025 sigmai = gmx::sixthroot(c12i / c6i);
4026 sigmaj = gmx::sixthroot(c12j / c6j);
4027 epsi = c6i * c6i /(4.0 * c12i);
4028 epsj = c6j * c6j /(4.0 * c12j);
4029 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4033 *bLBRulesPossible = FALSE;
4034 c6_LB = c6_geometric;
4036 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4039 if (FALSE == bCanDoLBRules)
4041 *bC6ParametersWorkWithLBRules = FALSE;
4044 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4046 if (FALSE == bCanDoGeometricRules)
4048 *bC6ParametersWorkWithGeometricRules = FALSE;
4056 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
4059 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4061 check_combination_rule_differences(mtop, 0,
4062 &bC6ParametersWorkWithGeometricRules,
4063 &bC6ParametersWorkWithLBRules,
4065 if (ir->ljpme_combination_rule == eljpmeLB)
4067 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
4069 warning(wi, "You are using arithmetic-geometric combination rules "
4070 "in LJ-PME, but your non-bonded C6 parameters do not "
4071 "follow these rules.");
4076 if (FALSE == bC6ParametersWorkWithGeometricRules)
4078 if (ir->eDispCorr != edispcNO)
4080 warning_note(wi, "You are using geometric combination rules in "
4081 "LJ-PME, but your non-bonded C6 parameters do "
4082 "not follow these rules. "
4083 "This will introduce very small errors in the forces and energies in "
4084 "your simulations. Dispersion correction will correct total energy "
4085 "and/or pressure for isotropic systems, but not forces or surface tensions.");
4089 warning_note(wi, "You are using geometric combination rules in "
4090 "LJ-PME, but your non-bonded C6 parameters do "
4091 "not follow these rules. "
4092 "This will introduce very small errors in the forces and energies in "
4093 "your simulations. If your system is homogeneous, consider using dispersion correction "
4094 "for the total energy and pressure.");
4100 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
4103 char err_buf[STRLEN];
4105 gmx_bool bCharge, bAcc;
4108 gmx_mtop_atomloop_block_t aloopb;
4109 gmx_mtop_atomloop_all_t aloop;
4111 char warn_buf[STRLEN];
4113 set_warning_line(wi, mdparin, -1);
4115 if (ir->cutoff_scheme == ecutsVERLET &&
4116 ir->verletbuf_tol > 0 &&
4118 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
4119 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4121 /* Check if a too small Verlet buffer might potentially
4122 * cause more drift than the thermostat can couple off.
4124 /* Temperature error fraction for warning and suggestion */
4125 const real T_error_warn = 0.002;
4126 const real T_error_suggest = 0.001;
4127 /* For safety: 2 DOF per atom (typical with constraints) */
4128 const real nrdf_at = 2;
4129 real T, tau, max_T_error;
4134 for (i = 0; i < ir->opts.ngtc; i++)
4136 T = std::max(T, ir->opts.ref_t[i]);
4137 tau = std::max(tau, ir->opts.tau_t[i]);
4141 /* This is a worst case estimate of the temperature error,
4142 * assuming perfect buffer estimation and no cancelation
4143 * of errors. The factor 0.5 is because energy distributes
4144 * equally over Ekin and Epot.
4146 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4147 if (max_T_error > T_error_warn)
4149 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.",
4150 ir->verletbuf_tol, T, tau,
4152 100*T_error_suggest,
4153 ir->verletbuf_tol*T_error_suggest/max_T_error);
4154 warning(wi, warn_buf);
4159 if (ETC_ANDERSEN(ir->etc))
4163 for (i = 0; i < ir->opts.ngtc; i++)
4165 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4166 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4167 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
4168 i, ir->opts.tau_t[i]);
4169 CHECK(ir->opts.tau_t[i] < 0);
4172 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
4174 for (i = 0; i < ir->opts.ngtc; i++)
4176 int nsteps = static_cast<int>(ir->opts.tau_t[i]/ir->delta_t + 0.5);
4177 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);
4178 CHECK(nsteps % ir->nstcomm != 0);
4183 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4184 ir->comm_mode == ecmNO &&
4185 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4186 !ETC_ANDERSEN(ir->etc))
4188 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");
4191 /* Check for pressure coupling with absolute position restraints */
4192 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4194 absolute_reference(ir, sys, TRUE, AbsRef);
4196 for (m = 0; m < DIM; m++)
4198 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4200 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4208 aloopb = gmx_mtop_atomloop_block_init(sys);
4210 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4212 if (atom->q != 0 || atom->qB != 0)
4220 if (EEL_FULL(ir->coulombtype))
4223 "You are using full electrostatics treatment %s for a system without charges.\n"
4224 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4225 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4226 warning(wi, err_buf);
4231 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4234 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4235 "You might want to consider using %s electrostatics.\n",
4237 warning_note(wi, err_buf);
4241 /* Check if combination rules used in LJ-PME are the same as in the force field */
4242 if (EVDW_PME(ir->vdwtype))
4244 check_combination_rules(ir, sys, wi);
4247 /* Generalized reaction field */
4248 if (ir->opts.ngtc == 0)
4250 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4252 CHECK(ir->coulombtype == eelGRF);
4256 sprintf(err_buf, "When using coulombtype = %s"
4257 " ref-t for temperature coupling should be > 0",
4259 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4263 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4265 for (m = 0; (m < DIM); m++)
4267 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4276 snew(mgrp, sys->groups.grps[egcACC].nr);
4277 aloop = gmx_mtop_atomloop_all_init(sys);
4279 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4281 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4284 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4286 for (m = 0; (m < DIM); m++)
4288 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4292 for (m = 0; (m < DIM); m++)
4294 if (fabs(acc[m]) > 1e-6)
4296 const char *dim[DIM] = { "X", "Y", "Z" };
4298 "Net Acceleration in %s direction, will %s be corrected\n",
4299 dim[m], ir->nstcomm != 0 ? "" : "not");
4300 if (ir->nstcomm != 0 && m < ndof_com(ir))
4303 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4305 ir->opts.acc[i][m] -= acc[m];
4313 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4314 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4316 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4324 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4326 if (ir->pull->coord[i].group[0] == 0 ||
4327 ir->pull->coord[i].group[1] == 0)
4329 absolute_reference(ir, sys, FALSE, AbsRef);
4330 for (m = 0; m < DIM; m++)
4332 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4334 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.");
4342 for (i = 0; i < 3; i++)
4344 for (m = 0; m <= i; m++)
4346 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4347 ir->deform[i][m] != 0)
4349 for (c = 0; c < ir->pull->ncoord; c++)
4351 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4352 ir->pull->coord[c].vec[m] != 0)
4354 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4365 void double_check(t_inputrec *ir, matrix box,
4366 gmx_bool bHasNormalConstraints,
4367 gmx_bool bHasAnyConstraints,
4371 char warn_buf[STRLEN];
4374 ptr = check_box(ir->ePBC, box);
4377 warning_error(wi, ptr);
4380 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4382 if (ir->shake_tol <= 0.0)
4384 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4386 warning_error(wi, warn_buf);
4390 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4392 /* If we have Lincs constraints: */
4393 if (ir->eI == eiMD && ir->etc == etcNO &&
4394 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4396 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4397 warning_note(wi, warn_buf);
4400 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4402 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4403 warning_note(wi, warn_buf);
4405 if (ir->epc == epcMTTK)
4407 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4411 if (bHasAnyConstraints && ir->epc == epcMTTK)
4413 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4416 if (ir->LincsWarnAngle > 90.0)
4418 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4419 warning(wi, warn_buf);
4420 ir->LincsWarnAngle = 90.0;
4423 if (ir->ePBC != epbcNONE)
4425 if (ir->nstlist == 0)
4427 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4429 if (ir->ns_type == ensGRID)
4431 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4433 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");
4434 warning_error(wi, warn_buf);
4439 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4440 if (2*ir->rlist >= min_size)
4442 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4443 warning_error(wi, warn_buf);
4446 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4453 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4457 real rvdw1, rvdw2, rcoul1, rcoul2;
4458 char warn_buf[STRLEN];
4460 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4464 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4469 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4475 if (rvdw1 + rvdw2 > ir->rlist ||
4476 rcoul1 + rcoul2 > ir->rlist)
4479 "The sum of the two largest charge group radii (%f) "
4480 "is larger than rlist (%f)\n",
4481 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4482 warning(wi, warn_buf);
4486 /* Here we do not use the zero at cut-off macro,
4487 * since user defined interactions might purposely
4488 * not be zero at the cut-off.
4490 if (ir_vdw_is_zero_at_cutoff(ir) &&
4491 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4493 sprintf(warn_buf, "The sum of the two largest charge group "
4494 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4495 "With exact cut-offs, better performance can be "
4496 "obtained with cutoff-scheme = %s, because it "
4497 "does not use charge groups at all.",
4499 ir->rlist, ir->rvdw,
4500 ecutscheme_names[ecutsVERLET]);
4503 warning(wi, warn_buf);
4507 warning_note(wi, warn_buf);
4510 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4511 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4513 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4514 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4516 ir->rlist, ir->rcoulomb,
4517 ecutscheme_names[ecutsVERLET]);
4520 warning(wi, warn_buf);
4524 warning_note(wi, warn_buf);