2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/fileio/readinp.h"
50 #include "gromacs/fileio/warninp.h"
51 #include "gromacs/gmxlib/chargegroup.h"
52 #include "gromacs/gmxlib/network.h"
53 #include "gromacs/gmxpreprocess/toputil.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdlib/calc_verletbuf.h"
58 #include "gromacs/mdrunutility/mdmodules.h"
59 #include "gromacs/mdtypes/inputrec.h"
60 #include "gromacs/mdtypes/md_enums.h"
61 #include "gromacs/mdtypes/pull-params.h"
62 #include "gromacs/options/options.h"
63 #include "gromacs/options/treesupport.h"
64 #include "gromacs/pbcutil/pbc.h"
65 #include "gromacs/topology/block.h"
66 #include "gromacs/topology/ifunc.h"
67 #include "gromacs/topology/index.h"
68 #include "gromacs/topology/mtop_util.h"
69 #include "gromacs/topology/symtab.h"
70 #include "gromacs/topology/topology.h"
71 #include "gromacs/utility/cstringutil.h"
72 #include "gromacs/utility/exceptions.h"
73 #include "gromacs/utility/fatalerror.h"
74 #include "gromacs/utility/filestream.h"
75 #include "gromacs/utility/gmxassert.h"
76 #include "gromacs/utility/ikeyvaluetreeerror.h"
77 #include "gromacs/utility/keyvaluetree.h"
78 #include "gromacs/utility/keyvaluetreetransform.h"
79 #include "gromacs/utility/smalloc.h"
80 #include "gromacs/utility/stringcompare.h"
81 #include "gromacs/utility/stringutil.h"
86 /* Resource parameters
87 * Do not change any of these until you read the instruction
88 * in readinp.h. Some cpp's do not take spaces after the backslash
89 * (like the c-shell), which will give you a very weird compiler
93 typedef struct t_inputrec_strings
95 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
96 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
97 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
98 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
99 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
101 char fep_lambda[efptNR][STRLEN];
102 char lambda_weights[STRLEN];
105 char anneal[STRLEN], anneal_npoints[STRLEN],
106 anneal_time[STRLEN], anneal_temp[STRLEN];
107 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
108 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
109 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
111 } gmx_inputrec_strings;
113 static gmx_inputrec_strings *is = nullptr;
115 void init_inputrec_strings()
119 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.");
124 void done_inputrec_strings()
132 egrptpALL, /* All particles have to be a member of a group. */
133 egrptpALL_GENREST, /* A rest group with name is generated for particles *
134 * that are not part of any group. */
135 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
136 * for the rest group. */
137 egrptpONE /* Merge all selected groups into one group, *
138 * make a rest group for the remaining particles. */
141 static const char *constraints[eshNR+1] = {
142 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
145 static const char *couple_lam[ecouplamNR+1] = {
146 "vdw-q", "vdw", "q", "none", nullptr
149 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
154 for (i = 0; i < ntemps; i++)
156 /* simple linear scaling -- allows more control */
157 if (simtemp->eSimTempScale == esimtempLINEAR)
159 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
161 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
163 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
165 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
167 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
172 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
173 gmx_fatal(FARGS, errorstr);
180 static void _low_check(gmx_bool b, const char *s, warninp_t wi)
184 warning_error(wi, s);
188 static void check_nst(const char *desc_nst, int nst,
189 const char *desc_p, int *p,
194 if (*p > 0 && *p % nst != 0)
196 /* Round up to the next multiple of nst */
197 *p = ((*p)/nst + 1)*nst;
198 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
199 desc_p, desc_nst, desc_p, *p);
204 static gmx_bool ir_NVE(const t_inputrec *ir)
206 return (EI_MD(ir->eI) && ir->etc == etcNO);
209 static int lcd(int n1, int n2)
214 for (i = 2; (i <= n1 && i <= n2); i++)
216 if (n1 % i == 0 && n2 % i == 0)
225 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
227 if (*eintmod == eintmodPOTSHIFT_VERLET)
229 if (ir->cutoff_scheme == ecutsVERLET)
231 *eintmod = eintmodPOTSHIFT;
235 *eintmod = eintmodNONE;
240 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
242 /* Check internal consistency.
243 * NOTE: index groups are not set here yet, don't check things
244 * like temperature coupling group options here, but in triple_check
247 /* Strange macro: first one fills the err_buf, and then one can check
248 * the condition, which will print the message and increase the error
251 #define CHECK(b) _low_check(b, err_buf, wi)
252 char err_buf[256], warn_buf[STRLEN];
255 t_lambda *fep = ir->fepvals;
256 t_expanded *expand = ir->expandedvals;
258 set_warning_line(wi, mdparin, -1);
260 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
262 sprintf(warn_buf, "%s electrostatics is no longer supported",
263 eel_names[eelRF_NEC_UNSUPPORTED]);
264 warning_error(wi, warn_buf);
267 /* BASIC CUT-OFF STUFF */
268 if (ir->rcoulomb < 0)
270 warning_error(wi, "rcoulomb should be >= 0");
274 warning_error(wi, "rvdw should be >= 0");
277 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
279 warning_error(wi, "rlist should be >= 0");
281 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.)");
282 CHECK(ir->nstlist < 0);
284 process_interaction_modifier(ir, &ir->coulomb_modifier);
285 process_interaction_modifier(ir, &ir->vdw_modifier);
287 if (ir->cutoff_scheme == ecutsGROUP)
290 "The group cutoff scheme is deprecated since GROMACS 5.0 and will be removed in a future "
291 "release when all interaction forms are supported for the verlet scheme. The verlet "
292 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
294 if (ir->rlist > 0 && ir->rlist < ir->rcoulomb)
296 gmx_fatal(FARGS, "rcoulomb must not be greater than rlist (twin-range schemes are not supported)");
298 if (ir->rlist > 0 && ir->rlist < ir->rvdw)
300 gmx_fatal(FARGS, "rvdw must not be greater than rlist (twin-range schemes are not supported)");
303 if (ir->rlist == 0 && ir->ePBC != epbcNONE)
305 warning_error(wi, "Can not have an infinite cut-off with PBC");
309 if (ir->cutoff_scheme == ecutsVERLET)
313 /* Normal Verlet type neighbor-list, currently only limited feature support */
314 if (inputrec2nboundeddim(ir) < 3)
316 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
319 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
320 if (ir->rcoulomb != ir->rvdw)
322 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
323 // for PME load balancing, we can support this exception.
324 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
325 ir->vdwtype == evdwCUT &&
326 ir->rcoulomb > ir->rvdw);
327 if (!bUsesPmeTwinRangeKernel)
329 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
333 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
335 if (ir->vdw_modifier == eintmodNONE ||
336 ir->vdw_modifier == eintmodPOTSHIFT)
338 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
340 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]);
341 warning_note(wi, warn_buf);
343 ir->vdwtype = evdwCUT;
347 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
348 warning_error(wi, warn_buf);
352 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
354 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
356 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
357 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
359 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
361 if (!(ir->coulomb_modifier == eintmodNONE ||
362 ir->coulomb_modifier == eintmodPOTSHIFT))
364 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
365 warning_error(wi, warn_buf);
368 if (ir->implicit_solvent != eisNO)
370 warning_error(wi, "Implicit solvent is not (yet) supported with the with Verlet lists.");
373 if (ir->nstlist <= 0)
375 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
378 if (ir->nstlist < 10)
380 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.");
383 rc_max = std::max(ir->rvdw, ir->rcoulomb);
385 if (ir->verletbuf_tol <= 0)
387 if (ir->verletbuf_tol == 0)
389 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
392 if (ir->rlist < rc_max)
394 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
397 if (ir->rlist == rc_max && ir->nstlist > 1)
399 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.");
404 if (ir->rlist > rc_max)
406 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.");
409 if (ir->nstlist == 1)
411 /* No buffer required */
416 if (EI_DYNAMICS(ir->eI))
418 if (inputrec2nboundeddim(ir) < 3)
420 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.");
422 /* Set rlist temporarily so we can continue processing */
427 /* Set the buffer to 5% of the cut-off */
428 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
434 /* GENERAL INTEGRATOR STUFF */
437 if (ir->etc != etcNO)
439 if (EI_RANDOM(ir->eI))
441 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]);
445 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
447 warning_note(wi, warn_buf);
451 if (ir->eI == eiVVAK)
453 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]);
454 warning_note(wi, warn_buf);
456 if (!EI_DYNAMICS(ir->eI))
458 if (ir->epc != epcNO)
460 sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
461 warning_note(wi, warn_buf);
465 if (EI_DYNAMICS(ir->eI))
467 if (ir->nstcalcenergy < 0)
469 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
470 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
472 /* nstcalcenergy larger than nstener does not make sense.
473 * We ideally want nstcalcenergy=nstener.
477 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
481 ir->nstcalcenergy = ir->nstenergy;
485 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
486 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
487 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
490 const char *nsten = "nstenergy";
491 const char *nstdh = "nstdhdl";
492 const char *min_name = nsten;
493 int min_nst = ir->nstenergy;
495 /* find the smallest of ( nstenergy, nstdhdl ) */
496 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
497 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
499 min_nst = ir->fepvals->nstdhdl;
502 /* If the user sets nstenergy small, we should respect that */
504 "Setting nstcalcenergy (%d) equal to %s (%d)",
505 ir->nstcalcenergy, min_name, min_nst);
506 warning_note(wi, warn_buf);
507 ir->nstcalcenergy = min_nst;
510 if (ir->epc != epcNO)
512 if (ir->nstpcouple < 0)
514 ir->nstpcouple = ir_optimal_nstpcouple(ir);
518 if (ir->nstcalcenergy > 0)
520 if (ir->efep != efepNO)
522 /* nstdhdl should be a multiple of nstcalcenergy */
523 check_nst("nstcalcenergy", ir->nstcalcenergy,
524 "nstdhdl", &ir->fepvals->nstdhdl, wi);
525 /* nstexpanded should be a multiple of nstcalcenergy */
526 check_nst("nstcalcenergy", ir->nstcalcenergy,
527 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
529 /* for storing exact averages nstenergy should be
530 * a multiple of nstcalcenergy
532 check_nst("nstcalcenergy", ir->nstcalcenergy,
533 "nstenergy", &ir->nstenergy, wi);
537 if (ir->nsteps == 0 && !ir->bContinuation)
539 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
543 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
544 ir->bContinuation && ir->ld_seed != -1)
546 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
552 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
553 CHECK(ir->ePBC != epbcXYZ);
554 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
555 CHECK(ir->ns_type != ensGRID);
556 sprintf(err_buf, "with TPI nstlist should be larger than zero");
557 CHECK(ir->nstlist <= 0);
558 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
559 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
560 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
561 CHECK(ir->cutoff_scheme == ecutsVERLET);
565 if ( (opts->nshake > 0) && (opts->bMorse) )
568 "Using morse bond-potentials while constraining bonds is useless");
569 warning(wi, warn_buf);
572 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
573 ir->bContinuation && ir->ld_seed != -1)
575 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)");
577 /* verify simulated tempering options */
581 gmx_bool bAllTempZero = TRUE;
582 for (i = 0; i < fep->n_lambda; i++)
584 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]);
585 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
586 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
588 bAllTempZero = FALSE;
591 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
592 CHECK(bAllTempZero == TRUE);
594 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
595 CHECK(ir->eI != eiVV);
597 /* check compatability of the temperature coupling with simulated tempering */
599 if (ir->etc == etcNOSEHOOVER)
601 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
602 warning_note(wi, warn_buf);
605 /* check that the temperatures make sense */
607 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);
608 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
610 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
611 CHECK(ir->simtempvals->simtemp_high <= 0);
613 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
614 CHECK(ir->simtempvals->simtemp_low <= 0);
617 /* verify free energy options */
619 if (ir->efep != efepNO)
622 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
624 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
626 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
627 (int)fep->sc_r_power);
628 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
630 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
631 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
633 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
634 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
636 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
637 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
639 sprintf(err_buf, "Free-energy not implemented for Ewald");
640 CHECK(ir->coulombtype == eelEWALD);
642 /* check validty of lambda inputs */
643 if (fep->n_lambda == 0)
645 /* Clear output in case of no states:*/
646 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
647 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
651 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
652 CHECK((fep->init_fep_state >= fep->n_lambda));
655 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
656 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
658 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",
659 fep->init_lambda, fep->init_fep_state);
660 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
664 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
668 for (i = 0; i < efptNR; i++)
670 if (fep->separate_dvdl[i])
675 if (n_lambda_terms > 1)
677 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.");
678 warning(wi, warn_buf);
681 if (n_lambda_terms < 2 && fep->n_lambda > 0)
684 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
688 for (j = 0; j < efptNR; j++)
690 for (i = 0; i < fep->n_lambda; i++)
692 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]);
693 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
697 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
699 for (i = 0; i < fep->n_lambda; i++)
701 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],
702 fep->all_lambda[efptCOUL][i]);
703 CHECK((fep->sc_alpha > 0) &&
704 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
705 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
706 ((fep->all_lambda[efptVDW][i] > 0.0) &&
707 (fep->all_lambda[efptVDW][i] < 1.0))));
711 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
713 real sigma, lambda, r_sc;
716 /* Maximum estimate for A and B charges equal with lambda power 1 */
718 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);
719 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.",
721 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
722 warning_note(wi, warn_buf);
725 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
726 be treated differently, but that's the next step */
728 for (i = 0; i < efptNR; i++)
730 for (j = 0; j < fep->n_lambda; j++)
732 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
733 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
738 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
742 /* checking equilibration of weights inputs for validity */
744 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
745 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
746 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
748 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
749 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
750 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
752 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
753 expand->equil_steps, elmceq_names[elmceqSTEPS]);
754 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
756 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
757 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
758 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
760 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
761 expand->equil_ratio, elmceq_names[elmceqRATIO]);
762 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
764 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
765 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
766 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
768 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
769 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
770 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
772 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
773 expand->equil_steps, elmceq_names[elmceqSTEPS]);
774 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
776 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
777 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
778 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
780 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
781 expand->equil_ratio, elmceq_names[elmceqRATIO]);
782 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
784 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
785 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
786 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
788 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
789 CHECK((expand->lmc_repeats <= 0));
790 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
791 CHECK((expand->minvarmin <= 0));
792 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
793 CHECK((expand->c_range < 0));
794 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
795 fep->init_fep_state, expand->lmc_forced_nstart);
796 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
797 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
798 CHECK((expand->lmc_forced_nstart < 0));
799 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
800 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
802 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
803 CHECK((expand->init_wl_delta < 0));
804 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
805 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
806 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
807 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
809 /* if there is no temperature control, we need to specify an MC temperature */
810 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
811 if (expand->nstTij > 0)
813 sprintf(err_buf, "nstlog must be non-zero");
814 CHECK(ir->nstlog == 0);
815 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
816 expand->nstTij, ir->nstlog);
817 CHECK((expand->nstTij % ir->nstlog) != 0);
822 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
823 CHECK(ir->nwall && ir->ePBC != epbcXY);
826 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
828 if (ir->ePBC == epbcNONE)
830 if (ir->epc != epcNO)
832 warning(wi, "Turning off pressure coupling for vacuum system");
838 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
839 epbc_names[ir->ePBC]);
840 CHECK(ir->epc != epcNO);
842 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
843 CHECK(EEL_FULL(ir->coulombtype));
845 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
846 epbc_names[ir->ePBC]);
847 CHECK(ir->eDispCorr != edispcNO);
850 if (ir->rlist == 0.0)
852 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
853 "with coulombtype = %s or coulombtype = %s\n"
854 "without periodic boundary conditions (pbc = %s) and\n"
855 "rcoulomb and rvdw set to zero",
856 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
857 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
858 (ir->ePBC != epbcNONE) ||
859 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
863 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
868 if (ir->nstcomm == 0)
870 ir->comm_mode = ecmNO;
872 if (ir->comm_mode != ecmNO)
876 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");
877 ir->nstcomm = abs(ir->nstcomm);
880 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
882 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
883 ir->nstcomm = ir->nstcalcenergy;
886 if (ir->comm_mode == ecmANGULAR)
888 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
889 CHECK(ir->bPeriodicMols);
890 if (ir->ePBC != epbcNONE)
892 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.");
897 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
899 warning_note(wi, "Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
902 /* TEMPERATURE COUPLING */
903 if (ir->etc == etcYES)
905 ir->etc = etcBERENDSEN;
906 warning_note(wi, "Old option for temperature coupling given: "
907 "changing \"yes\" to \"Berendsen\"\n");
910 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
912 if (ir->opts.nhchainlength < 1)
914 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
915 ir->opts.nhchainlength = 1;
916 warning(wi, warn_buf);
919 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
921 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
922 ir->opts.nhchainlength = 1;
927 ir->opts.nhchainlength = 0;
930 if (ir->eI == eiVVAK)
932 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
934 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
937 if (ETC_ANDERSEN(ir->etc))
939 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
940 CHECK(!(EI_VV(ir->eI)));
942 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
944 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]);
945 warning_note(wi, warn_buf);
948 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]);
949 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
952 if (ir->etc == etcBERENDSEN)
954 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
955 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
956 warning_note(wi, warn_buf);
959 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
960 && ir->epc == epcBERENDSEN)
962 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
963 "true ensemble for the thermostat");
964 warning(wi, warn_buf);
967 /* PRESSURE COUPLING */
968 if (ir->epc == epcISOTROPIC)
970 ir->epc = epcBERENDSEN;
971 warning_note(wi, "Old option for pressure coupling given: "
972 "changing \"Isotropic\" to \"Berendsen\"\n");
975 if (ir->epc != epcNO)
977 dt_pcoupl = ir->nstpcouple*ir->delta_t;
979 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
980 CHECK(ir->tau_p <= 0);
982 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
984 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
985 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
986 warning(wi, warn_buf);
989 sprintf(err_buf, "compressibility must be > 0 when using pressure"
990 " coupling %s\n", EPCOUPLTYPE(ir->epc));
991 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
992 ir->compress[ZZ][ZZ] < 0 ||
993 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
994 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
996 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
999 "You are generating velocities so I am assuming you "
1000 "are equilibrating a system. You are using "
1001 "%s pressure coupling, but this can be "
1002 "unstable for equilibration. If your system crashes, try "
1003 "equilibrating first with Berendsen pressure coupling. If "
1004 "you are not equilibrating the system, you can probably "
1005 "ignore this warning.",
1006 epcoupl_names[ir->epc]);
1007 warning(wi, warn_buf);
1013 if (ir->epc > epcNO)
1015 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1017 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.");
1023 if (ir->epc == epcMTTK)
1025 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1029 /* ELECTROSTATICS */
1030 /* More checks are in triple check (grompp.c) */
1032 if (ir->coulombtype == eelSWITCH)
1034 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1035 "artifacts, advice: use coulombtype = %s",
1036 eel_names[ir->coulombtype],
1037 eel_names[eelRF_ZERO]);
1038 warning(wi, warn_buf);
1041 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1043 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1044 warning_note(wi, warn_buf);
1047 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1049 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);
1050 warning(wi, warn_buf);
1051 ir->epsilon_rf = ir->epsilon_r;
1052 ir->epsilon_r = 1.0;
1055 if (ir->epsilon_r == 0)
1058 "It is pointless to use long-range or Generalized Born electrostatics with infinite relative permittivity."
1059 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1060 CHECK(EEL_FULL(ir->coulombtype) || ir->implicit_solvent == eisGBSA);
1063 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1065 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1066 CHECK(ir->epsilon_r < 0);
1069 if (EEL_RF(ir->coulombtype))
1071 /* reaction field (at the cut-off) */
1073 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1075 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1076 eel_names[ir->coulombtype]);
1077 warning(wi, warn_buf);
1078 ir->epsilon_rf = 0.0;
1081 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1082 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1083 (ir->epsilon_r == 0));
1084 if (ir->epsilon_rf == ir->epsilon_r)
1086 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1087 eel_names[ir->coulombtype]);
1088 warning(wi, warn_buf);
1091 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1092 * means the interaction is zero outside rcoulomb, but it helps to
1093 * provide accurate energy conservation.
1095 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1097 if (ir_coulomb_switched(ir))
1100 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1101 eel_names[ir->coulombtype]);
1102 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1105 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1107 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1109 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1110 eel_names[ir->coulombtype]);
1111 CHECK(ir->rlist > ir->rcoulomb);
1115 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1118 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1119 CHECK( ir->coulomb_modifier != eintmodNONE);
1121 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1124 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1125 CHECK( ir->vdw_modifier != eintmodNONE);
1128 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1129 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1132 "The switch/shift interaction settings are just for compatibility; you will get better "
1133 "performance from applying potential modifiers to your interactions!\n");
1134 warning_note(wi, warn_buf);
1137 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1139 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1141 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1142 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.",
1143 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1144 warning(wi, warn_buf);
1148 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1150 if (ir->rvdw_switch == 0)
1152 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.");
1153 warning(wi, warn_buf);
1157 if (EEL_FULL(ir->coulombtype))
1159 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1160 ir->coulombtype == eelPMEUSERSWITCH)
1162 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1163 eel_names[ir->coulombtype]);
1164 CHECK(ir->rcoulomb > ir->rlist);
1166 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1168 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1171 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist.\n"
1172 "For optimal energy conservation,consider using\n"
1173 "a potential modifier.", eel_names[ir->coulombtype]);
1174 CHECK(ir->rcoulomb != ir->rlist);
1179 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1181 // TODO: Move these checks into the ewald module with the options class
1183 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1185 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1187 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
1188 warning_error(wi, warn_buf);
1192 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1194 if (ir->ewald_geometry == eewg3D)
1196 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1197 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1198 warning(wi, warn_buf);
1200 /* This check avoids extra pbc coding for exclusion corrections */
1201 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1202 CHECK(ir->wall_ewald_zfac < 2);
1204 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1205 EEL_FULL(ir->coulombtype))
1207 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1208 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1209 warning(wi, warn_buf);
1211 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1213 if (ir->cutoff_scheme == ecutsVERLET)
1215 sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
1216 eel_names[ir->coulombtype]);
1217 warning(wi, warn_buf);
1221 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",
1222 eel_names[ir->coulombtype]);
1223 warning_note(wi, warn_buf);
1227 if (ir_vdw_switched(ir))
1229 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1230 CHECK(ir->rvdw_switch >= ir->rvdw);
1232 if (ir->rvdw_switch < 0.5*ir->rvdw)
1234 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.",
1235 ir->rvdw_switch, ir->rvdw);
1236 warning_note(wi, warn_buf);
1239 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1241 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1243 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1244 CHECK(ir->rlist > ir->rvdw);
1248 if (ir->vdwtype == evdwPME)
1250 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1252 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1253 evdw_names[ir->vdwtype],
1254 eintmod_names[eintmodPOTSHIFT],
1255 eintmod_names[eintmodNONE]);
1259 if (ir->cutoff_scheme == ecutsGROUP)
1261 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1262 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)))
1264 warning_note(wi, "With exact cut-offs, rlist should be "
1265 "larger than rcoulomb and rvdw, so that there "
1266 "is a buffer region for particle motion "
1267 "between neighborsearch steps");
1270 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlist <= ir->rcoulomb)
1272 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rcoulomb.");
1273 warning_note(wi, warn_buf);
1275 if (ir_vdw_switched(ir) && (ir->rlist <= ir->rvdw))
1277 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rvdw.");
1278 warning_note(wi, warn_buf);
1282 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1284 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.");
1287 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1290 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1293 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1295 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1298 /* ENERGY CONSERVATION */
1299 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1301 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1303 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1304 evdw_names[evdwSHIFT]);
1305 warning_note(wi, warn_buf);
1307 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1309 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1310 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1311 warning_note(wi, warn_buf);
1315 /* IMPLICIT SOLVENT */
1316 if (ir->coulombtype == eelGB_NOTUSED)
1318 sprintf(warn_buf, "Invalid option %s for coulombtype",
1319 eel_names[ir->coulombtype]);
1320 warning_error(wi, warn_buf);
1323 if (ir->sa_algorithm == esaSTILL)
1325 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1326 CHECK(ir->sa_algorithm == esaSTILL);
1329 if (ir->implicit_solvent == eisGBSA)
1331 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1332 CHECK(ir->rgbradii != ir->rlist);
1334 if (ir->coulombtype != eelCUT)
1336 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1337 CHECK(ir->coulombtype != eelCUT);
1339 if (ir->vdwtype != evdwCUT)
1341 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1342 CHECK(ir->vdwtype != evdwCUT);
1344 if (ir->nstgbradii < 1)
1346 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1347 warning_note(wi, warn_buf);
1350 if (ir->sa_algorithm == esaNO)
1352 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1353 warning_note(wi, warn_buf);
1355 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1357 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");
1358 warning_note(wi, warn_buf);
1360 if (ir->gb_algorithm == egbSTILL)
1362 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1366 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1369 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1371 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1372 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1379 if (ir->cutoff_scheme != ecutsGROUP)
1381 warning_error(wi, "QMMM is currently only supported with cutoff-scheme=group");
1383 if (!EI_DYNAMICS(ir->eI))
1386 sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
1387 warning_error(wi, buf);
1393 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1397 /* count the number of text elemets separated by whitespace in a string.
1398 str = the input string
1399 maxptr = the maximum number of allowed elements
1400 ptr = the output array of pointers to the first character of each element
1401 returns: the number of elements. */
1402 int str_nelem(const char *str, int maxptr, char *ptr[])
1407 copy0 = gmx_strdup(str);
1410 while (*copy != '\0')
1414 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1422 while ((*copy != '\0') && !isspace(*copy))
1441 /* interpret a number of doubles from a string and put them in an array,
1442 after allocating space for them.
1443 str = the input string
1444 n = the (pre-allocated) number of doubles read
1445 r = the output array of doubles. */
1446 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1451 char warn_buf[STRLEN];
1453 *n = str_nelem(str, MAXPTR, ptr);
1456 for (i = 0; i < *n; i++)
1458 (*r)[i] = strtod(ptr[i], &endptr);
1461 sprintf(warn_buf, "Invalid value %s in string in mdp file. Expected a real number.", ptr[i]);
1462 warning_error(wi, warn_buf);
1467 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1470 int i, j, max_n_lambda, nweights, nfep[efptNR];
1471 t_lambda *fep = ir->fepvals;
1472 t_expanded *expand = ir->expandedvals;
1473 real **count_fep_lambdas;
1474 gmx_bool bOneLambda = TRUE;
1476 snew(count_fep_lambdas, efptNR);
1478 /* FEP input processing */
1479 /* first, identify the number of lambda values for each type.
1480 All that are nonzero must have the same number */
1482 for (i = 0; i < efptNR; i++)
1484 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1487 /* now, determine the number of components. All must be either zero, or equal. */
1490 for (i = 0; i < efptNR; i++)
1492 if (nfep[i] > max_n_lambda)
1494 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1495 must have the same number if its not zero.*/
1500 for (i = 0; i < efptNR; i++)
1504 ir->fepvals->separate_dvdl[i] = FALSE;
1506 else if (nfep[i] == max_n_lambda)
1508 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1509 respect to the temperature currently */
1511 ir->fepvals->separate_dvdl[i] = TRUE;
1516 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1517 nfep[i], efpt_names[i], max_n_lambda);
1520 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1521 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1523 /* the number of lambdas is the number we've read in, which is either zero
1524 or the same for all */
1525 fep->n_lambda = max_n_lambda;
1527 /* allocate space for the array of lambda values */
1528 snew(fep->all_lambda, efptNR);
1529 /* if init_lambda is defined, we need to set lambda */
1530 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1532 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1534 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1535 for (i = 0; i < efptNR; i++)
1537 snew(fep->all_lambda[i], fep->n_lambda);
1538 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1541 for (j = 0; j < fep->n_lambda; j++)
1543 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1545 sfree(count_fep_lambdas[i]);
1548 sfree(count_fep_lambdas);
1550 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1551 bookkeeping -- for now, init_lambda */
1553 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1555 for (i = 0; i < fep->n_lambda; i++)
1557 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1561 /* check to see if only a single component lambda is defined, and soft core is defined.
1562 In this case, turn on coulomb soft core */
1564 if (max_n_lambda == 0)
1570 for (i = 0; i < efptNR; i++)
1572 if ((nfep[i] != 0) && (i != efptFEP))
1578 if ((bOneLambda) && (fep->sc_alpha > 0))
1580 fep->bScCoul = TRUE;
1583 /* Fill in the others with the efptFEP if they are not explicitly
1584 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1585 they are all zero. */
1587 for (i = 0; i < efptNR; i++)
1589 if ((nfep[i] == 0) && (i != efptFEP))
1591 for (j = 0; j < fep->n_lambda; j++)
1593 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1599 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1600 if (fep->sc_r_power == 48)
1602 if (fep->sc_alpha > 0.1)
1604 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1608 /* now read in the weights */
1609 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1612 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1614 else if (nweights != fep->n_lambda)
1616 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1617 nweights, fep->n_lambda);
1619 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1621 expand->nstexpanded = fep->nstdhdl;
1622 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1624 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1626 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1627 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1628 2*tau_t just to be careful so it's not to frequent */
1633 static void do_simtemp_params(t_inputrec *ir)
1636 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1637 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1642 static void do_wall_params(t_inputrec *ir,
1643 char *wall_atomtype, char *wall_density,
1647 char *names[MAXPTR];
1650 opts->wall_atomtype[0] = nullptr;
1651 opts->wall_atomtype[1] = nullptr;
1653 ir->wall_atomtype[0] = -1;
1654 ir->wall_atomtype[1] = -1;
1655 ir->wall_density[0] = 0;
1656 ir->wall_density[1] = 0;
1660 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1661 if (nstr != ir->nwall)
1663 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1666 for (i = 0; i < ir->nwall; i++)
1668 opts->wall_atomtype[i] = gmx_strdup(names[i]);
1671 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1673 nstr = str_nelem(wall_density, MAXPTR, names);
1674 if (nstr != ir->nwall)
1676 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1678 for (i = 0; i < ir->nwall; i++)
1680 if (sscanf(names[i], "%lf", &dbl) != 1)
1682 gmx_fatal(FARGS, "Could not parse wall-density value from string '%s'", names[i]);
1686 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1688 ir->wall_density[i] = dbl;
1694 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1702 srenew(groups->grpname, groups->ngrpname+nwall);
1703 grps = &(groups->grps[egcENER]);
1704 srenew(grps->nm_ind, grps->nr+nwall);
1705 for (i = 0; i < nwall; i++)
1707 sprintf(str, "wall%d", i);
1708 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1709 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1714 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1715 t_expanded *expand, warninp_t wi)
1723 /* read expanded ensemble parameters */
1724 CCTYPE ("expanded ensemble variables");
1725 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1726 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1727 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1728 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1729 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1730 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1731 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1732 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1733 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1734 CCTYPE("Seed for Monte Carlo in lambda space");
1735 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1736 RTYPE ("mc-temperature", expand->mc_temp, -1);
1737 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1738 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1739 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1740 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1741 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1742 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1743 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1744 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1745 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1746 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1747 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1755 /*! \brief Return whether an end state with the given coupling-lambda
1756 * value describes fully-interacting VDW.
1758 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1759 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1761 static gmx_bool couple_lambda_has_vdw_on(int couple_lambda_value)
1763 return (couple_lambda_value == ecouplamVDW ||
1764 couple_lambda_value == ecouplamVDWQ);
1770 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1773 explicit MdpErrorHandler(warninp_t wi)
1774 : wi_(wi), mapping_(nullptr)
1778 void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
1780 mapping_ = &mapping;
1783 virtual bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context)
1785 ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
1786 getOptionName(context).c_str()));
1787 std::string message = gmx::formatExceptionMessageToString(*ex);
1788 warning_error(wi_, message.c_str());
1793 std::string getOptionName(const gmx::KeyValueTreePath &context)
1795 if (mapping_ != nullptr)
1797 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1798 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1801 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1806 const gmx::IKeyValueTreeBackMapping *mapping_;
1811 void get_ir(const char *mdparin, const char *mdparout,
1812 gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
1813 WriteMdpHeader writeMdpHeader, warninp_t wi)
1816 double dumdub[2][6];
1820 char warn_buf[STRLEN];
1821 t_lambda *fep = ir->fepvals;
1822 t_expanded *expand = ir->expandedvals;
1824 init_inputrec_strings();
1825 gmx::TextInputFile stream(mdparin);
1826 inp = read_inpfile(&stream, mdparin, &ninp, wi);
1828 snew(dumstr[0], STRLEN);
1829 snew(dumstr[1], STRLEN);
1831 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1834 "%s did not specify a value for the .mdp option "
1835 "\"cutoff-scheme\". Probably it was first intended for use "
1836 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1837 "introduced, but the group scheme was still the default. "
1838 "The default is now the Verlet scheme, so you will observe "
1839 "different behaviour.", mdparin);
1840 warning_note(wi, warn_buf);
1843 /* ignore the following deprecated commands */
1846 REM_TYPE("domain-decomposition");
1847 REM_TYPE("andersen-seed");
1849 REM_TYPE("dihre-fc");
1850 REM_TYPE("dihre-tau");
1851 REM_TYPE("nstdihreout");
1852 REM_TYPE("nstcheckpoint");
1853 REM_TYPE("optimize-fft");
1854 REM_TYPE("adress_type");
1855 REM_TYPE("adress_const_wf");
1856 REM_TYPE("adress_ex_width");
1857 REM_TYPE("adress_hy_width");
1858 REM_TYPE("adress_ex_forcecap");
1859 REM_TYPE("adress_interface_correction");
1860 REM_TYPE("adress_site");
1861 REM_TYPE("adress_reference_coords");
1862 REM_TYPE("adress_tf_grp_names");
1863 REM_TYPE("adress_cg_grp_names");
1864 REM_TYPE("adress_do_hybridpairs");
1865 REM_TYPE("rlistlong");
1866 REM_TYPE("nstcalclr");
1867 REM_TYPE("pull-print-com2");
1869 /* replace the following commands with the clearer new versions*/
1870 REPL_TYPE("unconstrained-start", "continuation");
1871 REPL_TYPE("foreign-lambda", "fep-lambdas");
1872 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1873 REPL_TYPE("nstxtcout", "nstxout-compressed");
1874 REPL_TYPE("xtc-grps", "compressed-x-grps");
1875 REPL_TYPE("xtc-precision", "compressed-x-precision");
1876 REPL_TYPE("pull-print-com1", "pull-print-com");
1878 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1879 CTYPE ("Preprocessor information: use cpp syntax.");
1880 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1881 STYPE ("include", opts->include, nullptr);
1882 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1883 STYPE ("define", opts->define, nullptr);
1885 CCTYPE ("RUN CONTROL PARAMETERS");
1886 EETYPE("integrator", ir->eI, ei_names);
1887 CTYPE ("Start time and timestep in ps");
1888 RTYPE ("tinit", ir->init_t, 0.0);
1889 RTYPE ("dt", ir->delta_t, 0.001);
1890 STEPTYPE ("nsteps", ir->nsteps, 0);
1891 CTYPE ("For exact run continuation or redoing part of a run");
1892 STEPTYPE ("init-step", ir->init_step, 0);
1893 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1894 ITYPE ("simulation-part", ir->simulation_part, 1);
1895 CTYPE ("mode for center of mass motion removal");
1896 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1897 CTYPE ("number of steps for center of mass motion removal");
1898 ITYPE ("nstcomm", ir->nstcomm, 100);
1899 CTYPE ("group(s) for center of mass motion removal");
1900 STYPE ("comm-grps", is->vcm, nullptr);
1902 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1903 CTYPE ("Friction coefficient (amu/ps) and random seed");
1904 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1905 STEPTYPE ("ld-seed", ir->ld_seed, -1);
1908 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1909 CTYPE ("Force tolerance and initial step-size");
1910 RTYPE ("emtol", ir->em_tol, 10.0);
1911 RTYPE ("emstep", ir->em_stepsize, 0.01);
1912 CTYPE ("Max number of iterations in relax-shells");
1913 ITYPE ("niter", ir->niter, 20);
1914 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1915 RTYPE ("fcstep", ir->fc_stepsize, 0);
1916 CTYPE ("Frequency of steepest descents steps when doing CG");
1917 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1918 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1920 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1921 RTYPE ("rtpi", ir->rtpi, 0.05);
1923 /* Output options */
1924 CCTYPE ("OUTPUT CONTROL OPTIONS");
1925 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1926 ITYPE ("nstxout", ir->nstxout, 0);
1927 ITYPE ("nstvout", ir->nstvout, 0);
1928 ITYPE ("nstfout", ir->nstfout, 0);
1929 CTYPE ("Output frequency for energies to log file and energy file");
1930 ITYPE ("nstlog", ir->nstlog, 1000);
1931 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1932 ITYPE ("nstenergy", ir->nstenergy, 1000);
1933 CTYPE ("Output frequency and precision for .xtc file");
1934 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
1935 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1936 CTYPE ("This selects the subset of atoms for the compressed");
1937 CTYPE ("trajectory file. You can select multiple groups. By");
1938 CTYPE ("default, all atoms will be written.");
1939 STYPE ("compressed-x-grps", is->x_compressed_groups, nullptr);
1940 CTYPE ("Selection of energy groups");
1941 STYPE ("energygrps", is->energy, nullptr);
1943 /* Neighbor searching */
1944 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1945 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1946 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1947 CTYPE ("nblist update frequency");
1948 ITYPE ("nstlist", ir->nstlist, 10);
1949 CTYPE ("ns algorithm (simple or grid)");
1950 EETYPE("ns-type", ir->ns_type, ens_names);
1951 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1952 EETYPE("pbc", ir->ePBC, epbc_names);
1953 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1954 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1955 CTYPE ("a value of -1 means: use rlist");
1956 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1957 CTYPE ("nblist cut-off");
1958 RTYPE ("rlist", ir->rlist, 1.0);
1959 CTYPE ("long-range cut-off for switched potentials");
1961 /* Electrostatics */
1962 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1963 CTYPE ("Method for doing electrostatics");
1964 EETYPE("coulombtype", ir->coulombtype, eel_names);
1965 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1966 CTYPE ("cut-off lengths");
1967 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1968 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1969 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1970 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1971 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1972 CTYPE ("Method for doing Van der Waals");
1973 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1974 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1975 CTYPE ("cut-off lengths");
1976 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1977 RTYPE ("rvdw", ir->rvdw, 1.0);
1978 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1979 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1980 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1981 RTYPE ("table-extension", ir->tabext, 1.0);
1982 CTYPE ("Separate tables between energy group pairs");
1983 STYPE ("energygrp-table", is->egptable, nullptr);
1984 CTYPE ("Spacing for the PME/PPPM FFT grid");
1985 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1986 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1987 ITYPE ("fourier-nx", ir->nkx, 0);
1988 ITYPE ("fourier-ny", ir->nky, 0);
1989 ITYPE ("fourier-nz", ir->nkz, 0);
1990 CTYPE ("EWALD/PME/PPPM parameters");
1991 ITYPE ("pme-order", ir->pme_order, 4);
1992 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1993 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1994 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1995 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1996 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1998 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1999 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
2001 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
2002 CTYPE ("Algorithm for calculating Born radii");
2003 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
2004 CTYPE ("Frequency of calculating the Born radii inside rlist");
2005 ITYPE ("nstgbradii", ir->nstgbradii, 1);
2006 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
2007 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
2008 RTYPE ("rgbradii", ir->rgbradii, 1.0);
2009 CTYPE ("Dielectric coefficient of the implicit solvent");
2010 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
2011 CTYPE ("Salt concentration in M for Generalized Born models");
2012 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
2013 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
2014 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
2015 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
2016 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
2017 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
2018 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
2019 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
2020 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
2021 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
2023 /* Coupling stuff */
2024 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
2025 CTYPE ("Temperature coupling");
2026 EETYPE("tcoupl", ir->etc, etcoupl_names);
2027 ITYPE ("nsttcouple", ir->nsttcouple, -1);
2028 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
2029 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
2030 CTYPE ("Groups to couple separately");
2031 STYPE ("tc-grps", is->tcgrps, nullptr);
2032 CTYPE ("Time constant (ps) and reference temperature (K)");
2033 STYPE ("tau-t", is->tau_t, nullptr);
2034 STYPE ("ref-t", is->ref_t, nullptr);
2035 CTYPE ("pressure coupling");
2036 EETYPE("pcoupl", ir->epc, epcoupl_names);
2037 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
2038 ITYPE ("nstpcouple", ir->nstpcouple, -1);
2039 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
2040 RTYPE ("tau-p", ir->tau_p, 1.0);
2041 STYPE ("compressibility", dumstr[0], nullptr);
2042 STYPE ("ref-p", dumstr[1], nullptr);
2043 CTYPE ("Scaling of reference coordinates, No, All or COM");
2044 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
2047 CCTYPE ("OPTIONS FOR QMMM calculations");
2048 EETYPE("QMMM", ir->bQMMM, yesno_names);
2049 CTYPE ("Groups treated Quantum Mechanically");
2050 STYPE ("QMMM-grps", is->QMMM, nullptr);
2051 CTYPE ("QM method");
2052 STYPE("QMmethod", is->QMmethod, nullptr);
2053 CTYPE ("QMMM scheme");
2054 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
2055 CTYPE ("QM basisset");
2056 STYPE("QMbasis", is->QMbasis, nullptr);
2057 CTYPE ("QM charge");
2058 STYPE ("QMcharge", is->QMcharge, nullptr);
2059 CTYPE ("QM multiplicity");
2060 STYPE ("QMmult", is->QMmult, nullptr);
2061 CTYPE ("Surface Hopping");
2062 STYPE ("SH", is->bSH, nullptr);
2063 CTYPE ("CAS space options");
2064 STYPE ("CASorbitals", is->CASorbitals, nullptr);
2065 STYPE ("CASelectrons", is->CASelectrons, nullptr);
2066 STYPE ("SAon", is->SAon, nullptr);
2067 STYPE ("SAoff", is->SAoff, nullptr);
2068 STYPE ("SAsteps", is->SAsteps, nullptr);
2069 CTYPE ("Scale factor for MM charges");
2070 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
2071 CTYPE ("Optimization of QM subsystem");
2072 STYPE ("bOPT", is->bOPT, nullptr);
2073 STYPE ("bTS", is->bTS, nullptr);
2075 /* Simulated annealing */
2076 CCTYPE("SIMULATED ANNEALING");
2077 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
2078 STYPE ("annealing", is->anneal, nullptr);
2079 CTYPE ("Number of time points to use for specifying annealing in each group");
2080 STYPE ("annealing-npoints", is->anneal_npoints, nullptr);
2081 CTYPE ("List of times at the annealing points for each group");
2082 STYPE ("annealing-time", is->anneal_time, nullptr);
2083 CTYPE ("Temp. at each annealing point, for each group.");
2084 STYPE ("annealing-temp", is->anneal_temp, nullptr);
2087 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
2088 EETYPE("gen-vel", opts->bGenVel, yesno_names);
2089 RTYPE ("gen-temp", opts->tempi, 300.0);
2090 ITYPE ("gen-seed", opts->seed, -1);
2093 CCTYPE ("OPTIONS FOR BONDS");
2094 EETYPE("constraints", opts->nshake, constraints);
2095 CTYPE ("Type of constraint algorithm");
2096 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
2097 CTYPE ("Do not constrain the start configuration");
2098 EETYPE("continuation", ir->bContinuation, yesno_names);
2099 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2100 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2101 CTYPE ("Relative tolerance of shake");
2102 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2103 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2104 ITYPE ("lincs-order", ir->nProjOrder, 4);
2105 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2106 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2107 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2108 ITYPE ("lincs-iter", ir->nLincsIter, 1);
2109 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2110 CTYPE ("rotates over more degrees than");
2111 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2112 CTYPE ("Convert harmonic bonds to morse potentials");
2113 EETYPE("morse", opts->bMorse, yesno_names);
2115 /* Energy group exclusions */
2116 CCTYPE ("ENERGY GROUP EXCLUSIONS");
2117 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2118 STYPE ("energygrp-excl", is->egpexcl, nullptr);
2122 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2123 ITYPE ("nwall", ir->nwall, 0);
2124 EETYPE("wall-type", ir->wall_type, ewt_names);
2125 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2126 STYPE ("wall-atomtype", is->wall_atomtype, nullptr);
2127 STYPE ("wall-density", is->wall_density, nullptr);
2128 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2131 CCTYPE("COM PULLING");
2132 EETYPE("pull", ir->bPull, yesno_names);
2136 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, wi);
2139 /* Enforced rotation */
2140 CCTYPE("ENFORCED ROTATION");
2141 CTYPE("Enforced rotation: No or Yes");
2142 EETYPE("rotation", ir->bRot, yesno_names);
2146 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2149 /* Interactive MD */
2151 CCTYPE("Group to display and/or manipulate in interactive MD session");
2152 STYPE ("IMD-group", is->imd_grp, nullptr);
2153 if (is->imd_grp[0] != '\0')
2160 CCTYPE("NMR refinement stuff");
2161 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2162 EETYPE("disre", ir->eDisre, edisre_names);
2163 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2164 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2165 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2166 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2167 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2168 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2169 CTYPE ("Output frequency for pair distances to energy file");
2170 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2171 CTYPE ("Orientation restraints: No or Yes");
2172 EETYPE("orire", opts->bOrire, yesno_names);
2173 CTYPE ("Orientation restraints force constant and tau for time averaging");
2174 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2175 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2176 STYPE ("orire-fitgrp", is->orirefitgrp, nullptr);
2177 CTYPE ("Output frequency for trace(SD) and S to energy file");
2178 ITYPE ("nstorireout", ir->nstorireout, 100);
2180 /* free energy variables */
2181 CCTYPE ("Free energy variables");
2182 EETYPE("free-energy", ir->efep, efep_names);
2183 STYPE ("couple-moltype", is->couple_moltype, nullptr);
2184 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2185 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2186 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2188 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2190 it was not entered */
2191 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2192 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2193 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2194 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], nullptr);
2195 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], nullptr);
2196 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
2197 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
2198 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
2199 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
2200 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
2201 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2202 STYPE ("init-lambda-weights", is->lambda_weights, nullptr);
2203 EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
2204 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2205 ITYPE ("sc-power", fep->sc_power, 1);
2206 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2207 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2208 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2209 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2210 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2211 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2212 separate_dhdl_file_names);
2213 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2214 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2215 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2217 /* Non-equilibrium MD stuff */
2218 CCTYPE("Non-equilibrium MD stuff");
2219 STYPE ("acc-grps", is->accgrps, nullptr);
2220 STYPE ("accelerate", is->acc, nullptr);
2221 STYPE ("freezegrps", is->freeze, nullptr);
2222 STYPE ("freezedim", is->frdim, nullptr);
2223 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2224 STYPE ("deform", is->deform, nullptr);
2226 /* simulated tempering variables */
2227 CCTYPE("simulated tempering variables");
2228 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2229 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2230 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2231 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2233 /* expanded ensemble variables */
2234 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2236 read_expandedparams(&ninp, &inp, expand, wi);
2239 /* Electric fields */
2241 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(ninp, inp);
2242 gmx::KeyValueTreeTransformer transform;
2243 transform.rules()->addRule()
2244 .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2245 mdModules->initMdpTransform(transform.rules());
2246 for (const auto &path : transform.mappedPaths())
2248 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2249 mark_einp_set(ninp, inp, path[0].c_str());
2251 MdpErrorHandler errorHandler(wi);
2253 = transform.transform(convertedValues, &errorHandler);
2254 ir->params = new gmx::KeyValueTreeObject(result.object());
2255 mdModules->adjustInputrecBasedOnModules(ir);
2256 errorHandler.setBackMapping(result.backMapping());
2257 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2260 /* Ion/water position swapping ("computational electrophysiology") */
2261 CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2262 CTYPE("Swap positions along direction: no, X, Y, Z");
2263 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2264 if (ir->eSwapCoords != eswapNO)
2271 CTYPE("Swap attempt frequency");
2272 ITYPE("swap-frequency", ir->swap->nstswap, 1);
2273 CTYPE("Number of ion types to be controlled");
2274 ITYPE("iontypes", nIonTypes, 1);
2277 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2279 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2280 snew(ir->swap->grp, ir->swap->ngrp);
2281 for (i = 0; i < ir->swap->ngrp; i++)
2283 snew(ir->swap->grp[i].molname, STRLEN);
2285 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2286 STYPE("split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2287 STYPE("split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2288 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2289 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2290 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2292 CTYPE("Name of solvent molecules");
2293 STYPE("solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2295 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2296 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2297 CTYPE("however, if correctly defined, the permeation events are recorded per channel");
2298 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2299 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2300 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2301 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2302 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2303 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2305 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2306 ITYPE("coupl-steps", ir->swap->nAverage, 10);
2308 CTYPE("Names of the ion types that can be exchanged with solvent molecules,");
2309 CTYPE("and the requested number of ions of this type in compartments A and B");
2310 CTYPE("-1 means fix the numbers as found in step 0");
2311 for (i = 0; i < nIonTypes; i++)
2313 int ig = eSwapFixedGrpNR + i;
2315 sprintf(buf, "iontype%d-name", i);
2316 STYPE(buf, ir->swap->grp[ig].molname, nullptr);
2317 sprintf(buf, "iontype%d-in-A", i);
2318 ITYPE(buf, ir->swap->grp[ig].nmolReq[0], -1);
2319 sprintf(buf, "iontype%d-in-B", i);
2320 ITYPE(buf, ir->swap->grp[ig].nmolReq[1], -1);
2323 CTYPE("By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2324 CTYPE("at maximum distance (= bulk concentration) to the split group layers. However,");
2325 CTYPE("an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2326 CTYPE("towards one of the compartment-partitioning layers (at +/- 1.0).");
2327 RTYPE("bulk-offsetA", ir->swap->bulkOffset[0], 0.0);
2328 RTYPE("bulk-offsetB", ir->swap->bulkOffset[1], 0.0);
2329 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2330 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2332 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2335 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2336 RTYPE("threshold", ir->swap->threshold, 1.0);
2339 /* AdResS is no longer supported, but we need mdrun to be able to refuse to run old AdResS .tpr files */
2340 EETYPE("adress", ir->bAdress, yesno_names);
2342 /* User defined thingies */
2343 CCTYPE ("User defined thingies");
2344 STYPE ("user1-grps", is->user1, nullptr);
2345 STYPE ("user2-grps", is->user2, nullptr);
2346 ITYPE ("userint1", ir->userint1, 0);
2347 ITYPE ("userint2", ir->userint2, 0);
2348 ITYPE ("userint3", ir->userint3, 0);
2349 ITYPE ("userint4", ir->userint4, 0);
2350 RTYPE ("userreal1", ir->userreal1, 0);
2351 RTYPE ("userreal2", ir->userreal2, 0);
2352 RTYPE ("userreal3", ir->userreal3, 0);
2353 RTYPE ("userreal4", ir->userreal4, 0);
2357 gmx::TextOutputFile stream(mdparout);
2358 write_inpfile(&stream, mdparout, ninp, inp, FALSE, writeMdpHeader, wi);
2362 for (i = 0; (i < ninp); i++)
2365 sfree(inp[i].value);
2369 /* Process options if necessary */
2370 for (m = 0; m < 2; m++)
2372 for (i = 0; i < 2*DIM; i++)
2381 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2383 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
2385 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2387 case epctSEMIISOTROPIC:
2388 case epctSURFACETENSION:
2389 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2391 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
2393 dumdub[m][YY] = dumdub[m][XX];
2395 case epctANISOTROPIC:
2396 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2397 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2398 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2400 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
2404 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2405 epcoupltype_names[ir->epct]);
2409 clear_mat(ir->ref_p);
2410 clear_mat(ir->compress);
2411 for (i = 0; i < DIM; i++)
2413 ir->ref_p[i][i] = dumdub[1][i];
2414 ir->compress[i][i] = dumdub[0][i];
2416 if (ir->epct == epctANISOTROPIC)
2418 ir->ref_p[XX][YY] = dumdub[1][3];
2419 ir->ref_p[XX][ZZ] = dumdub[1][4];
2420 ir->ref_p[YY][ZZ] = dumdub[1][5];
2421 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2423 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2425 ir->compress[XX][YY] = dumdub[0][3];
2426 ir->compress[XX][ZZ] = dumdub[0][4];
2427 ir->compress[YY][ZZ] = dumdub[0][5];
2428 for (i = 0; i < DIM; i++)
2430 for (m = 0; m < i; m++)
2432 ir->ref_p[i][m] = ir->ref_p[m][i];
2433 ir->compress[i][m] = ir->compress[m][i];
2438 if (ir->comm_mode == ecmNO)
2443 opts->couple_moltype = nullptr;
2444 if (strlen(is->couple_moltype) > 0)
2446 if (ir->efep != efepNO)
2448 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2449 if (opts->couple_lam0 == opts->couple_lam1)
2451 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2453 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2454 opts->couple_lam1 == ecouplamNONE))
2456 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2461 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2464 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2465 if (ir->efep != efepNO)
2467 if (fep->delta_lambda > 0)
2469 ir->efep = efepSLOWGROWTH;
2473 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2475 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2476 warning_note(wi, "Old option for dhdl-print-energy given: "
2477 "changing \"yes\" to \"total\"\n");
2480 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2482 /* always print out the energy to dhdl if we are doing
2483 expanded ensemble, since we need the total energy for
2484 analysis if the temperature is changing. In some
2485 conditions one may only want the potential energy, so
2486 we will allow that if the appropriate mdp setting has
2487 been enabled. Otherwise, total it is:
2489 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2492 if ((ir->efep != efepNO) || ir->bSimTemp)
2494 ir->bExpanded = FALSE;
2495 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2497 ir->bExpanded = TRUE;
2499 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2500 if (ir->bSimTemp) /* done after fep params */
2502 do_simtemp_params(ir);
2505 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2506 * setup and not on the old way of specifying the free-energy setup,
2507 * we should check for using soft-core when not needed, since that
2508 * can complicate the sampling significantly.
2509 * Note that we only check for the automated coupling setup.
2510 * If the (advanced) user does FEP through manual topology changes,
2511 * this check will not be triggered.
2513 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2514 ir->fepvals->sc_alpha != 0 &&
2515 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2516 couple_lambda_has_vdw_on(opts->couple_lam1)))
2518 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.");
2523 ir->fepvals->n_lambda = 0;
2526 /* WALL PARAMETERS */
2528 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2530 /* ORIENTATION RESTRAINT PARAMETERS */
2532 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, nullptr) != 1)
2534 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2537 /* DEFORMATION PARAMETERS */
2539 clear_mat(ir->deform);
2540 for (i = 0; i < 6; i++)
2545 double gmx_unused canary;
2546 int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
2547 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2548 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
2550 if (strlen(is->deform) > 0 && ndeform != 6)
2552 sprintf(warn_buf, "Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform);
2553 warning_error(wi, warn_buf);
2555 for (i = 0; i < 3; i++)
2557 ir->deform[i][i] = dumdub[0][i];
2559 ir->deform[YY][XX] = dumdub[0][3];
2560 ir->deform[ZZ][XX] = dumdub[0][4];
2561 ir->deform[ZZ][YY] = dumdub[0][5];
2562 if (ir->epc != epcNO)
2564 for (i = 0; i < 3; i++)
2566 for (j = 0; j <= i; j++)
2568 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2570 warning_error(wi, "A box element has deform set and compressibility > 0");
2574 for (i = 0; i < 3; i++)
2576 for (j = 0; j < i; j++)
2578 if (ir->deform[i][j] != 0)
2580 for (m = j; m < DIM; m++)
2582 if (ir->compress[m][j] != 0)
2584 sprintf(warn_buf, "An off-diagonal box element has deform set while compressibility > 0 for the same component of another box vector, this might lead to spurious periodicity effects.");
2585 warning(wi, warn_buf);
2593 /* Ion/water position swapping checks */
2594 if (ir->eSwapCoords != eswapNO)
2596 if (ir->swap->nstswap < 1)
2598 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2600 if (ir->swap->nAverage < 1)
2602 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2604 if (ir->swap->threshold < 1.0)
2606 warning_error(wi, "Ion count threshold must be at least 1.\n");
2614 static int search_QMstring(const char *s, int ng, const char *gn[])
2616 /* same as normal search_string, but this one searches QM strings */
2619 for (i = 0; (i < ng); i++)
2621 if (gmx_strcasecmp(s, gn[i]) == 0)
2627 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2631 } /* search_QMstring */
2633 /* We would like gn to be const as well, but C doesn't allow this */
2634 /* TODO this is utility functionality (search for the index of a
2635 string in a collection), so should be refactored and located more
2637 int search_string(const char *s, int ng, char *gn[])
2641 for (i = 0; (i < ng); i++)
2643 if (gmx_strcasecmp(s, gn[i]) == 0)
2650 "Group %s referenced in the .mdp file was not found in the index file.\n"
2651 "Group names must match either [moleculetype] names or custom index group\n"
2652 "names, in which case you must supply an index file to the '-n' option\n"
2659 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2660 t_blocka *block, char *gnames[],
2661 int gtype, int restnm,
2662 int grptp, gmx_bool bVerbose,
2665 unsigned short *cbuf;
2666 t_grps *grps = &(groups->grps[gtype]);
2667 int i, j, gid, aj, ognr, ntot = 0;
2670 char warn_buf[STRLEN];
2674 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2677 title = gtypes[gtype];
2680 /* Mark all id's as not set */
2681 for (i = 0; (i < natoms); i++)
2686 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2687 for (i = 0; (i < ng); i++)
2689 /* Lookup the group name in the block structure */
2690 gid = search_string(ptrs[i], block->nr, gnames);
2691 if ((grptp != egrptpONE) || (i == 0))
2693 grps->nm_ind[grps->nr++] = gid;
2697 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2700 /* Now go over the atoms in the group */
2701 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2706 /* Range checking */
2707 if ((aj < 0) || (aj >= natoms))
2709 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2711 /* Lookup up the old group number */
2715 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2716 aj+1, title, ognr+1, i+1);
2720 /* Store the group number in buffer */
2721 if (grptp == egrptpONE)
2734 /* Now check whether we have done all atoms */
2738 if (grptp == egrptpALL)
2740 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2741 natoms-ntot, title);
2743 else if (grptp == egrptpPART)
2745 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2746 natoms-ntot, title);
2747 warning_note(wi, warn_buf);
2749 /* Assign all atoms currently unassigned to a rest group */
2750 for (j = 0; (j < natoms); j++)
2752 if (cbuf[j] == NOGID)
2758 if (grptp != egrptpPART)
2763 "Making dummy/rest group for %s containing %d elements\n",
2764 title, natoms-ntot);
2766 /* Add group name "rest" */
2767 grps->nm_ind[grps->nr] = restnm;
2769 /* Assign the rest name to all atoms not currently assigned to a group */
2770 for (j = 0; (j < natoms); j++)
2772 if (cbuf[j] == NOGID)
2781 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2783 /* All atoms are part of one (or no) group, no index required */
2784 groups->ngrpnr[gtype] = 0;
2785 groups->grpnr[gtype] = nullptr;
2789 groups->ngrpnr[gtype] = natoms;
2790 snew(groups->grpnr[gtype], natoms);
2791 for (j = 0; (j < natoms); j++)
2793 groups->grpnr[gtype][j] = cbuf[j];
2799 return (bRest && grptp == egrptpPART);
2802 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2805 gmx_groups_t *groups;
2806 pull_params_t *pull;
2807 int natoms, ai, aj, i, j, d, g, imin, jmin;
2809 int *nrdf2, *na_vcm, na_tot;
2810 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2812 gmx_mtop_atomloop_all_t aloop;
2813 int mb, mol, ftype, as;
2814 gmx_molblock_t *molb;
2815 gmx_moltype_t *molt;
2818 * First calc 3xnr-atoms for each group
2819 * then subtract half a degree of freedom for each constraint
2821 * Only atoms and nuclei contribute to the degrees of freedom...
2826 groups = &mtop->groups;
2827 natoms = mtop->natoms;
2829 /* Allocate one more for a possible rest group */
2830 /* We need to sum degrees of freedom into doubles,
2831 * since floats give too low nrdf's above 3 million atoms.
2833 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2834 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2835 snew(dof_vcm, groups->grps[egcVCM].nr+1);
2836 snew(na_vcm, groups->grps[egcVCM].nr+1);
2837 snew(nrdf_vcm_sub, groups->grps[egcVCM].nr+1);
2839 for (i = 0; i < groups->grps[egcTC].nr; i++)
2843 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2846 clear_ivec(dof_vcm[i]);
2848 nrdf_vcm_sub[i] = 0;
2851 snew(nrdf2, natoms);
2852 aloop = gmx_mtop_atomloop_all_init(mtop);
2854 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2857 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2859 g = ggrpnr(groups, egcFREEZE, i);
2860 for (d = 0; d < DIM; d++)
2862 if (opts->nFreeze[g][d] == 0)
2864 /* Add one DOF for particle i (counted as 2*1) */
2866 /* VCM group i has dim d as a DOF */
2867 dof_vcm[ggrpnr(groups, egcVCM, i)][d] = 1;
2870 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2871 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2876 for (mb = 0; mb < mtop->nmolblock; mb++)
2878 molb = &mtop->molblock[mb];
2879 molt = &mtop->moltype[molb->type];
2880 atom = molt->atoms.atom;
2881 for (mol = 0; mol < molb->nmol; mol++)
2883 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2885 ia = molt->ilist[ftype].iatoms;
2886 for (i = 0; i < molt->ilist[ftype].nr; )
2888 /* Subtract degrees of freedom for the constraints,
2889 * if the particles still have degrees of freedom left.
2890 * If one of the particles is a vsite or a shell, then all
2891 * constraint motion will go there, but since they do not
2892 * contribute to the constraints the degrees of freedom do not
2897 if (((atom[ia[1]].ptype == eptNucleus) ||
2898 (atom[ia[1]].ptype == eptAtom)) &&
2899 ((atom[ia[2]].ptype == eptNucleus) ||
2900 (atom[ia[2]].ptype == eptAtom)))
2918 imin = std::min(imin, nrdf2[ai]);
2919 jmin = std::min(jmin, nrdf2[aj]);
2922 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2923 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2924 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2925 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2927 ia += interaction_function[ftype].nratoms+1;
2928 i += interaction_function[ftype].nratoms+1;
2931 ia = molt->ilist[F_SETTLE].iatoms;
2932 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2934 /* Subtract 1 dof from every atom in the SETTLE */
2935 for (j = 0; j < 3; j++)
2938 imin = std::min(2, nrdf2[ai]);
2940 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2941 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2946 as += molt->atoms.nr;
2952 /* Correct nrdf for the COM constraints.
2953 * We correct using the TC and VCM group of the first atom
2954 * in the reference and pull group. If atoms in one pull group
2955 * belong to different TC or VCM groups it is anyhow difficult
2956 * to determine the optimal nrdf assignment.
2960 for (i = 0; i < pull->ncoord; i++)
2962 if (pull->coord[i].eType != epullCONSTRAINT)
2969 for (j = 0; j < 2; j++)
2971 const t_pull_group *pgrp;
2973 pgrp = &pull->group[pull->coord[i].group[j]];
2977 /* Subtract 1/2 dof from each group */
2979 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2980 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2981 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2983 gmx_fatal(FARGS, "Center of mass pulling constraints caused the number of degrees of freedom for temperature coupling group %s to be negative", gnames[groups->grps[egcTC].nm_ind[ggrpnr(groups, egcTC, ai)]]);
2988 /* We need to subtract the whole DOF from group j=1 */
2995 if (ir->nstcomm != 0)
2999 /* We remove COM motion up to dim ndof_com() */
3000 ndim_rm_vcm = ndof_com(ir);
3002 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
3003 * the number of degrees of freedom in each vcm group when COM
3004 * translation is removed and 6 when rotation is removed as well.
3006 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3008 switch (ir->comm_mode)
3011 nrdf_vcm_sub[j] = 0;
3012 for (d = 0; d < ndim_rm_vcm; d++)
3021 nrdf_vcm_sub[j] = 6;
3024 gmx_incons("Checking comm_mode");
3028 for (i = 0; i < groups->grps[egcTC].nr; i++)
3030 /* Count the number of atoms of TC group i for every VCM group */
3031 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3036 for (ai = 0; ai < natoms; ai++)
3038 if (ggrpnr(groups, egcTC, ai) == i)
3040 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
3044 /* Correct for VCM removal according to the fraction of each VCM
3045 * group present in this TC group.
3047 nrdf_uc = nrdf_tc[i];
3050 fprintf(debug, "T-group[%d] nrdf_uc = %g\n", i, nrdf_uc);
3053 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
3055 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
3057 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
3058 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
3062 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
3063 j, nrdf_vcm[j], nrdf_tc[i]);
3068 for (i = 0; (i < groups->grps[egcTC].nr); i++)
3070 opts->nrdf[i] = nrdf_tc[i];
3071 if (opts->nrdf[i] < 0)
3076 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
3077 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
3085 sfree(nrdf_vcm_sub);
3088 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
3089 const char *option, const char *val, int flag)
3091 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3092 * But since this is much larger than STRLEN, such a line can not be parsed.
3093 * The real maximum is the number of names that fit in a string: STRLEN/2.
3095 #define EGP_MAX (STRLEN/2)
3096 int nelem, i, j, k, nr;
3097 char *names[EGP_MAX];
3101 gnames = groups->grpname;
3103 nelem = str_nelem(val, EGP_MAX, names);
3106 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3108 nr = groups->grps[egcENER].nr;
3110 for (i = 0; i < nelem/2; i++)
3114 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
3120 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3121 names[2*i], option);
3125 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3131 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3132 names[2*i+1], option);
3134 if ((j < nr) && (k < nr))
3136 ir->opts.egp_flags[nr*j+k] |= flag;
3137 ir->opts.egp_flags[nr*k+j] |= flag;
3146 static void make_swap_groups(
3151 int ig = -1, i = 0, gind;
3155 /* Just a quick check here, more thorough checks are in mdrun */
3156 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3158 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3161 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3162 for (ig = 0; ig < swap->ngrp; ig++)
3164 swapg = &swap->grp[ig];
3165 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3166 swapg->nat = grps->index[gind+1] - grps->index[gind];
3170 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3171 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3172 swap->grp[ig].molname, swapg->nat);
3173 snew(swapg->ind, swapg->nat);
3174 for (i = 0; i < swapg->nat; i++)
3176 swapg->ind[i] = grps->a[grps->index[gind]+i];
3181 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3187 void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3192 ig = search_string(IMDgname, grps->nr, gnames);
3193 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3195 if (IMDgroup->nat > 0)
3197 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3198 IMDgname, IMDgroup->nat);
3199 snew(IMDgroup->ind, IMDgroup->nat);
3200 for (i = 0; i < IMDgroup->nat; i++)
3202 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3208 void do_index(const char* mdparin, const char *ndx,
3215 gmx_groups_t *groups;
3219 char warnbuf[STRLEN], **gnames;
3220 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3223 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3224 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3225 int i, j, k, restnm;
3226 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3227 int nQMmethod, nQMbasis, nQMg;
3228 char warn_buf[STRLEN];
3233 fprintf(stderr, "processing index file...\n");
3238 snew(grps->index, 1);
3240 atoms_all = gmx_mtop_global_atoms(mtop);
3241 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3242 done_atom(&atoms_all);
3246 grps = init_index(ndx, &gnames);
3249 groups = &mtop->groups;
3250 natoms = mtop->natoms;
3251 symtab = &mtop->symtab;
3253 snew(groups->grpname, grps->nr+1);
3255 for (i = 0; (i < grps->nr); i++)
3257 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3259 groups->grpname[i] = put_symtab(symtab, "rest");
3261 srenew(gnames, grps->nr+1);
3262 gnames[restnm] = *(groups->grpname[i]);
3263 groups->ngrpname = grps->nr+1;
3265 set_warning_line(wi, mdparin, -1);
3267 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3268 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3269 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3270 if ((ntau_t != ntcg) || (nref_t != ntcg))
3272 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3273 "%d tau-t values", ntcg, nref_t, ntau_t);
3276 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3277 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3278 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3279 nr = groups->grps[egcTC].nr;
3281 snew(ir->opts.nrdf, nr);
3282 snew(ir->opts.tau_t, nr);
3283 snew(ir->opts.ref_t, nr);
3284 if (ir->eI == eiBD && ir->bd_fric == 0)
3286 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3293 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3297 for (i = 0; (i < nr); i++)
3299 ir->opts.tau_t[i] = strtod(ptr1[i], &endptr);
3302 warning_error(wi, "Invalid value for mdp option tau-t. tau-t should only consist of real numbers separated by spaces.");
3304 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3306 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3307 warning_error(wi, warn_buf);
3310 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3312 warning_note(wi, "tau-t = -1 is the value to signal that a group should not have temperature coupling. Treating your use of tau-t = 0 as if you used -1.");
3315 if (ir->opts.tau_t[i] >= 0)
3317 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3320 if (ir->etc != etcNO && ir->nsttcouple == -1)
3322 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3327 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3329 gmx_fatal(FARGS, "Cannot do Nose-Hoover temperature with Berendsen pressure control with md-vv; use either vrescale temperature with berendsen pressure or Nose-Hoover temperature with MTTK pressure");
3331 if (ir->epc == epcMTTK)
3333 if (ir->etc != etcNOSEHOOVER)
3335 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3339 if (ir->nstpcouple != ir->nsttcouple)
3341 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3342 ir->nstpcouple = ir->nsttcouple = mincouple;
3343 sprintf(warn_buf, "for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple) = %d", mincouple);
3344 warning_note(wi, warn_buf);
3349 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3350 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3352 if (ETC_ANDERSEN(ir->etc))
3354 if (ir->nsttcouple != 1)
3357 sprintf(warn_buf, "Andersen temperature control methods assume nsttcouple = 1; there is no need for larger nsttcouple > 1, since no global parameters are computed. nsttcouple has been reset to 1");
3358 warning_note(wi, warn_buf);
3361 nstcmin = tcouple_min_integration_steps(ir->etc);
3364 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3366 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3367 ETCOUPLTYPE(ir->etc),
3369 ir->nsttcouple*ir->delta_t);
3370 warning(wi, warn_buf);
3373 for (i = 0; (i < nr); i++)
3375 ir->opts.ref_t[i] = strtod(ptr2[i], &endptr);
3378 warning_error(wi, "Invalid value for mdp option ref-t. ref-t should only consist of real numbers separated by spaces.");
3380 if (ir->opts.ref_t[i] < 0)
3382 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3385 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3386 if we are in this conditional) if mc_temp is negative */
3387 if (ir->expandedvals->mc_temp < 0)
3389 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3393 /* Simulated annealing for each group. There are nr groups */
3394 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3395 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3399 if (nSA > 0 && nSA != nr)
3401 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3405 snew(ir->opts.annealing, nr);
3406 snew(ir->opts.anneal_npoints, nr);
3407 snew(ir->opts.anneal_time, nr);
3408 snew(ir->opts.anneal_temp, nr);
3409 for (i = 0; i < nr; i++)
3411 ir->opts.annealing[i] = eannNO;
3412 ir->opts.anneal_npoints[i] = 0;
3413 ir->opts.anneal_time[i] = nullptr;
3414 ir->opts.anneal_temp[i] = nullptr;
3419 for (i = 0; i < nr; i++)
3421 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3423 ir->opts.annealing[i] = eannNO;
3425 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3427 ir->opts.annealing[i] = eannSINGLE;
3430 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3432 ir->opts.annealing[i] = eannPERIODIC;
3438 /* Read the other fields too */
3439 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3440 if (nSA_points != nSA)
3442 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3444 for (k = 0, i = 0; i < nr; i++)
3446 ir->opts.anneal_npoints[i] = strtol(ptr1[i], &endptr, 10);
3449 warning_error(wi, "Invalid value for mdp option annealing-npoints. annealing should only consist of integers separated by spaces.");
3451 if (ir->opts.anneal_npoints[i] == 1)
3453 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3455 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3456 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3457 k += ir->opts.anneal_npoints[i];
3460 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3463 gmx_fatal(FARGS, "Found %d annealing-time values, wanted %d\n", nSA_time, k);
3465 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3468 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3471 for (i = 0, k = 0; i < nr; i++)
3474 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3476 ir->opts.anneal_time[i][j] = strtod(ptr1[k], &endptr);
3479 warning_error(wi, "Invalid value for mdp option anneal-time. anneal-time should only consist of real numbers separated by spaces.");
3481 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], &endptr);
3484 warning_error(wi, "Invalid value for anneal-temp. anneal-temp should only consist of real numbers separated by spaces.");
3488 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3490 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3496 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3498 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3499 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3502 if (ir->opts.anneal_temp[i][j] < 0)
3504 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3509 /* Print out some summary information, to make sure we got it right */
3510 for (i = 0, k = 0; i < nr; i++)
3512 if (ir->opts.annealing[i] != eannNO)
3514 j = groups->grps[egcTC].nm_ind[i];
3515 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3516 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3517 ir->opts.anneal_npoints[i]);
3518 fprintf(stderr, "Time (ps) Temperature (K)\n");
3519 /* All terms except the last one */
3520 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3522 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3525 /* Finally the last one */
3526 j = ir->opts.anneal_npoints[i]-1;
3527 if (ir->opts.annealing[i] == eannSINGLE)
3529 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3533 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3534 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3536 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3547 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3549 make_pull_coords(ir->pull);
3554 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3557 if (ir->eSwapCoords != eswapNO)
3559 make_swap_groups(ir->swap, grps, gnames);
3562 /* Make indices for IMD session */
3565 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3568 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3569 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3570 if (nacg*DIM != nacc)
3572 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3575 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3576 restnm, egrptpALL_GENREST, bVerbose, wi);
3577 nr = groups->grps[egcACC].nr;
3578 snew(ir->opts.acc, nr);
3579 ir->opts.ngacc = nr;
3581 for (i = k = 0; (i < nacg); i++)
3583 for (j = 0; (j < DIM); j++, k++)
3585 ir->opts.acc[i][j] = strtod(ptr1[k], &endptr);
3588 warning_error(wi, "Invalid value for mdp option accelerate. accelerate should only consist of real numbers separated by spaces.");
3592 for (; (i < nr); i++)
3594 for (j = 0; (j < DIM); j++)
3596 ir->opts.acc[i][j] = 0;
3600 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3601 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3602 if (nfrdim != DIM*nfreeze)
3604 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3607 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3608 restnm, egrptpALL_GENREST, bVerbose, wi);
3609 nr = groups->grps[egcFREEZE].nr;
3610 ir->opts.ngfrz = nr;
3611 snew(ir->opts.nFreeze, nr);
3612 for (i = k = 0; (i < nfreeze); i++)
3614 for (j = 0; (j < DIM); j++, k++)
3616 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3617 if (!ir->opts.nFreeze[i][j])
3619 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3621 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3622 "(not %s)", ptr1[k]);
3623 warning(wi, warn_buf);
3628 for (; (i < nr); i++)
3630 for (j = 0; (j < DIM); j++)
3632 ir->opts.nFreeze[i][j] = 0;
3636 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3637 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3638 restnm, egrptpALL_GENREST, bVerbose, wi);
3639 add_wall_energrps(groups, ir->nwall, symtab);
3640 ir->opts.ngener = groups->grps[egcENER].nr;
3641 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3643 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3644 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3647 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3648 "This may lead to artifacts.\n"
3649 "In most cases one should use one group for the whole system.");
3652 /* Now we have filled the freeze struct, so we can calculate NRDF */
3653 calc_nrdf(mtop, ir, gnames);
3655 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3656 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3657 restnm, egrptpALL_GENREST, bVerbose, wi);
3658 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3659 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3660 restnm, egrptpALL_GENREST, bVerbose, wi);
3661 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3662 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3663 restnm, egrptpONE, bVerbose, wi);
3664 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3665 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3666 restnm, egrptpALL_GENREST, bVerbose, wi);
3668 /* QMMM input processing */
3669 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3670 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3671 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3672 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3674 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3675 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3677 /* group rest, if any, is always MM! */
3678 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3679 restnm, egrptpALL_GENREST, bVerbose, wi);
3680 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3681 ir->opts.ngQM = nQMg;
3682 snew(ir->opts.QMmethod, nr);
3683 snew(ir->opts.QMbasis, nr);
3684 for (i = 0; i < nr; i++)
3686 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3687 * converted to the corresponding enum in names.c
3689 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3691 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3695 str_nelem(is->QMmult, MAXPTR, ptr1);
3696 str_nelem(is->QMcharge, MAXPTR, ptr2);
3697 str_nelem(is->bSH, MAXPTR, ptr3);
3698 snew(ir->opts.QMmult, nr);
3699 snew(ir->opts.QMcharge, nr);
3700 snew(ir->opts.bSH, nr);
3702 for (i = 0; i < nr; i++)
3704 ir->opts.QMmult[i] = strtol(ptr1[i], &endptr, 10);
3707 warning_error(wi, "Invalid value for mdp option QMmult. QMmult should only consist of integers separated by spaces.");
3709 ir->opts.QMcharge[i] = strtol(ptr2[i], &endptr, 10);
3712 warning_error(wi, "Invalid value for mdp option QMcharge. QMcharge should only consist of integers separated by spaces.");
3714 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3717 str_nelem(is->CASelectrons, MAXPTR, ptr1);
3718 str_nelem(is->CASorbitals, MAXPTR, ptr2);
3719 snew(ir->opts.CASelectrons, nr);
3720 snew(ir->opts.CASorbitals, nr);
3721 for (i = 0; i < nr; i++)
3723 ir->opts.CASelectrons[i] = strtol(ptr1[i], &endptr, 10);
3726 warning_error(wi, "Invalid value for mdp option CASelectrons. CASelectrons should only consist of integers separated by spaces.");
3728 ir->opts.CASorbitals[i] = strtol(ptr2[i], &endptr, 10);
3731 warning_error(wi, "Invalid value for mdp option CASorbitals. CASorbitals should only consist of integers separated by spaces.");
3734 /* special optimization options */
3736 str_nelem(is->bOPT, MAXPTR, ptr1);
3737 str_nelem(is->bTS, MAXPTR, ptr2);
3738 snew(ir->opts.bOPT, nr);
3739 snew(ir->opts.bTS, nr);
3740 for (i = 0; i < nr; i++)
3742 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3743 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3745 str_nelem(is->SAon, MAXPTR, ptr1);
3746 str_nelem(is->SAoff, MAXPTR, ptr2);
3747 str_nelem(is->SAsteps, MAXPTR, ptr3);
3748 snew(ir->opts.SAon, nr);
3749 snew(ir->opts.SAoff, nr);
3750 snew(ir->opts.SAsteps, nr);
3752 for (i = 0; i < nr; i++)
3754 ir->opts.SAon[i] = strtod(ptr1[i], &endptr);
3757 warning_error(wi, "Invalid value for mdp option SAon. SAon should only consist of real numbers separated by spaces.");
3759 ir->opts.SAoff[i] = strtod(ptr2[i], &endptr);
3762 warning_error(wi, "Invalid value for mdp option SAoff. SAoff should only consist of real numbers separated by spaces.");
3764 ir->opts.SAsteps[i] = strtol(ptr3[i], &endptr, 10);
3767 warning_error(wi, "Invalid value for mdp option SAsteps. SAsteps should only consist of integers separated by spaces.");
3770 /* end of QMMM input */
3774 for (i = 0; (i < egcNR); i++)
3776 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3777 for (j = 0; (j < groups->grps[i].nr); j++)
3779 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3781 fprintf(stderr, "\n");
3785 nr = groups->grps[egcENER].nr;
3786 snew(ir->opts.egp_flags, nr*nr);
3788 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3789 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3791 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3793 if (bExcl && EEL_FULL(ir->coulombtype))
3795 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3798 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3799 if (bTable && !(ir->vdwtype == evdwUSER) &&
3800 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3801 !(ir->coulombtype == eelPMEUSERSWITCH))
3803 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3806 for (i = 0; (i < grps->nr); i++)
3818 static void check_disre(gmx_mtop_t *mtop)
3820 gmx_ffparams_t *ffparams;
3821 t_functype *functype;
3823 int i, ndouble, ftype;
3824 int label, old_label;
3826 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3828 ffparams = &mtop->ffparams;
3829 functype = ffparams->functype;
3830 ip = ffparams->iparams;
3833 for (i = 0; i < ffparams->ntypes; i++)
3835 ftype = functype[i];
3836 if (ftype == F_DISRES)
3838 label = ip[i].disres.label;
3839 if (label == old_label)
3841 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3849 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3850 "probably the parameters for multiple pairs in one restraint "
3851 "are not identical\n", ndouble);
3856 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3857 gmx_bool posres_only,
3861 gmx_mtop_ilistloop_t iloop;
3871 for (d = 0; d < DIM; d++)
3873 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3875 /* Check for freeze groups */
3876 for (g = 0; g < ir->opts.ngfrz; g++)
3878 for (d = 0; d < DIM; d++)
3880 if (ir->opts.nFreeze[g][d] != 0)
3888 /* Check for position restraints */
3889 iloop = gmx_mtop_ilistloop_init(sys);
3890 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3893 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3895 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3897 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3898 for (d = 0; d < DIM; d++)
3900 if (pr->posres.fcA[d] != 0)
3906 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3908 /* Check for flat-bottom posres */
3909 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3910 if (pr->fbposres.k != 0)
3912 switch (pr->fbposres.geom)
3914 case efbposresSPHERE:
3915 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3917 case efbposresCYLINDERX:
3918 AbsRef[YY] = AbsRef[ZZ] = 1;
3920 case efbposresCYLINDERY:
3921 AbsRef[XX] = AbsRef[ZZ] = 1;
3923 case efbposresCYLINDER:
3924 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3925 case efbposresCYLINDERZ:
3926 AbsRef[XX] = AbsRef[YY] = 1;
3928 case efbposresX: /* d=XX */
3929 case efbposresY: /* d=YY */
3930 case efbposresZ: /* d=ZZ */
3931 d = pr->fbposres.geom - efbposresX;
3935 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3936 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3944 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3948 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3949 gmx_bool *bC6ParametersWorkWithGeometricRules,
3950 gmx_bool *bC6ParametersWorkWithLBRules,
3951 gmx_bool *bLBRulesPossible)
3953 int ntypes, tpi, tpj;
3956 double c6i, c6j, c12i, c12j;
3957 double c6, c6_geometric, c6_LB;
3958 double sigmai, sigmaj, epsi, epsj;
3959 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3962 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3963 * force-field floating point parameters.
3966 ptr = getenv("GMX_LJCOMB_TOL");
3970 double gmx_unused canary;
3972 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3974 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
3979 *bC6ParametersWorkWithLBRules = TRUE;
3980 *bC6ParametersWorkWithGeometricRules = TRUE;
3981 bCanDoLBRules = TRUE;
3982 ntypes = mtop->ffparams.atnr;
3983 snew(typecount, ntypes);
3984 gmx_mtop_count_atomtypes(mtop, state, typecount);
3985 *bLBRulesPossible = TRUE;
3986 for (tpi = 0; tpi < ntypes; ++tpi)
3988 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3989 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3990 for (tpj = tpi; tpj < ntypes; ++tpj)
3992 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3993 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3994 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3995 c6_geometric = std::sqrt(c6i * c6j);
3996 if (!gmx_numzero(c6_geometric))
3998 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
4000 sigmai = gmx::sixthroot(c12i / c6i);
4001 sigmaj = gmx::sixthroot(c12j / c6j);
4002 epsi = c6i * c6i /(4.0 * c12i);
4003 epsj = c6j * c6j /(4.0 * c12j);
4004 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
4008 *bLBRulesPossible = FALSE;
4009 c6_LB = c6_geometric;
4011 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
4014 if (FALSE == bCanDoLBRules)
4016 *bC6ParametersWorkWithLBRules = FALSE;
4019 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
4021 if (FALSE == bCanDoGeometricRules)
4023 *bC6ParametersWorkWithGeometricRules = FALSE;
4031 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
4034 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
4036 check_combination_rule_differences(mtop, 0,
4037 &bC6ParametersWorkWithGeometricRules,
4038 &bC6ParametersWorkWithLBRules,
4040 if (ir->ljpme_combination_rule == eljpmeLB)
4042 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
4044 warning(wi, "You are using arithmetic-geometric combination rules "
4045 "in LJ-PME, but your non-bonded C6 parameters do not "
4046 "follow these rules.");
4051 if (FALSE == bC6ParametersWorkWithGeometricRules)
4053 if (ir->eDispCorr != edispcNO)
4055 warning_note(wi, "You are using geometric combination rules in "
4056 "LJ-PME, but your non-bonded C6 parameters do "
4057 "not follow these rules. "
4058 "This will introduce very small errors in the forces and energies in "
4059 "your simulations. Dispersion correction will correct total energy "
4060 "and/or pressure for isotropic systems, but not forces or surface tensions.");
4064 warning_note(wi, "You are using geometric combination rules in "
4065 "LJ-PME, but your non-bonded C6 parameters do "
4066 "not follow these rules. "
4067 "This will introduce very small errors in the forces and energies in "
4068 "your simulations. If your system is homogeneous, consider using dispersion correction "
4069 "for the total energy and pressure.");
4075 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
4078 char err_buf[STRLEN];
4080 gmx_bool bCharge, bAcc;
4083 gmx_mtop_atomloop_block_t aloopb;
4084 gmx_mtop_atomloop_all_t aloop;
4086 char warn_buf[STRLEN];
4088 set_warning_line(wi, mdparin, -1);
4090 if (ir->cutoff_scheme == ecutsVERLET &&
4091 ir->verletbuf_tol > 0 &&
4093 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
4094 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4096 /* Check if a too small Verlet buffer might potentially
4097 * cause more drift than the thermostat can couple off.
4099 /* Temperature error fraction for warning and suggestion */
4100 const real T_error_warn = 0.002;
4101 const real T_error_suggest = 0.001;
4102 /* For safety: 2 DOF per atom (typical with constraints) */
4103 const real nrdf_at = 2;
4104 real T, tau, max_T_error;
4109 for (i = 0; i < ir->opts.ngtc; i++)
4111 T = std::max(T, ir->opts.ref_t[i]);
4112 tau = std::max(tau, ir->opts.tau_t[i]);
4116 /* This is a worst case estimate of the temperature error,
4117 * assuming perfect buffer estimation and no cancelation
4118 * of errors. The factor 0.5 is because energy distributes
4119 * equally over Ekin and Epot.
4121 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4122 if (max_T_error > T_error_warn)
4124 sprintf(warn_buf, "With a verlet-buffer-tolerance of %g kJ/mol/ps, a reference temperature of %g and a tau_t of %g, your temperature might be off by up to %.1f%%. To ensure the error is below %.1f%%, decrease verlet-buffer-tolerance to %.0e or decrease tau_t.",
4125 ir->verletbuf_tol, T, tau,
4127 100*T_error_suggest,
4128 ir->verletbuf_tol*T_error_suggest/max_T_error);
4129 warning(wi, warn_buf);
4134 if (ETC_ANDERSEN(ir->etc))
4138 for (i = 0; i < ir->opts.ngtc; i++)
4140 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4141 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4142 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
4143 i, ir->opts.tau_t[i]);
4144 CHECK(ir->opts.tau_t[i] < 0);
4147 for (i = 0; i < ir->opts.ngtc; i++)
4149 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
4150 sprintf(err_buf, "tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization", i, etcoupl_names[ir->etc], ir->nstcomm, ir->opts.tau_t[i], nsteps);
4151 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
4155 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4156 ir->comm_mode == ecmNO &&
4157 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4158 !ETC_ANDERSEN(ir->etc))
4160 warning(wi, "You are not using center of mass motion removal (mdp option comm-mode), numerical rounding errors can lead to build up of kinetic energy of the center of mass");
4163 /* Check for pressure coupling with absolute position restraints */
4164 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4166 absolute_reference(ir, sys, TRUE, AbsRef);
4168 for (m = 0; m < DIM; m++)
4170 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4172 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4180 aloopb = gmx_mtop_atomloop_block_init(sys);
4182 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4184 if (atom->q != 0 || atom->qB != 0)
4192 if (EEL_FULL(ir->coulombtype))
4195 "You are using full electrostatics treatment %s for a system without charges.\n"
4196 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4197 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4198 warning(wi, err_buf);
4203 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4206 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4207 "You might want to consider using %s electrostatics.\n",
4209 warning_note(wi, err_buf);
4213 /* Check if combination rules used in LJ-PME are the same as in the force field */
4214 if (EVDW_PME(ir->vdwtype))
4216 check_combination_rules(ir, sys, wi);
4219 /* Generalized reaction field */
4220 if (ir->opts.ngtc == 0)
4222 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4224 CHECK(ir->coulombtype == eelGRF);
4228 sprintf(err_buf, "When using coulombtype = %s"
4229 " ref-t for temperature coupling should be > 0",
4231 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4235 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4237 for (m = 0; (m < DIM); m++)
4239 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4248 snew(mgrp, sys->groups.grps[egcACC].nr);
4249 aloop = gmx_mtop_atomloop_all_init(sys);
4251 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4253 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4256 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4258 for (m = 0; (m < DIM); m++)
4260 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4264 for (m = 0; (m < DIM); m++)
4266 if (fabs(acc[m]) > 1e-6)
4268 const char *dim[DIM] = { "X", "Y", "Z" };
4270 "Net Acceleration in %s direction, will %s be corrected\n",
4271 dim[m], ir->nstcomm != 0 ? "" : "not");
4272 if (ir->nstcomm != 0 && m < ndof_com(ir))
4275 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4277 ir->opts.acc[i][m] -= acc[m];
4285 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4286 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4288 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4296 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4298 if (ir->pull->coord[i].group[0] == 0 ||
4299 ir->pull->coord[i].group[1] == 0)
4301 absolute_reference(ir, sys, FALSE, AbsRef);
4302 for (m = 0; m < DIM; m++)
4304 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4306 warning(wi, "You are using an absolute reference for pulling, but the rest of the system does not have an absolute reference. This will lead to artifacts.");
4314 for (i = 0; i < 3; i++)
4316 for (m = 0; m <= i; m++)
4318 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4319 ir->deform[i][m] != 0)
4321 for (c = 0; c < ir->pull->ncoord; c++)
4323 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4324 ir->pull->coord[c].vec[m] != 0)
4326 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4337 void double_check(t_inputrec *ir, matrix box,
4338 gmx_bool bHasNormalConstraints,
4339 gmx_bool bHasAnyConstraints,
4343 char warn_buf[STRLEN];
4346 ptr = check_box(ir->ePBC, box);
4349 warning_error(wi, ptr);
4352 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4354 if (ir->shake_tol <= 0.0)
4356 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4358 warning_error(wi, warn_buf);
4362 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4364 /* If we have Lincs constraints: */
4365 if (ir->eI == eiMD && ir->etc == etcNO &&
4366 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4368 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4369 warning_note(wi, warn_buf);
4372 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4374 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4375 warning_note(wi, warn_buf);
4377 if (ir->epc == epcMTTK)
4379 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4383 if (bHasAnyConstraints && ir->epc == epcMTTK)
4385 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4388 if (ir->LincsWarnAngle > 90.0)
4390 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4391 warning(wi, warn_buf);
4392 ir->LincsWarnAngle = 90.0;
4395 if (ir->ePBC != epbcNONE)
4397 if (ir->nstlist == 0)
4399 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4401 if (ir->ns_type == ensGRID)
4403 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4405 sprintf(warn_buf, "ERROR: The cut-off length is longer than half the shortest box vector or longer than the smallest box diagonal element. Increase the box size or decrease rlist.\n");
4406 warning_error(wi, warn_buf);
4411 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4412 if (2*ir->rlist >= min_size)
4414 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4415 warning_error(wi, warn_buf);
4418 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4425 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4429 real rvdw1, rvdw2, rcoul1, rcoul2;
4430 char warn_buf[STRLEN];
4432 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4436 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4441 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4447 if (rvdw1 + rvdw2 > ir->rlist ||
4448 rcoul1 + rcoul2 > ir->rlist)
4451 "The sum of the two largest charge group radii (%f) "
4452 "is larger than rlist (%f)\n",
4453 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4454 warning(wi, warn_buf);
4458 /* Here we do not use the zero at cut-off macro,
4459 * since user defined interactions might purposely
4460 * not be zero at the cut-off.
4462 if (ir_vdw_is_zero_at_cutoff(ir) &&
4463 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4465 sprintf(warn_buf, "The sum of the two largest charge group "
4466 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4467 "With exact cut-offs, better performance can be "
4468 "obtained with cutoff-scheme = %s, because it "
4469 "does not use charge groups at all.",
4471 ir->rlist, ir->rvdw,
4472 ecutscheme_names[ecutsVERLET]);
4475 warning(wi, warn_buf);
4479 warning_note(wi, warn_buf);
4482 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4483 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4485 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4486 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4488 ir->rlist, ir->rcoulomb,
4489 ecutscheme_names[ecutsVERLET]);
4492 warning(wi, warn_buf);
4496 warning_note(wi, warn_buf);