2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/awh/read-params.h"
50 #include "gromacs/fileio/readinp.h"
51 #include "gromacs/fileio/warninp.h"
52 #include "gromacs/gmxlib/chargegroup.h"
53 #include "gromacs/gmxlib/network.h"
54 #include "gromacs/gmxpreprocess/keyvaluetreemdpwriter.h"
55 #include "gromacs/gmxpreprocess/toputil.h"
56 #include "gromacs/math/functions.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/mdlib/calc_verletbuf.h"
60 #include "gromacs/mdrunutility/mdmodules.h"
61 #include "gromacs/mdtypes/inputrec.h"
62 #include "gromacs/mdtypes/md_enums.h"
63 #include "gromacs/mdtypes/pull-params.h"
64 #include "gromacs/options/options.h"
65 #include "gromacs/options/treesupport.h"
66 #include "gromacs/pbcutil/pbc.h"
67 #include "gromacs/topology/block.h"
68 #include "gromacs/topology/ifunc.h"
69 #include "gromacs/topology/index.h"
70 #include "gromacs/topology/mtop_util.h"
71 #include "gromacs/topology/symtab.h"
72 #include "gromacs/topology/topology.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/exceptions.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/filestream.h"
77 #include "gromacs/utility/gmxassert.h"
78 #include "gromacs/utility/ikeyvaluetreeerror.h"
79 #include "gromacs/utility/keyvaluetree.h"
80 #include "gromacs/utility/keyvaluetreebuilder.h"
81 #include "gromacs/utility/keyvaluetreetransform.h"
82 #include "gromacs/utility/smalloc.h"
83 #include "gromacs/utility/stringcompare.h"
84 #include "gromacs/utility/stringutil.h"
85 #include "gromacs/utility/textwriter.h"
90 /* Resource parameters
91 * Do not change any of these until you read the instruction
92 * in readinp.h. Some cpp's do not take spaces after the backslash
93 * (like the c-shell), which will give you a very weird compiler
97 typedef struct t_inputrec_strings
99 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
100 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
101 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
102 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
103 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
105 char fep_lambda[efptNR][STRLEN];
106 char lambda_weights[STRLEN];
109 char anneal[STRLEN], anneal_npoints[STRLEN],
110 anneal_time[STRLEN], anneal_temp[STRLEN];
111 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
112 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
113 SAoff[STRLEN], SAsteps[STRLEN];
115 } gmx_inputrec_strings;
117 static gmx_inputrec_strings *is = nullptr;
119 void init_inputrec_strings()
123 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.");
128 void done_inputrec_strings()
136 egrptpALL, /* All particles have to be a member of a group. */
137 egrptpALL_GENREST, /* A rest group with name is generated for particles *
138 * that are not part of any group. */
139 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
140 * for the rest group. */
141 egrptpONE /* Merge all selected groups into one group, *
142 * make a rest group for the remaining particles. */
145 static const char *constraints[eshNR+1] = {
146 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", nullptr
149 static const char *couple_lam[ecouplamNR+1] = {
150 "vdw-q", "vdw", "q", "none", nullptr
153 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
158 for (i = 0; i < ntemps; i++)
160 /* simple linear scaling -- allows more control */
161 if (simtemp->eSimTempScale == esimtempLINEAR)
163 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
165 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
167 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
169 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
171 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
176 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
177 gmx_fatal(FARGS, "%s", errorstr);
184 static void _low_check(bool b, const char *s, warninp_t wi)
188 warning_error(wi, s);
192 static void check_nst(const char *desc_nst, int nst,
193 const char *desc_p, int *p,
198 if (*p > 0 && *p % nst != 0)
200 /* Round up to the next multiple of nst */
201 *p = ((*p)/nst + 1)*nst;
202 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
203 desc_p, desc_nst, desc_p, *p);
208 static bool ir_NVE(const t_inputrec *ir)
210 return (EI_MD(ir->eI) && ir->etc == etcNO);
213 static int lcd(int n1, int n2)
218 for (i = 2; (i <= n1 && i <= n2); i++)
220 if (n1 % i == 0 && n2 % i == 0)
229 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
231 if (*eintmod == eintmodPOTSHIFT_VERLET)
233 if (ir->cutoff_scheme == ecutsVERLET)
235 *eintmod = eintmodPOTSHIFT;
239 *eintmod = eintmodNONE;
244 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
246 /* Check internal consistency.
247 * NOTE: index groups are not set here yet, don't check things
248 * like temperature coupling group options here, but in triple_check
251 /* Strange macro: first one fills the err_buf, and then one can check
252 * the condition, which will print the message and increase the error
255 #define CHECK(b) _low_check(b, err_buf, wi)
256 char err_buf[256], warn_buf[STRLEN];
259 t_lambda *fep = ir->fepvals;
260 t_expanded *expand = ir->expandedvals;
262 set_warning_line(wi, mdparin, -1);
264 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
266 sprintf(warn_buf, "%s electrostatics is no longer supported",
267 eel_names[eelRF_NEC_UNSUPPORTED]);
268 warning_error(wi, warn_buf);
271 /* BASIC CUT-OFF STUFF */
272 if (ir->rcoulomb < 0)
274 warning_error(wi, "rcoulomb should be >= 0");
278 warning_error(wi, "rvdw should be >= 0");
281 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
283 warning_error(wi, "rlist should be >= 0");
285 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.)");
286 CHECK(ir->nstlist < 0);
288 process_interaction_modifier(ir, &ir->coulomb_modifier);
289 process_interaction_modifier(ir, &ir->vdw_modifier);
291 if (ir->cutoff_scheme == ecutsGROUP)
294 "The group cutoff scheme is deprecated since GROMACS 5.0 and will be removed in a future "
295 "release when all interaction forms are supported for the verlet scheme. The verlet "
296 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
298 if (ir->rlist > 0 && ir->rlist < ir->rcoulomb)
300 gmx_fatal(FARGS, "rcoulomb must not be greater than rlist (twin-range schemes are not supported)");
302 if (ir->rlist > 0 && ir->rlist < ir->rvdw)
304 gmx_fatal(FARGS, "rvdw must not be greater than rlist (twin-range schemes are not supported)");
307 if (ir->rlist == 0 && ir->ePBC != epbcNONE)
309 warning_error(wi, "Can not have an infinite cut-off with PBC");
313 if (ir->cutoff_scheme == ecutsVERLET)
317 /* Normal Verlet type neighbor-list, currently only limited feature support */
318 if (inputrec2nboundeddim(ir) < 3)
320 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
323 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
324 if (ir->rcoulomb != ir->rvdw)
326 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
327 // for PME load balancing, we can support this exception.
328 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
329 ir->vdwtype == evdwCUT &&
330 ir->rcoulomb > ir->rvdw);
331 if (!bUsesPmeTwinRangeKernel)
333 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
337 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
339 if (ir->vdw_modifier == eintmodNONE ||
340 ir->vdw_modifier == eintmodPOTSHIFT)
342 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
344 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]);
345 warning_note(wi, warn_buf);
347 ir->vdwtype = evdwCUT;
351 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
352 warning_error(wi, warn_buf);
356 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
358 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
360 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
361 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
363 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
365 if (!(ir->coulomb_modifier == eintmodNONE ||
366 ir->coulomb_modifier == eintmodPOTSHIFT))
368 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
369 warning_error(wi, warn_buf);
372 if (EEL_USER(ir->coulombtype))
374 sprintf(warn_buf, "Coulomb type %s is not supported with the verlet scheme", eel_names[ir->coulombtype]);
375 warning_error(wi, warn_buf);
378 if (ir->nstlist <= 0)
380 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
383 if (ir->nstlist < 10)
385 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.");
388 rc_max = std::max(ir->rvdw, ir->rcoulomb);
390 if (ir->verletbuf_tol <= 0)
392 if (ir->verletbuf_tol == 0)
394 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
397 if (ir->rlist < rc_max)
399 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
402 if (ir->rlist == rc_max && ir->nstlist > 1)
404 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.");
409 if (ir->rlist > rc_max)
411 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.");
414 if (ir->nstlist == 1)
416 /* No buffer required */
421 if (EI_DYNAMICS(ir->eI))
423 if (inputrec2nboundeddim(ir) < 3)
425 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.");
427 /* Set rlist temporarily so we can continue processing */
432 /* Set the buffer to 5% of the cut-off */
433 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
439 /* GENERAL INTEGRATOR STUFF */
442 if (ir->etc != etcNO)
444 if (EI_RANDOM(ir->eI))
446 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]);
450 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
452 warning_note(wi, warn_buf);
456 if (ir->eI == eiVVAK)
458 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]);
459 warning_note(wi, warn_buf);
461 if (!EI_DYNAMICS(ir->eI))
463 if (ir->epc != epcNO)
465 sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
466 warning_note(wi, warn_buf);
470 if (EI_DYNAMICS(ir->eI))
472 if (ir->nstcalcenergy < 0)
474 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
475 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
477 /* nstcalcenergy larger than nstener does not make sense.
478 * We ideally want nstcalcenergy=nstener.
482 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
486 ir->nstcalcenergy = ir->nstenergy;
490 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
491 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
492 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
495 const char *nsten = "nstenergy";
496 const char *nstdh = "nstdhdl";
497 const char *min_name = nsten;
498 int min_nst = ir->nstenergy;
500 /* find the smallest of ( nstenergy, nstdhdl ) */
501 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
502 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
504 min_nst = ir->fepvals->nstdhdl;
507 /* If the user sets nstenergy small, we should respect that */
509 "Setting nstcalcenergy (%d) equal to %s (%d)",
510 ir->nstcalcenergy, min_name, min_nst);
511 warning_note(wi, warn_buf);
512 ir->nstcalcenergy = min_nst;
515 if (ir->epc != epcNO)
517 if (ir->nstpcouple < 0)
519 ir->nstpcouple = ir_optimal_nstpcouple(ir);
523 if (ir->nstcalcenergy > 0)
525 if (ir->efep != efepNO)
527 /* nstdhdl should be a multiple of nstcalcenergy */
528 check_nst("nstcalcenergy", ir->nstcalcenergy,
529 "nstdhdl", &ir->fepvals->nstdhdl, wi);
530 /* nstexpanded should be a multiple of nstcalcenergy */
531 check_nst("nstcalcenergy", ir->nstcalcenergy,
532 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
534 /* for storing exact averages nstenergy should be
535 * a multiple of nstcalcenergy
537 check_nst("nstcalcenergy", ir->nstcalcenergy,
538 "nstenergy", &ir->nstenergy, wi);
542 if (ir->nsteps == 0 && !ir->bContinuation)
544 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
548 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
549 ir->bContinuation && ir->ld_seed != -1)
551 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)");
557 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
558 CHECK(ir->ePBC != epbcXYZ);
559 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
560 CHECK(ir->ns_type != ensGRID);
561 sprintf(err_buf, "with TPI nstlist should be larger than zero");
562 CHECK(ir->nstlist <= 0);
563 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
564 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
565 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
566 CHECK(ir->cutoff_scheme == ecutsVERLET);
570 if ( (opts->nshake > 0) && (opts->bMorse) )
573 "Using morse bond-potentials while constraining bonds is useless");
574 warning(wi, warn_buf);
577 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
578 ir->bContinuation && ir->ld_seed != -1)
580 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)");
582 /* verify simulated tempering options */
586 bool bAllTempZero = TRUE;
587 for (i = 0; i < fep->n_lambda; i++)
589 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]);
590 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
591 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
593 bAllTempZero = FALSE;
596 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
597 CHECK(bAllTempZero == TRUE);
599 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
600 CHECK(ir->eI != eiVV);
602 /* check compatability of the temperature coupling with simulated tempering */
604 if (ir->etc == etcNOSEHOOVER)
606 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
607 warning_note(wi, warn_buf);
610 /* check that the temperatures make sense */
612 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);
613 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
615 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
616 CHECK(ir->simtempvals->simtemp_high <= 0);
618 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
619 CHECK(ir->simtempvals->simtemp_low <= 0);
622 /* verify free energy options */
624 if (ir->efep != efepNO)
627 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
629 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
631 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
632 static_cast<int>(fep->sc_r_power));
633 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
635 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
636 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
638 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
639 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
641 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
642 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
644 sprintf(err_buf, "Free-energy not implemented for Ewald");
645 CHECK(ir->coulombtype == eelEWALD);
647 /* check validty of lambda inputs */
648 if (fep->n_lambda == 0)
650 /* Clear output in case of no states:*/
651 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
652 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
656 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
657 CHECK((fep->init_fep_state >= fep->n_lambda));
660 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
661 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
663 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",
664 fep->init_lambda, fep->init_fep_state);
665 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
669 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
673 for (i = 0; i < efptNR; i++)
675 if (fep->separate_dvdl[i])
680 if (n_lambda_terms > 1)
682 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.");
683 warning(wi, warn_buf);
686 if (n_lambda_terms < 2 && fep->n_lambda > 0)
689 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
693 for (j = 0; j < efptNR; j++)
695 for (i = 0; i < fep->n_lambda; i++)
697 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]);
698 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
702 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
704 for (i = 0; i < fep->n_lambda; i++)
706 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],
707 fep->all_lambda[efptCOUL][i]);
708 CHECK((fep->sc_alpha > 0) &&
709 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
710 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
711 ((fep->all_lambda[efptVDW][i] > 0.0) &&
712 (fep->all_lambda[efptVDW][i] < 1.0))));
716 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
718 real sigma, lambda, r_sc;
721 /* Maximum estimate for A and B charges equal with lambda power 1 */
723 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);
724 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.",
726 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
727 warning_note(wi, warn_buf);
730 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
731 be treated differently, but that's the next step */
733 for (i = 0; i < efptNR; i++)
735 for (j = 0; j < fep->n_lambda; j++)
737 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
738 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
743 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
747 /* checking equilibration of weights inputs for validity */
749 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
750 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
751 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
753 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
754 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
755 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
757 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
758 expand->equil_steps, elmceq_names[elmceqSTEPS]);
759 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
761 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
762 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
763 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
765 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
766 expand->equil_ratio, elmceq_names[elmceqRATIO]);
767 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
769 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
770 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
771 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
773 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
774 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
775 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
777 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
778 expand->equil_steps, elmceq_names[elmceqSTEPS]);
779 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
781 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
782 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
783 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
785 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
786 expand->equil_ratio, elmceq_names[elmceqRATIO]);
787 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
789 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
790 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
791 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
793 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
794 CHECK((expand->lmc_repeats <= 0));
795 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
796 CHECK((expand->minvarmin <= 0));
797 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
798 CHECK((expand->c_range < 0));
799 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
800 fep->init_fep_state, expand->lmc_forced_nstart);
801 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
802 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
803 CHECK((expand->lmc_forced_nstart < 0));
804 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
805 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
807 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
808 CHECK((expand->init_wl_delta < 0));
809 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
810 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
811 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
812 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
814 /* if there is no temperature control, we need to specify an MC temperature */
815 if (!integratorHasReferenceTemperature(ir) && (expand->elmcmove != elmcmoveNO) && (expand->mc_temp <= 0.0))
817 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!='no', mc_temp must be set to a positive number");
818 warning_error(wi, err_buf);
820 if (expand->nstTij > 0)
822 sprintf(err_buf, "nstlog must be non-zero");
823 CHECK(ir->nstlog == 0);
824 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
825 expand->nstTij, ir->nstlog);
826 CHECK((expand->nstTij % ir->nstlog) != 0);
831 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
832 CHECK(ir->nwall && ir->ePBC != epbcXY);
835 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
837 if (ir->ePBC == epbcNONE)
839 if (ir->epc != epcNO)
841 warning(wi, "Turning off pressure coupling for vacuum system");
847 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
848 epbc_names[ir->ePBC]);
849 CHECK(ir->epc != epcNO);
851 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
852 CHECK(EEL_FULL(ir->coulombtype));
854 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
855 epbc_names[ir->ePBC]);
856 CHECK(ir->eDispCorr != edispcNO);
859 if (ir->rlist == 0.0)
861 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
862 "with coulombtype = %s or coulombtype = %s\n"
863 "without periodic boundary conditions (pbc = %s) and\n"
864 "rcoulomb and rvdw set to zero",
865 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
866 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
867 (ir->ePBC != epbcNONE) ||
868 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
872 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
877 if (ir->nstcomm == 0)
879 ir->comm_mode = ecmNO;
881 if (ir->comm_mode != ecmNO)
885 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");
886 ir->nstcomm = abs(ir->nstcomm);
889 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
891 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
892 ir->nstcomm = ir->nstcalcenergy;
895 if (ir->comm_mode == ecmANGULAR)
897 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
898 CHECK(ir->bPeriodicMols);
899 if (ir->ePBC != epbcNONE)
901 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.");
906 if (EI_STATE_VELOCITY(ir->eI) && !EI_SD(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
908 sprintf(warn_buf, "Tumbling and flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR or use integrator = %s.", ei_names[eiSD1]);
909 warning_note(wi, warn_buf);
912 /* TEMPERATURE COUPLING */
913 if (ir->etc == etcYES)
915 ir->etc = etcBERENDSEN;
916 warning_note(wi, "Old option for temperature coupling given: "
917 "changing \"yes\" to \"Berendsen\"\n");
920 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
922 if (ir->opts.nhchainlength < 1)
924 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
925 ir->opts.nhchainlength = 1;
926 warning(wi, warn_buf);
929 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
931 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
932 ir->opts.nhchainlength = 1;
937 ir->opts.nhchainlength = 0;
940 if (ir->eI == eiVVAK)
942 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
944 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
947 if (ETC_ANDERSEN(ir->etc))
949 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
950 CHECK(!(EI_VV(ir->eI)));
952 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
954 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]);
955 warning_note(wi, warn_buf);
958 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]);
959 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
962 if (ir->etc == etcBERENDSEN)
964 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
965 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
966 warning_note(wi, warn_buf);
969 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
970 && ir->epc == epcBERENDSEN)
972 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
973 "true ensemble for the thermostat");
974 warning(wi, warn_buf);
977 /* PRESSURE COUPLING */
978 if (ir->epc == epcISOTROPIC)
980 ir->epc = epcBERENDSEN;
981 warning_note(wi, "Old option for pressure coupling given: "
982 "changing \"Isotropic\" to \"Berendsen\"\n");
985 if (ir->epc != epcNO)
987 dt_pcoupl = ir->nstpcouple*ir->delta_t;
989 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
990 CHECK(ir->tau_p <= 0);
992 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
994 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
995 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
996 warning(wi, warn_buf);
999 sprintf(err_buf, "compressibility must be > 0 when using pressure"
1000 " coupling %s\n", EPCOUPLTYPE(ir->epc));
1001 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
1002 ir->compress[ZZ][ZZ] < 0 ||
1003 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
1004 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1006 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1009 "You are generating velocities so I am assuming you "
1010 "are equilibrating a system. You are using "
1011 "%s pressure coupling, but this can be "
1012 "unstable for equilibration. If your system crashes, try "
1013 "equilibrating first with Berendsen pressure coupling. If "
1014 "you are not equilibrating the system, you can probably "
1015 "ignore this warning.",
1016 epcoupl_names[ir->epc]);
1017 warning(wi, warn_buf);
1023 if (ir->epc > epcNO)
1025 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1027 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.");
1033 if (ir->epc == epcMTTK)
1035 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1039 /* ELECTROSTATICS */
1040 /* More checks are in triple check (grompp.c) */
1042 if (ir->coulombtype == eelSWITCH)
1044 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1045 "artifacts, advice: use coulombtype = %s",
1046 eel_names[ir->coulombtype],
1047 eel_names[eelRF_ZERO]);
1048 warning(wi, warn_buf);
1051 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1053 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);
1054 warning(wi, warn_buf);
1055 ir->epsilon_rf = ir->epsilon_r;
1056 ir->epsilon_r = 1.0;
1059 if (ir->epsilon_r == 0)
1062 "It is pointless to use long-range electrostatics with infinite relative permittivity."
1063 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1064 CHECK(EEL_FULL(ir->coulombtype));
1067 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == nullptr)
1069 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1070 CHECK(ir->epsilon_r < 0);
1073 if (EEL_RF(ir->coulombtype))
1075 /* reaction field (at the cut-off) */
1077 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1079 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1080 eel_names[ir->coulombtype]);
1081 warning(wi, warn_buf);
1082 ir->epsilon_rf = 0.0;
1085 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1086 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1087 (ir->epsilon_r == 0));
1088 if (ir->epsilon_rf == ir->epsilon_r)
1090 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1091 eel_names[ir->coulombtype]);
1092 warning(wi, warn_buf);
1095 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1096 * means the interaction is zero outside rcoulomb, but it helps to
1097 * provide accurate energy conservation.
1099 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1101 if (ir_coulomb_switched(ir))
1104 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1105 eel_names[ir->coulombtype]);
1106 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1109 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1111 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1113 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1114 eel_names[ir->coulombtype]);
1115 CHECK(ir->rlist > ir->rcoulomb);
1119 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1122 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1123 CHECK( ir->coulomb_modifier != eintmodNONE);
1125 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1128 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1129 CHECK( ir->vdw_modifier != eintmodNONE);
1132 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1133 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1136 "The switch/shift interaction settings are just for compatibility; you will get better "
1137 "performance from applying potential modifiers to your interactions!\n");
1138 warning_note(wi, warn_buf);
1141 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1143 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1145 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1146 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.",
1147 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1148 warning(wi, warn_buf);
1152 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1154 if (ir->rvdw_switch == 0)
1156 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.");
1157 warning(wi, warn_buf);
1161 if (EEL_FULL(ir->coulombtype))
1163 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1164 ir->coulombtype == eelPMEUSERSWITCH)
1166 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1167 eel_names[ir->coulombtype]);
1168 CHECK(ir->rcoulomb > ir->rlist);
1170 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1172 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1175 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist.\n"
1176 "For optimal energy conservation,consider using\n"
1177 "a potential modifier.", eel_names[ir->coulombtype]);
1178 CHECK(ir->rcoulomb != ir->rlist);
1183 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1185 // TODO: Move these checks into the ewald module with the options class
1187 int orderMax = (ir->coulombtype == eelP3M_AD ? 8 : 12);
1189 if (ir->pme_order < orderMin || ir->pme_order > orderMax)
1191 sprintf(warn_buf, "With coulombtype = %s, you should have %d <= pme-order <= %d", eel_names[ir->coulombtype], orderMin, orderMax);
1192 warning_error(wi, warn_buf);
1196 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1198 if (ir->ewald_geometry == eewg3D)
1200 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1201 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1202 warning(wi, warn_buf);
1204 /* This check avoids extra pbc coding for exclusion corrections */
1205 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1206 CHECK(ir->wall_ewald_zfac < 2);
1208 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1209 EEL_FULL(ir->coulombtype))
1211 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1212 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1213 warning(wi, warn_buf);
1215 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1217 if (ir->cutoff_scheme == ecutsVERLET)
1219 sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
1220 eel_names[ir->coulombtype]);
1221 warning(wi, warn_buf);
1225 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",
1226 eel_names[ir->coulombtype]);
1227 warning_note(wi, warn_buf);
1231 if (ir_vdw_switched(ir))
1233 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1234 CHECK(ir->rvdw_switch >= ir->rvdw);
1236 if (ir->rvdw_switch < 0.5*ir->rvdw)
1238 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.",
1239 ir->rvdw_switch, ir->rvdw);
1240 warning_note(wi, warn_buf);
1243 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1245 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1247 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1248 CHECK(ir->rlist > ir->rvdw);
1252 if (ir->vdwtype == evdwPME)
1254 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1256 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1257 evdw_names[ir->vdwtype],
1258 eintmod_names[eintmodPOTSHIFT],
1259 eintmod_names[eintmodNONE]);
1260 warning_error(wi, err_buf);
1264 if (ir->cutoff_scheme == ecutsGROUP)
1266 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1267 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)))
1269 warning_note(wi, "With exact cut-offs, rlist should be "
1270 "larger than rcoulomb and rvdw, so that there "
1271 "is a buffer region for particle motion "
1272 "between neighborsearch steps");
1275 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlist <= ir->rcoulomb)
1277 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rcoulomb.");
1278 warning_note(wi, warn_buf);
1280 if (ir_vdw_switched(ir) && (ir->rlist <= ir->rvdw))
1282 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rvdw.");
1283 warning_note(wi, warn_buf);
1287 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1289 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.");
1292 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1295 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1298 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1300 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1303 /* ENERGY CONSERVATION */
1304 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1306 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1308 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1309 evdw_names[evdwSHIFT]);
1310 warning_note(wi, warn_buf);
1312 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1314 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1315 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1316 warning_note(wi, warn_buf);
1320 /* IMPLICIT SOLVENT */
1321 if (ir->coulombtype == eelGB_NOTUSED)
1323 sprintf(warn_buf, "Invalid option %s for coulombtype",
1324 eel_names[ir->coulombtype]);
1325 warning_error(wi, warn_buf);
1330 if (ir->cutoff_scheme != ecutsGROUP)
1332 warning_error(wi, "QMMM is currently only supported with cutoff-scheme=group");
1334 if (!EI_DYNAMICS(ir->eI))
1337 sprintf(buf, "QMMM is only supported with dynamics, not with integrator %s", ei_names[ir->eI]);
1338 warning_error(wi, buf);
1344 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1348 /* count the number of text elemets separated by whitespace in a string.
1349 str = the input string
1350 maxptr = the maximum number of allowed elements
1351 ptr = the output array of pointers to the first character of each element
1352 returns: the number of elements. */
1353 int str_nelem(const char *str, int maxptr, char *ptr[])
1358 copy0 = gmx_strdup(str);
1361 while (*copy != '\0')
1365 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1373 while ((*copy != '\0') && !isspace(*copy))
1392 /* interpret a number of doubles from a string and put them in an array,
1393 after allocating space for them.
1394 str = the input string
1395 n = the (pre-allocated) number of doubles read
1396 r = the output array of doubles. */
1397 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1402 char warn_buf[STRLEN];
1404 *n = str_nelem(str, MAXPTR, ptr);
1407 for (i = 0; i < *n; i++)
1409 (*r)[i] = strtod(ptr[i], &endptr);
1412 sprintf(warn_buf, "Invalid value %s in string in mdp file. Expected a real number.", ptr[i]);
1413 warning_error(wi, warn_buf);
1418 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1421 int i, j, max_n_lambda, nweights, nfep[efptNR];
1422 t_lambda *fep = ir->fepvals;
1423 t_expanded *expand = ir->expandedvals;
1424 real **count_fep_lambdas;
1425 bool bOneLambda = TRUE;
1427 snew(count_fep_lambdas, efptNR);
1429 /* FEP input processing */
1430 /* first, identify the number of lambda values for each type.
1431 All that are nonzero must have the same number */
1433 for (i = 0; i < efptNR; i++)
1435 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1438 /* now, determine the number of components. All must be either zero, or equal. */
1441 for (i = 0; i < efptNR; i++)
1443 if (nfep[i] > max_n_lambda)
1445 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1446 must have the same number if its not zero.*/
1451 for (i = 0; i < efptNR; i++)
1455 ir->fepvals->separate_dvdl[i] = FALSE;
1457 else if (nfep[i] == max_n_lambda)
1459 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1460 respect to the temperature currently */
1462 ir->fepvals->separate_dvdl[i] = TRUE;
1467 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1468 nfep[i], efpt_names[i], max_n_lambda);
1471 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1472 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1474 /* the number of lambdas is the number we've read in, which is either zero
1475 or the same for all */
1476 fep->n_lambda = max_n_lambda;
1478 /* allocate space for the array of lambda values */
1479 snew(fep->all_lambda, efptNR);
1480 /* if init_lambda is defined, we need to set lambda */
1481 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1483 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1485 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1486 for (i = 0; i < efptNR; i++)
1488 snew(fep->all_lambda[i], fep->n_lambda);
1489 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1492 for (j = 0; j < fep->n_lambda; j++)
1494 fep->all_lambda[i][j] = static_cast<double>(count_fep_lambdas[i][j]);
1496 sfree(count_fep_lambdas[i]);
1499 sfree(count_fep_lambdas);
1501 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1502 bookkeeping -- for now, init_lambda */
1504 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1506 for (i = 0; i < fep->n_lambda; i++)
1508 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1512 /* check to see if only a single component lambda is defined, and soft core is defined.
1513 In this case, turn on coulomb soft core */
1515 if (max_n_lambda == 0)
1521 for (i = 0; i < efptNR; i++)
1523 if ((nfep[i] != 0) && (i != efptFEP))
1529 if ((bOneLambda) && (fep->sc_alpha > 0))
1531 fep->bScCoul = TRUE;
1534 /* Fill in the others with the efptFEP if they are not explicitly
1535 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1536 they are all zero. */
1538 for (i = 0; i < efptNR; i++)
1540 if ((nfep[i] == 0) && (i != efptFEP))
1542 for (j = 0; j < fep->n_lambda; j++)
1544 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1550 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1551 if (fep->sc_r_power == 48)
1553 if (fep->sc_alpha > 0.1)
1555 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1559 /* now read in the weights */
1560 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1563 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1565 else if (nweights != fep->n_lambda)
1567 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1568 nweights, fep->n_lambda);
1570 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1572 expand->nstexpanded = fep->nstdhdl;
1573 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1575 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1577 expand->nstexpanded = 2*static_cast<int>(ir->opts.tau_t[0]/ir->delta_t);
1578 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1579 2*tau_t just to be careful so it's not to frequent */
1584 static void do_simtemp_params(t_inputrec *ir)
1587 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1588 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1591 static void do_wall_params(t_inputrec *ir,
1592 char *wall_atomtype, char *wall_density,
1596 char *names[MAXPTR];
1599 opts->wall_atomtype[0] = nullptr;
1600 opts->wall_atomtype[1] = nullptr;
1602 ir->wall_atomtype[0] = -1;
1603 ir->wall_atomtype[1] = -1;
1604 ir->wall_density[0] = 0;
1605 ir->wall_density[1] = 0;
1609 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1610 if (nstr != ir->nwall)
1612 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1615 for (i = 0; i < ir->nwall; i++)
1617 opts->wall_atomtype[i] = gmx_strdup(names[i]);
1620 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1622 nstr = str_nelem(wall_density, MAXPTR, names);
1623 if (nstr != ir->nwall)
1625 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1627 for (i = 0; i < ir->nwall; i++)
1629 if (sscanf(names[i], "%lf", &dbl) != 1)
1631 gmx_fatal(FARGS, "Could not parse wall-density value from string '%s'", names[i]);
1635 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1637 ir->wall_density[i] = dbl;
1643 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1651 srenew(groups->grpname, groups->ngrpname+nwall);
1652 grps = &(groups->grps[egcENER]);
1653 srenew(grps->nm_ind, grps->nr+nwall);
1654 for (i = 0; i < nwall; i++)
1656 sprintf(str, "wall%d", i);
1657 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1658 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1663 static void read_expandedparams(std::vector<t_inpfile> *inp,
1664 t_expanded *expand, warninp_t wi)
1666 /* read expanded ensemble parameters */
1667 printStringNewline(inp, "expanded ensemble variables");
1668 expand->nstexpanded = get_eint(inp, "nstexpanded", -1, wi);
1669 expand->elamstats = get_eeenum(inp, "lmc-stats", elamstats_names, wi);
1670 expand->elmcmove = get_eeenum(inp, "lmc-move", elmcmove_names, wi);
1671 expand->elmceq = get_eeenum(inp, "lmc-weights-equil", elmceq_names, wi);
1672 expand->equil_n_at_lam = get_eint(inp, "weight-equil-number-all-lambda", -1, wi);
1673 expand->equil_samples = get_eint(inp, "weight-equil-number-samples", -1, wi);
1674 expand->equil_steps = get_eint(inp, "weight-equil-number-steps", -1, wi);
1675 expand->equil_wl_delta = get_ereal(inp, "weight-equil-wl-delta", -1, wi);
1676 expand->equil_ratio = get_ereal(inp, "weight-equil-count-ratio", -1, wi);
1677 printStringNewline(inp, "Seed for Monte Carlo in lambda space");
1678 expand->lmc_seed = get_eint(inp, "lmc-seed", -1, wi);
1679 expand->mc_temp = get_ereal(inp, "mc-temperature", -1, wi);
1680 expand->lmc_repeats = get_eint(inp, "lmc-repeats", 1, wi);
1681 expand->gibbsdeltalam = get_eint(inp, "lmc-gibbsdelta", -1, wi);
1682 expand->lmc_forced_nstart = get_eint(inp, "lmc-forced-nstart", 0, wi);
1683 expand->bSymmetrizedTMatrix = (get_eeenum(inp, "symmetrized-transition-matrix", yesno_names, wi) != 0);
1684 expand->nstTij = get_eint(inp, "nst-transition-matrix", -1, wi);
1685 expand->minvarmin = get_eint(inp, "mininum-var-min", 100, wi); /*default is reasonable */
1686 expand->c_range = get_eint(inp, "weight-c-range", 0, wi); /* default is just C=0 */
1687 expand->wl_scale = get_ereal(inp, "wl-scale", 0.8, wi);
1688 expand->wl_ratio = get_ereal(inp, "wl-ratio", 0.8, wi);
1689 expand->init_wl_delta = get_ereal(inp, "init-wl-delta", 1.0, wi);
1690 expand->bWLoneovert = (get_eeenum(inp, "wl-oneovert", yesno_names, wi) != 0);
1693 /*! \brief Return whether an end state with the given coupling-lambda
1694 * value describes fully-interacting VDW.
1696 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1697 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1699 static bool couple_lambda_has_vdw_on(int couple_lambda_value)
1701 return (couple_lambda_value == ecouplamVDW ||
1702 couple_lambda_value == ecouplamVDWQ);
1708 class MdpErrorHandler : public gmx::IKeyValueTreeErrorHandler
1711 explicit MdpErrorHandler(warninp_t wi)
1712 : wi_(wi), mapping_(nullptr)
1716 void setBackMapping(const gmx::IKeyValueTreeBackMapping &mapping)
1718 mapping_ = &mapping;
1721 virtual bool onError(gmx::UserInputError *ex, const gmx::KeyValueTreePath &context)
1723 ex->prependContext(gmx::formatString("Error in mdp option \"%s\":",
1724 getOptionName(context).c_str()));
1725 std::string message = gmx::formatExceptionMessageToString(*ex);
1726 warning_error(wi_, message.c_str());
1731 std::string getOptionName(const gmx::KeyValueTreePath &context)
1733 if (mapping_ != nullptr)
1735 gmx::KeyValueTreePath path = mapping_->originalPath(context);
1736 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
1739 GMX_ASSERT(context.size() == 1, "Inconsistent context for mdp option parsing");
1744 const gmx::IKeyValueTreeBackMapping *mapping_;
1749 void get_ir(const char *mdparin, const char *mdparout,
1750 gmx::MDModules *mdModules, t_inputrec *ir, t_gromppopts *opts,
1751 WriteMdpHeader writeMdpHeader, warninp_t wi)
1754 double dumdub[2][6];
1756 char warn_buf[STRLEN];
1757 t_lambda *fep = ir->fepvals;
1758 t_expanded *expand = ir->expandedvals;
1760 const char *no_names[] = { "no", nullptr };
1762 init_inputrec_strings();
1763 gmx::TextInputFile stream(mdparin);
1764 std::vector<t_inpfile> inp = read_inpfile(&stream, mdparin, wi);
1766 snew(dumstr[0], STRLEN);
1767 snew(dumstr[1], STRLEN);
1769 if (-1 == search_einp(inp, "cutoff-scheme"))
1772 "%s did not specify a value for the .mdp option "
1773 "\"cutoff-scheme\". Probably it was first intended for use "
1774 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1775 "introduced, but the group scheme was still the default. "
1776 "The default is now the Verlet scheme, so you will observe "
1777 "different behaviour.", mdparin);
1778 warning_note(wi, warn_buf);
1781 /* ignore the following deprecated commands */
1782 replace_inp_entry(inp, "title", nullptr);
1783 replace_inp_entry(inp, "cpp", nullptr);
1784 replace_inp_entry(inp, "domain-decomposition", nullptr);
1785 replace_inp_entry(inp, "andersen-seed", nullptr);
1786 replace_inp_entry(inp, "dihre", nullptr);
1787 replace_inp_entry(inp, "dihre-fc", nullptr);
1788 replace_inp_entry(inp, "dihre-tau", nullptr);
1789 replace_inp_entry(inp, "nstdihreout", nullptr);
1790 replace_inp_entry(inp, "nstcheckpoint", nullptr);
1791 replace_inp_entry(inp, "optimize-fft", nullptr);
1792 replace_inp_entry(inp, "adress_type", nullptr);
1793 replace_inp_entry(inp, "adress_const_wf", nullptr);
1794 replace_inp_entry(inp, "adress_ex_width", nullptr);
1795 replace_inp_entry(inp, "adress_hy_width", nullptr);
1796 replace_inp_entry(inp, "adress_ex_forcecap", nullptr);
1797 replace_inp_entry(inp, "adress_interface_correction", nullptr);
1798 replace_inp_entry(inp, "adress_site", nullptr);
1799 replace_inp_entry(inp, "adress_reference_coords", nullptr);
1800 replace_inp_entry(inp, "adress_tf_grp_names", nullptr);
1801 replace_inp_entry(inp, "adress_cg_grp_names", nullptr);
1802 replace_inp_entry(inp, "adress_do_hybridpairs", nullptr);
1803 replace_inp_entry(inp, "rlistlong", nullptr);
1804 replace_inp_entry(inp, "nstcalclr", nullptr);
1805 replace_inp_entry(inp, "pull-print-com2", nullptr);
1806 replace_inp_entry(inp, "gb-algorithm", nullptr);
1807 replace_inp_entry(inp, "nstgbradii", nullptr);
1808 replace_inp_entry(inp, "rgbradii", nullptr);
1809 replace_inp_entry(inp, "gb-epsilon-solvent", nullptr);
1810 replace_inp_entry(inp, "gb-saltconc", nullptr);
1811 replace_inp_entry(inp, "gb-obc-alpha", nullptr);
1812 replace_inp_entry(inp, "gb-obc-beta", nullptr);
1813 replace_inp_entry(inp, "gb-obc-gamma", nullptr);
1814 replace_inp_entry(inp, "gb-dielectric-offset", nullptr);
1815 replace_inp_entry(inp, "sa-algorithm", nullptr);
1816 replace_inp_entry(inp, "sa-surface-tension", nullptr);
1818 /* replace the following commands with the clearer new versions*/
1819 replace_inp_entry(inp, "unconstrained-start", "continuation");
1820 replace_inp_entry(inp, "foreign-lambda", "fep-lambdas");
1821 replace_inp_entry(inp, "verlet-buffer-drift", "verlet-buffer-tolerance");
1822 replace_inp_entry(inp, "nstxtcout", "nstxout-compressed");
1823 replace_inp_entry(inp, "xtc-grps", "compressed-x-grps");
1824 replace_inp_entry(inp, "xtc-precision", "compressed-x-precision");
1825 replace_inp_entry(inp, "pull-print-com1", "pull-print-com");
1827 printStringNewline(&inp, "VARIOUS PREPROCESSING OPTIONS");
1828 printStringNoNewline(&inp, "Preprocessor information: use cpp syntax.");
1829 printStringNoNewline(&inp, "e.g.: -I/home/joe/doe -I/home/mary/roe");
1830 setStringEntry(&inp, "include", opts->include, nullptr);
1831 printStringNoNewline(&inp, "e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1832 setStringEntry(&inp, "define", opts->define, nullptr);
1834 printStringNewline(&inp, "RUN CONTROL PARAMETERS");
1835 ir->eI = get_eeenum(&inp, "integrator", ei_names, wi);
1836 printStringNoNewline(&inp, "Start time and timestep in ps");
1837 ir->init_t = get_ereal(&inp, "tinit", 0.0, wi);
1838 ir->delta_t = get_ereal(&inp, "dt", 0.001, wi);
1839 ir->nsteps = get_eint64(&inp, "nsteps", 0, wi);
1840 printStringNoNewline(&inp, "For exact run continuation or redoing part of a run");
1841 ir->init_step = get_eint64(&inp, "init-step", 0, wi);
1842 printStringNoNewline(&inp, "Part index is updated automatically on checkpointing (keeps files separate)");
1843 ir->simulation_part = get_eint(&inp, "simulation-part", 1, wi);
1844 printStringNoNewline(&inp, "mode for center of mass motion removal");
1845 ir->comm_mode = get_eeenum(&inp, "comm-mode", ecm_names, wi);
1846 printStringNoNewline(&inp, "number of steps for center of mass motion removal");
1847 ir->nstcomm = get_eint(&inp, "nstcomm", 100, wi);
1848 printStringNoNewline(&inp, "group(s) for center of mass motion removal");
1849 setStringEntry(&inp, "comm-grps", is->vcm, nullptr);
1851 printStringNewline(&inp, "LANGEVIN DYNAMICS OPTIONS");
1852 printStringNoNewline(&inp, "Friction coefficient (amu/ps) and random seed");
1853 ir->bd_fric = get_ereal(&inp, "bd-fric", 0.0, wi);
1854 ir->ld_seed = get_eint64(&inp, "ld-seed", -1, wi);
1857 printStringNewline(&inp, "ENERGY MINIMIZATION OPTIONS");
1858 printStringNoNewline(&inp, "Force tolerance and initial step-size");
1859 ir->em_tol = get_ereal(&inp, "emtol", 10.0, wi);
1860 ir->em_stepsize = get_ereal(&inp, "emstep", 0.01, wi);
1861 printStringNoNewline(&inp, "Max number of iterations in relax-shells");
1862 ir->niter = get_eint(&inp, "niter", 20, wi);
1863 printStringNoNewline(&inp, "Step size (ps^2) for minimization of flexible constraints");
1864 ir->fc_stepsize = get_ereal(&inp, "fcstep", 0, wi);
1865 printStringNoNewline(&inp, "Frequency of steepest descents steps when doing CG");
1866 ir->nstcgsteep = get_eint(&inp, "nstcgsteep", 1000, wi);
1867 ir->nbfgscorr = get_eint(&inp, "nbfgscorr", 10, wi);
1869 printStringNewline(&inp, "TEST PARTICLE INSERTION OPTIONS");
1870 ir->rtpi = get_ereal(&inp, "rtpi", 0.05, wi);
1872 /* Output options */
1873 printStringNewline(&inp, "OUTPUT CONTROL OPTIONS");
1874 printStringNoNewline(&inp, "Output frequency for coords (x), velocities (v) and forces (f)");
1875 ir->nstxout = get_eint(&inp, "nstxout", 0, wi);
1876 ir->nstvout = get_eint(&inp, "nstvout", 0, wi);
1877 ir->nstfout = get_eint(&inp, "nstfout", 0, wi);
1878 printStringNoNewline(&inp, "Output frequency for energies to log file and energy file");
1879 ir->nstlog = get_eint(&inp, "nstlog", 1000, wi);
1880 ir->nstcalcenergy = get_eint(&inp, "nstcalcenergy", 100, wi);
1881 ir->nstenergy = get_eint(&inp, "nstenergy", 1000, wi);
1882 printStringNoNewline(&inp, "Output frequency and precision for .xtc file");
1883 ir->nstxout_compressed = get_eint(&inp, "nstxout-compressed", 0, wi);
1884 ir->x_compression_precision = get_ereal(&inp, "compressed-x-precision", 1000.0, wi);
1885 printStringNoNewline(&inp, "This selects the subset of atoms for the compressed");
1886 printStringNoNewline(&inp, "trajectory file. You can select multiple groups. By");
1887 printStringNoNewline(&inp, "default, all atoms will be written.");
1888 setStringEntry(&inp, "compressed-x-grps", is->x_compressed_groups, nullptr);
1889 printStringNoNewline(&inp, "Selection of energy groups");
1890 setStringEntry(&inp, "energygrps", is->energy, nullptr);
1892 /* Neighbor searching */
1893 printStringNewline(&inp, "NEIGHBORSEARCHING PARAMETERS");
1894 printStringNoNewline(&inp, "cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1895 ir->cutoff_scheme = get_eeenum(&inp, "cutoff-scheme", ecutscheme_names, wi);
1896 printStringNoNewline(&inp, "nblist update frequency");
1897 ir->nstlist = get_eint(&inp, "nstlist", 10, wi);
1898 printStringNoNewline(&inp, "ns algorithm (simple or grid)");
1899 ir->ns_type = get_eeenum(&inp, "ns-type", ens_names, wi);
1900 printStringNoNewline(&inp, "Periodic boundary conditions: xyz, no, xy");
1901 ir->ePBC = get_eeenum(&inp, "pbc", epbc_names, wi);
1902 ir->bPeriodicMols = get_eeenum(&inp, "periodic-molecules", yesno_names, wi) != 0;
1903 printStringNoNewline(&inp, "Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1904 printStringNoNewline(&inp, "a value of -1 means: use rlist");
1905 ir->verletbuf_tol = get_ereal(&inp, "verlet-buffer-tolerance", 0.005, wi);
1906 printStringNoNewline(&inp, "nblist cut-off");
1907 ir->rlist = get_ereal(&inp, "rlist", 1.0, wi);
1908 printStringNoNewline(&inp, "long-range cut-off for switched potentials");
1910 /* Electrostatics */
1911 printStringNewline(&inp, "OPTIONS FOR ELECTROSTATICS AND VDW");
1912 printStringNoNewline(&inp, "Method for doing electrostatics");
1913 ir->coulombtype = get_eeenum(&inp, "coulombtype", eel_names, wi);
1914 ir->coulomb_modifier = get_eeenum(&inp, "coulomb-modifier", eintmod_names, wi);
1915 printStringNoNewline(&inp, "cut-off lengths");
1916 ir->rcoulomb_switch = get_ereal(&inp, "rcoulomb-switch", 0.0, wi);
1917 ir->rcoulomb = get_ereal(&inp, "rcoulomb", 1.0, wi);
1918 printStringNoNewline(&inp, "Relative dielectric constant for the medium and the reaction field");
1919 ir->epsilon_r = get_ereal(&inp, "epsilon-r", 1.0, wi);
1920 ir->epsilon_rf = get_ereal(&inp, "epsilon-rf", 0.0, wi);
1921 printStringNoNewline(&inp, "Method for doing Van der Waals");
1922 ir->vdwtype = get_eeenum(&inp, "vdw-type", evdw_names, wi);
1923 ir->vdw_modifier = get_eeenum(&inp, "vdw-modifier", eintmod_names, wi);
1924 printStringNoNewline(&inp, "cut-off lengths");
1925 ir->rvdw_switch = get_ereal(&inp, "rvdw-switch", 0.0, wi);
1926 ir->rvdw = get_ereal(&inp, "rvdw", 1.0, wi);
1927 printStringNoNewline(&inp, "Apply long range dispersion corrections for Energy and Pressure");
1928 ir->eDispCorr = get_eeenum(&inp, "DispCorr", edispc_names, wi);
1929 printStringNoNewline(&inp, "Extension of the potential lookup tables beyond the cut-off");
1930 ir->tabext = get_ereal(&inp, "table-extension", 1.0, wi);
1931 printStringNoNewline(&inp, "Separate tables between energy group pairs");
1932 setStringEntry(&inp, "energygrp-table", is->egptable, nullptr);
1933 printStringNoNewline(&inp, "Spacing for the PME/PPPM FFT grid");
1934 ir->fourier_spacing = get_ereal(&inp, "fourierspacing", 0.12, wi);
1935 printStringNoNewline(&inp, "FFT grid size, when a value is 0 fourierspacing will be used");
1936 ir->nkx = get_eint(&inp, "fourier-nx", 0, wi);
1937 ir->nky = get_eint(&inp, "fourier-ny", 0, wi);
1938 ir->nkz = get_eint(&inp, "fourier-nz", 0, wi);
1939 printStringNoNewline(&inp, "EWALD/PME/PPPM parameters");
1940 ir->pme_order = get_eint(&inp, "pme-order", 4, wi);
1941 ir->ewald_rtol = get_ereal(&inp, "ewald-rtol", 0.00001, wi);
1942 ir->ewald_rtol_lj = get_ereal(&inp, "ewald-rtol-lj", 0.001, wi);
1943 ir->ljpme_combination_rule = get_eeenum(&inp, "lj-pme-comb-rule", eljpme_names, wi);
1944 ir->ewald_geometry = get_eeenum(&inp, "ewald-geometry", eewg_names, wi);
1945 ir->epsilon_surface = get_ereal(&inp, "epsilon-surface", 0.0, wi);
1947 /* Implicit solvation is no longer supported, but we need grompp
1948 to be able to refuse old .mdp files that would have built a tpr
1949 to run it. Thus, only "no" is accepted. */
1950 ir->implicit_solvent = (get_eeenum(&inp, "implicit-solvent", no_names, wi) != 0);
1952 /* Coupling stuff */
1953 printStringNewline(&inp, "OPTIONS FOR WEAK COUPLING ALGORITHMS");
1954 printStringNoNewline(&inp, "Temperature coupling");
1955 ir->etc = get_eeenum(&inp, "tcoupl", etcoupl_names, wi);
1956 ir->nsttcouple = get_eint(&inp, "nsttcouple", -1, wi);
1957 ir->opts.nhchainlength = get_eint(&inp, "nh-chain-length", 10, wi);
1958 ir->bPrintNHChains = (get_eeenum(&inp, "print-nose-hoover-chain-variables", yesno_names, wi) != 0);
1959 printStringNoNewline(&inp, "Groups to couple separately");
1960 setStringEntry(&inp, "tc-grps", is->tcgrps, nullptr);
1961 printStringNoNewline(&inp, "Time constant (ps) and reference temperature (K)");
1962 setStringEntry(&inp, "tau-t", is->tau_t, nullptr);
1963 setStringEntry(&inp, "ref-t", is->ref_t, nullptr);
1964 printStringNoNewline(&inp, "pressure coupling");
1965 ir->epc = get_eeenum(&inp, "pcoupl", epcoupl_names, wi);
1966 ir->epct = get_eeenum(&inp, "pcoupltype", epcoupltype_names, wi);
1967 ir->nstpcouple = get_eint(&inp, "nstpcouple", -1, wi);
1968 printStringNoNewline(&inp, "Time constant (ps), compressibility (1/bar) and reference P (bar)");
1969 ir->tau_p = get_ereal(&inp, "tau-p", 1.0, wi);
1970 setStringEntry(&inp, "compressibility", dumstr[0], nullptr);
1971 setStringEntry(&inp, "ref-p", dumstr[1], nullptr);
1972 printStringNoNewline(&inp, "Scaling of reference coordinates, No, All or COM");
1973 ir->refcoord_scaling = get_eeenum(&inp, "refcoord-scaling", erefscaling_names, wi);
1976 printStringNewline(&inp, "OPTIONS FOR QMMM calculations");
1977 ir->bQMMM = (get_eeenum(&inp, "QMMM", yesno_names, wi) != 0);
1978 printStringNoNewline(&inp, "Groups treated Quantum Mechanically");
1979 setStringEntry(&inp, "QMMM-grps", is->QMMM, nullptr);
1980 printStringNoNewline(&inp, "QM method");
1981 setStringEntry(&inp, "QMmethod", is->QMmethod, nullptr);
1982 printStringNoNewline(&inp, "QMMM scheme");
1983 ir->QMMMscheme = get_eeenum(&inp, "QMMMscheme", eQMMMscheme_names, wi);
1984 printStringNoNewline(&inp, "QM basisset");
1985 setStringEntry(&inp, "QMbasis", is->QMbasis, nullptr);
1986 printStringNoNewline(&inp, "QM charge");
1987 setStringEntry(&inp, "QMcharge", is->QMcharge, nullptr);
1988 printStringNoNewline(&inp, "QM multiplicity");
1989 setStringEntry(&inp, "QMmult", is->QMmult, nullptr);
1990 printStringNoNewline(&inp, "Surface Hopping");
1991 setStringEntry(&inp, "SH", is->bSH, nullptr);
1992 printStringNoNewline(&inp, "CAS space options");
1993 setStringEntry(&inp, "CASorbitals", is->CASorbitals, nullptr);
1994 setStringEntry(&inp, "CASelectrons", is->CASelectrons, nullptr);
1995 setStringEntry(&inp, "SAon", is->SAon, nullptr);
1996 setStringEntry(&inp, "SAoff", is->SAoff, nullptr);
1997 setStringEntry(&inp, "SAsteps", is->SAsteps, nullptr);
1998 printStringNoNewline(&inp, "Scale factor for MM charges");
1999 ir->scalefactor = get_ereal(&inp, "MMChargeScaleFactor", 1.0, wi);
2001 /* Simulated annealing */
2002 printStringNewline(&inp, "SIMULATED ANNEALING");
2003 printStringNoNewline(&inp, "Type of annealing for each temperature group (no/single/periodic)");
2004 setStringEntry(&inp, "annealing", is->anneal, nullptr);
2005 printStringNoNewline(&inp, "Number of time points to use for specifying annealing in each group");
2006 setStringEntry(&inp, "annealing-npoints", is->anneal_npoints, nullptr);
2007 printStringNoNewline(&inp, "List of times at the annealing points for each group");
2008 setStringEntry(&inp, "annealing-time", is->anneal_time, nullptr);
2009 printStringNoNewline(&inp, "Temp. at each annealing point, for each group.");
2010 setStringEntry(&inp, "annealing-temp", is->anneal_temp, nullptr);
2013 printStringNewline(&inp, "GENERATE VELOCITIES FOR STARTUP RUN");
2014 opts->bGenVel = (get_eeenum(&inp, "gen-vel", yesno_names, wi) != 0);
2015 opts->tempi = get_ereal(&inp, "gen-temp", 300.0, wi);
2016 opts->seed = get_eint(&inp, "gen-seed", -1, wi);
2019 printStringNewline(&inp, "OPTIONS FOR BONDS");
2020 opts->nshake = get_eeenum(&inp, "constraints", constraints, wi);
2021 printStringNoNewline(&inp, "Type of constraint algorithm");
2022 ir->eConstrAlg = get_eeenum(&inp, "constraint-algorithm", econstr_names, wi);
2023 printStringNoNewline(&inp, "Do not constrain the start configuration");
2024 ir->bContinuation = (get_eeenum(&inp, "continuation", yesno_names, wi) != 0);
2025 printStringNoNewline(&inp, "Use successive overrelaxation to reduce the number of shake iterations");
2026 ir->bShakeSOR = (get_eeenum(&inp, "Shake-SOR", yesno_names, wi) != 0);
2027 printStringNoNewline(&inp, "Relative tolerance of shake");
2028 ir->shake_tol = get_ereal(&inp, "shake-tol", 0.0001, wi);
2029 printStringNoNewline(&inp, "Highest order in the expansion of the constraint coupling matrix");
2030 ir->nProjOrder = get_eint(&inp, "lincs-order", 4, wi);
2031 printStringNoNewline(&inp, "Number of iterations in the final step of LINCS. 1 is fine for");
2032 printStringNoNewline(&inp, "normal simulations, but use 2 to conserve energy in NVE runs.");
2033 printStringNoNewline(&inp, "For energy minimization with constraints it should be 4 to 8.");
2034 ir->nLincsIter = get_eint(&inp, "lincs-iter", 1, wi);
2035 printStringNoNewline(&inp, "Lincs will write a warning to the stderr if in one step a bond");
2036 printStringNoNewline(&inp, "rotates over more degrees than");
2037 ir->LincsWarnAngle = get_ereal(&inp, "lincs-warnangle", 30.0, wi);
2038 printStringNoNewline(&inp, "Convert harmonic bonds to morse potentials");
2039 opts->bMorse = (get_eeenum(&inp, "morse", yesno_names, wi) != 0);
2041 /* Energy group exclusions */
2042 printStringNewline(&inp, "ENERGY GROUP EXCLUSIONS");
2043 printStringNoNewline(&inp, "Pairs of energy groups for which all non-bonded interactions are excluded");
2044 setStringEntry(&inp, "energygrp-excl", is->egpexcl, nullptr);
2047 printStringNewline(&inp, "WALLS");
2048 printStringNoNewline(&inp, "Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2049 ir->nwall = get_eint(&inp, "nwall", 0, wi);
2050 ir->wall_type = get_eeenum(&inp, "wall-type", ewt_names, wi);
2051 ir->wall_r_linpot = get_ereal(&inp, "wall-r-linpot", -1, wi);
2052 setStringEntry(&inp, "wall-atomtype", is->wall_atomtype, nullptr);
2053 setStringEntry(&inp, "wall-density", is->wall_density, nullptr);
2054 ir->wall_ewald_zfac = get_ereal(&inp, "wall-ewald-zfac", 3, wi);
2057 printStringNewline(&inp, "COM PULLING");
2058 ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
2062 is->pull_grp = read_pullparams(&inp, ir->pull, wi);
2066 NOTE: needs COM pulling input */
2067 printStringNewline(&inp, "AWH biasing");
2068 ir->bDoAwh = (get_eeenum(&inp, "awh", yesno_names, wi) != 0);
2073 ir->awhParams = gmx::readAndCheckAwhParams(&inp, ir, wi);
2077 gmx_fatal(FARGS, "AWH biasing is only compatible with COM pulling turned on");
2081 /* Enforced rotation */
2082 printStringNewline(&inp, "ENFORCED ROTATION");
2083 printStringNoNewline(&inp, "Enforced rotation: No or Yes");
2084 ir->bRot = (get_eeenum(&inp, "rotation", yesno_names, wi) != 0);
2088 is->rot_grp = read_rotparams(&inp, ir->rot, wi);
2091 /* Interactive MD */
2093 printStringNewline(&inp, "Group to display and/or manipulate in interactive MD session");
2094 setStringEntry(&inp, "IMD-group", is->imd_grp, nullptr);
2095 if (is->imd_grp[0] != '\0')
2102 printStringNewline(&inp, "NMR refinement stuff");
2103 printStringNoNewline(&inp, "Distance restraints type: No, Simple or Ensemble");
2104 ir->eDisre = get_eeenum(&inp, "disre", edisre_names, wi);
2105 printStringNoNewline(&inp, "Force weighting of pairs in one distance restraint: Conservative or Equal");
2106 ir->eDisreWeighting = get_eeenum(&inp, "disre-weighting", edisreweighting_names, wi);
2107 printStringNoNewline(&inp, "Use sqrt of the time averaged times the instantaneous violation");
2108 ir->bDisreMixed = (get_eeenum(&inp, "disre-mixed", yesno_names, wi) != 0);
2109 ir->dr_fc = get_ereal(&inp, "disre-fc", 1000.0, wi);
2110 ir->dr_tau = get_ereal(&inp, "disre-tau", 0.0, wi);
2111 printStringNoNewline(&inp, "Output frequency for pair distances to energy file");
2112 ir->nstdisreout = get_eint(&inp, "nstdisreout", 100, wi);
2113 printStringNoNewline(&inp, "Orientation restraints: No or Yes");
2114 opts->bOrire = (get_eeenum(&inp, "orire", yesno_names, wi) != 0);
2115 printStringNoNewline(&inp, "Orientation restraints force constant and tau for time averaging");
2116 ir->orires_fc = get_ereal(&inp, "orire-fc", 0.0, wi);
2117 ir->orires_tau = get_ereal(&inp, "orire-tau", 0.0, wi);
2118 setStringEntry(&inp, "orire-fitgrp", is->orirefitgrp, nullptr);
2119 printStringNoNewline(&inp, "Output frequency for trace(SD) and S to energy file");
2120 ir->nstorireout = get_eint(&inp, "nstorireout", 100, wi);
2122 /* free energy variables */
2123 printStringNewline(&inp, "Free energy variables");
2124 ir->efep = get_eeenum(&inp, "free-energy", efep_names, wi);
2125 setStringEntry(&inp, "couple-moltype", is->couple_moltype, nullptr);
2126 opts->couple_lam0 = get_eeenum(&inp, "couple-lambda0", couple_lam, wi);
2127 opts->couple_lam1 = get_eeenum(&inp, "couple-lambda1", couple_lam, wi);
2128 opts->bCoupleIntra = (get_eeenum(&inp, "couple-intramol", yesno_names, wi) != 0);
2130 fep->init_lambda = get_ereal(&inp, "init-lambda", -1, wi); /* start with -1 so
2132 it was not entered */
2133 fep->init_fep_state = get_eint(&inp, "init-lambda-state", -1, wi);
2134 fep->delta_lambda = get_ereal(&inp, "delta-lambda", 0.0, wi);
2135 fep->nstdhdl = get_eint(&inp, "nstdhdl", 50, wi);
2136 setStringEntry(&inp, "fep-lambdas", is->fep_lambda[efptFEP], nullptr);
2137 setStringEntry(&inp, "mass-lambdas", is->fep_lambda[efptMASS], nullptr);
2138 setStringEntry(&inp, "coul-lambdas", is->fep_lambda[efptCOUL], nullptr);
2139 setStringEntry(&inp, "vdw-lambdas", is->fep_lambda[efptVDW], nullptr);
2140 setStringEntry(&inp, "bonded-lambdas", is->fep_lambda[efptBONDED], nullptr);
2141 setStringEntry(&inp, "restraint-lambdas", is->fep_lambda[efptRESTRAINT], nullptr);
2142 setStringEntry(&inp, "temperature-lambdas", is->fep_lambda[efptTEMPERATURE], nullptr);
2143 fep->lambda_neighbors = get_eint(&inp, "calc-lambda-neighbors", 1, wi);
2144 setStringEntry(&inp, "init-lambda-weights", is->lambda_weights, nullptr);
2145 fep->edHdLPrintEnergy = get_eeenum(&inp, "dhdl-print-energy", edHdLPrintEnergy_names, wi);
2146 fep->sc_alpha = get_ereal(&inp, "sc-alpha", 0.0, wi);
2147 fep->sc_power = get_eint(&inp, "sc-power", 1, wi);
2148 fep->sc_r_power = get_ereal(&inp, "sc-r-power", 6.0, wi);
2149 fep->sc_sigma = get_ereal(&inp, "sc-sigma", 0.3, wi);
2150 fep->bScCoul = (get_eeenum(&inp, "sc-coul", yesno_names, wi) != 0);
2151 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2152 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2153 fep->separate_dhdl_file = get_eeenum(&inp, "separate-dhdl-file", separate_dhdl_file_names, wi);
2154 fep->dhdl_derivatives = get_eeenum(&inp, "dhdl-derivatives", dhdl_derivatives_names, wi);
2155 fep->dh_hist_size = get_eint(&inp, "dh_hist_size", 0, wi);
2156 fep->dh_hist_spacing = get_ereal(&inp, "dh_hist_spacing", 0.1, wi);
2158 /* Non-equilibrium MD stuff */
2159 printStringNewline(&inp, "Non-equilibrium MD stuff");
2160 setStringEntry(&inp, "acc-grps", is->accgrps, nullptr);
2161 setStringEntry(&inp, "accelerate", is->acc, nullptr);
2162 setStringEntry(&inp, "freezegrps", is->freeze, nullptr);
2163 setStringEntry(&inp, "freezedim", is->frdim, nullptr);
2164 ir->cos_accel = get_ereal(&inp, "cos-acceleration", 0, wi);
2165 setStringEntry(&inp, "deform", is->deform, nullptr);
2167 /* simulated tempering variables */
2168 printStringNewline(&inp, "simulated tempering variables");
2169 ir->bSimTemp = (get_eeenum(&inp, "simulated-tempering", yesno_names, wi) != 0);
2170 ir->simtempvals->eSimTempScale = get_eeenum(&inp, "simulated-tempering-scaling", esimtemp_names, wi);
2171 ir->simtempvals->simtemp_low = get_ereal(&inp, "sim-temp-low", 300.0, wi);
2172 ir->simtempvals->simtemp_high = get_ereal(&inp, "sim-temp-high", 300.0, wi);
2174 /* expanded ensemble variables */
2175 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2177 read_expandedparams(&inp, expand, wi);
2180 /* Electric fields */
2182 gmx::KeyValueTreeObject convertedValues = flatKeyValueTreeFromInpFile(inp);
2183 gmx::KeyValueTreeTransformer transform;
2184 transform.rules()->addRule()
2185 .keyMatchType("/", gmx::StringCompareType::CaseAndDashInsensitive);
2186 mdModules->initMdpTransform(transform.rules());
2187 for (const auto &path : transform.mappedPaths())
2189 GMX_ASSERT(path.size() == 1, "Inconsistent mapping back to mdp options");
2190 mark_einp_set(inp, path[0].c_str());
2192 MdpErrorHandler errorHandler(wi);
2194 = transform.transform(convertedValues, &errorHandler);
2195 ir->params = new gmx::KeyValueTreeObject(result.object());
2196 mdModules->adjustInputrecBasedOnModules(ir);
2197 errorHandler.setBackMapping(result.backMapping());
2198 mdModules->assignOptionsToModules(*ir->params, &errorHandler);
2201 /* Ion/water position swapping ("computational electrophysiology") */
2202 printStringNewline(&inp, "Ion/water position swapping for computational electrophysiology setups");
2203 printStringNoNewline(&inp, "Swap positions along direction: no, X, Y, Z");
2204 ir->eSwapCoords = get_eeenum(&inp, "swapcoords", eSwapTypes_names, wi);
2205 if (ir->eSwapCoords != eswapNO)
2212 printStringNoNewline(&inp, "Swap attempt frequency");
2213 ir->swap->nstswap = get_eint(&inp, "swap-frequency", 1, wi);
2214 printStringNoNewline(&inp, "Number of ion types to be controlled");
2215 nIonTypes = get_eint(&inp, "iontypes", 1, wi);
2218 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2220 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2221 snew(ir->swap->grp, ir->swap->ngrp);
2222 for (i = 0; i < ir->swap->ngrp; i++)
2224 snew(ir->swap->grp[i].molname, STRLEN);
2226 printStringNoNewline(&inp, "Two index groups that contain the compartment-partitioning atoms");
2227 setStringEntry(&inp, "split-group0", ir->swap->grp[eGrpSplit0].molname, nullptr);
2228 setStringEntry(&inp, "split-group1", ir->swap->grp[eGrpSplit1].molname, nullptr);
2229 printStringNoNewline(&inp, "Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2230 ir->swap->massw_split[0] = (get_eeenum(&inp, "massw-split0", yesno_names, wi) != 0);
2231 ir->swap->massw_split[1] = (get_eeenum(&inp, "massw-split1", yesno_names, wi) != 0);
2233 printStringNoNewline(&inp, "Name of solvent molecules");
2234 setStringEntry(&inp, "solvent-group", ir->swap->grp[eGrpSolvent].molname, nullptr);
2236 printStringNoNewline(&inp, "Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2237 printStringNoNewline(&inp, "Note that the split cylinder settings do not have an influence on the swapping protocol,");
2238 printStringNoNewline(&inp, "however, if correctly defined, the permeation events are recorded per channel");
2239 ir->swap->cyl0r = get_ereal(&inp, "cyl0-r", 2.0, wi);
2240 ir->swap->cyl0u = get_ereal(&inp, "cyl0-up", 1.0, wi);
2241 ir->swap->cyl0l = get_ereal(&inp, "cyl0-down", 1.0, wi);
2242 ir->swap->cyl1r = get_ereal(&inp, "cyl1-r", 2.0, wi);
2243 ir->swap->cyl1u = get_ereal(&inp, "cyl1-up", 1.0, wi);
2244 ir->swap->cyl1l = get_ereal(&inp, "cyl1-down", 1.0, wi);
2246 printStringNoNewline(&inp, "Average the number of ions per compartment over these many swap attempt steps");
2247 ir->swap->nAverage = get_eint(&inp, "coupl-steps", 10, wi);
2249 printStringNoNewline(&inp, "Names of the ion types that can be exchanged with solvent molecules,");
2250 printStringNoNewline(&inp, "and the requested number of ions of this type in compartments A and B");
2251 printStringNoNewline(&inp, "-1 means fix the numbers as found in step 0");
2252 for (i = 0; i < nIonTypes; i++)
2254 int ig = eSwapFixedGrpNR + i;
2256 sprintf(buf, "iontype%d-name", i);
2257 setStringEntry(&inp, buf, ir->swap->grp[ig].molname, nullptr);
2258 sprintf(buf, "iontype%d-in-A", i);
2259 ir->swap->grp[ig].nmolReq[0] = get_eint(&inp, buf, -1, wi);
2260 sprintf(buf, "iontype%d-in-B", i);
2261 ir->swap->grp[ig].nmolReq[1] = get_eint(&inp, buf, -1, wi);
2264 printStringNoNewline(&inp, "By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2265 printStringNoNewline(&inp, "at maximum distance (= bulk concentration) to the split group layers. However,");
2266 printStringNoNewline(&inp, "an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2267 printStringNoNewline(&inp, "towards one of the compartment-partitioning layers (at +/- 1.0).");
2268 ir->swap->bulkOffset[0] = get_ereal(&inp, "bulk-offsetA", 0.0, wi);
2269 ir->swap->bulkOffset[1] = get_ereal(&inp, "bulk-offsetB", 0.0, wi);
2270 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2271 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2273 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2276 printStringNoNewline(&inp, "Start to swap ions if threshold difference to requested count is reached");
2277 ir->swap->threshold = get_ereal(&inp, "threshold", 1.0, wi);
2280 /* AdResS is no longer supported, but we need grompp to be able to
2281 refuse to process old .mdp files that used it. */
2282 ir->bAdress = (get_eeenum(&inp, "adress", no_names, wi) != 0);
2284 /* User defined thingies */
2285 printStringNewline(&inp, "User defined thingies");
2286 setStringEntry(&inp, "user1-grps", is->user1, nullptr);
2287 setStringEntry(&inp, "user2-grps", is->user2, nullptr);
2288 ir->userint1 = get_eint(&inp, "userint1", 0, wi);
2289 ir->userint2 = get_eint(&inp, "userint2", 0, wi);
2290 ir->userint3 = get_eint(&inp, "userint3", 0, wi);
2291 ir->userint4 = get_eint(&inp, "userint4", 0, wi);
2292 ir->userreal1 = get_ereal(&inp, "userreal1", 0, wi);
2293 ir->userreal2 = get_ereal(&inp, "userreal2", 0, wi);
2294 ir->userreal3 = get_ereal(&inp, "userreal3", 0, wi);
2295 ir->userreal4 = get_ereal(&inp, "userreal4", 0, wi);
2299 gmx::TextOutputFile stream(mdparout);
2300 write_inpfile(&stream, mdparout, &inp, FALSE, writeMdpHeader, wi);
2302 // Transform module data into a flat key-value tree for output.
2303 gmx::KeyValueTreeBuilder builder;
2304 gmx::KeyValueTreeObjectBuilder builderObject = builder.rootObject();
2305 mdModules->buildMdpOutput(&builderObject);
2307 gmx::TextWriter writer(&stream);
2308 writeKeyValueTreeAsMdp(&writer, builder.build());
2313 /* Process options if necessary */
2314 for (m = 0; m < 2; m++)
2316 for (i = 0; i < 2*DIM; i++)
2325 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2327 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 1)");
2329 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2331 case epctSEMIISOTROPIC:
2332 case epctSURFACETENSION:
2333 if (sscanf(dumstr[m], "%lf%lf", &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2335 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 2)");
2337 dumdub[m][YY] = dumdub[m][XX];
2339 case epctANISOTROPIC:
2340 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2341 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2342 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2344 warning_error(wi, "Pressure coupling incorrect number of values (I need exactly 6)");
2348 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2349 epcoupltype_names[ir->epct]);
2353 clear_mat(ir->ref_p);
2354 clear_mat(ir->compress);
2355 for (i = 0; i < DIM; i++)
2357 ir->ref_p[i][i] = dumdub[1][i];
2358 ir->compress[i][i] = dumdub[0][i];
2360 if (ir->epct == epctANISOTROPIC)
2362 ir->ref_p[XX][YY] = dumdub[1][3];
2363 ir->ref_p[XX][ZZ] = dumdub[1][4];
2364 ir->ref_p[YY][ZZ] = dumdub[1][5];
2365 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2367 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2369 ir->compress[XX][YY] = dumdub[0][3];
2370 ir->compress[XX][ZZ] = dumdub[0][4];
2371 ir->compress[YY][ZZ] = dumdub[0][5];
2372 for (i = 0; i < DIM; i++)
2374 for (m = 0; m < i; m++)
2376 ir->ref_p[i][m] = ir->ref_p[m][i];
2377 ir->compress[i][m] = ir->compress[m][i];
2382 if (ir->comm_mode == ecmNO)
2387 opts->couple_moltype = nullptr;
2388 if (strlen(is->couple_moltype) > 0)
2390 if (ir->efep != efepNO)
2392 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2393 if (opts->couple_lam0 == opts->couple_lam1)
2395 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2397 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2398 opts->couple_lam1 == ecouplamNONE))
2400 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2405 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2408 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2409 if (ir->efep != efepNO)
2411 if (fep->delta_lambda > 0)
2413 ir->efep = efepSLOWGROWTH;
2417 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2419 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2420 warning_note(wi, "Old option for dhdl-print-energy given: "
2421 "changing \"yes\" to \"total\"\n");
2424 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2426 /* always print out the energy to dhdl if we are doing
2427 expanded ensemble, since we need the total energy for
2428 analysis if the temperature is changing. In some
2429 conditions one may only want the potential energy, so
2430 we will allow that if the appropriate mdp setting has
2431 been enabled. Otherwise, total it is:
2433 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2436 if ((ir->efep != efepNO) || ir->bSimTemp)
2438 ir->bExpanded = FALSE;
2439 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2441 ir->bExpanded = TRUE;
2443 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2444 if (ir->bSimTemp) /* done after fep params */
2446 do_simtemp_params(ir);
2449 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2450 * setup and not on the old way of specifying the free-energy setup,
2451 * we should check for using soft-core when not needed, since that
2452 * can complicate the sampling significantly.
2453 * Note that we only check for the automated coupling setup.
2454 * If the (advanced) user does FEP through manual topology changes,
2455 * this check will not be triggered.
2457 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2458 ir->fepvals->sc_alpha != 0 &&
2459 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2460 couple_lambda_has_vdw_on(opts->couple_lam1)))
2462 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.");
2467 ir->fepvals->n_lambda = 0;
2470 /* WALL PARAMETERS */
2472 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2474 /* ORIENTATION RESTRAINT PARAMETERS */
2476 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, nullptr) != 1)
2478 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2481 /* DEFORMATION PARAMETERS */
2483 clear_mat(ir->deform);
2484 for (i = 0; i < 6; i++)
2489 double gmx_unused canary;
2490 int ndeform = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf %lf",
2491 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2492 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]), &canary);
2494 if (strlen(is->deform) > 0 && ndeform != 6)
2496 warning_error(wi, gmx::formatString("Cannot parse exactly 6 box deformation velocities from string '%s'", is->deform).c_str());
2498 for (i = 0; i < 3; i++)
2500 ir->deform[i][i] = dumdub[0][i];
2502 ir->deform[YY][XX] = dumdub[0][3];
2503 ir->deform[ZZ][XX] = dumdub[0][4];
2504 ir->deform[ZZ][YY] = dumdub[0][5];
2505 if (ir->epc != epcNO)
2507 for (i = 0; i < 3; i++)
2509 for (j = 0; j <= i; j++)
2511 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2513 warning_error(wi, "A box element has deform set and compressibility > 0");
2517 for (i = 0; i < 3; i++)
2519 for (j = 0; j < i; j++)
2521 if (ir->deform[i][j] != 0)
2523 for (m = j; m < DIM; m++)
2525 if (ir->compress[m][j] != 0)
2527 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.");
2528 warning(wi, warn_buf);
2536 /* Ion/water position swapping checks */
2537 if (ir->eSwapCoords != eswapNO)
2539 if (ir->swap->nstswap < 1)
2541 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2543 if (ir->swap->nAverage < 1)
2545 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2547 if (ir->swap->threshold < 1.0)
2549 warning_error(wi, "Ion count threshold must be at least 1.\n");
2557 static int search_QMstring(const char *s, int ng, const char *gn[])
2559 /* same as normal search_string, but this one searches QM strings */
2562 for (i = 0; (i < ng); i++)
2564 if (gmx_strcasecmp(s, gn[i]) == 0)
2570 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2571 } /* search_QMstring */
2573 /* We would like gn to be const as well, but C doesn't allow this */
2574 /* TODO this is utility functionality (search for the index of a
2575 string in a collection), so should be refactored and located more
2577 int search_string(const char *s, int ng, char *gn[])
2581 for (i = 0; (i < ng); i++)
2583 if (gmx_strcasecmp(s, gn[i]) == 0)
2590 "Group %s referenced in the .mdp file was not found in the index file.\n"
2591 "Group names must match either [moleculetype] names or custom index group\n"
2592 "names, in which case you must supply an index file to the '-n' option\n"
2597 static bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2598 t_blocka *block, char *gnames[],
2599 int gtype, int restnm,
2600 int grptp, bool bVerbose,
2603 unsigned short *cbuf;
2604 t_grps *grps = &(groups->grps[gtype]);
2605 int i, j, gid, aj, ognr, ntot = 0;
2608 char warn_buf[STRLEN];
2610 title = gtypes[gtype];
2613 /* Mark all id's as not set */
2614 for (i = 0; (i < natoms); i++)
2619 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2620 for (i = 0; (i < ng); i++)
2622 /* Lookup the group name in the block structure */
2623 gid = search_string(ptrs[i], block->nr, gnames);
2624 if ((grptp != egrptpONE) || (i == 0))
2626 grps->nm_ind[grps->nr++] = gid;
2629 /* Now go over the atoms in the group */
2630 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2635 /* Range checking */
2636 if ((aj < 0) || (aj >= natoms))
2638 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2640 /* Lookup up the old group number */
2644 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2645 aj+1, title, ognr+1, i+1);
2649 /* Store the group number in buffer */
2650 if (grptp == egrptpONE)
2663 /* Now check whether we have done all atoms */
2667 if (grptp == egrptpALL)
2669 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2670 natoms-ntot, title);
2672 else if (grptp == egrptpPART)
2674 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2675 natoms-ntot, title);
2676 warning_note(wi, warn_buf);
2678 /* Assign all atoms currently unassigned to a rest group */
2679 for (j = 0; (j < natoms); j++)
2681 if (cbuf[j] == NOGID)
2687 if (grptp != egrptpPART)
2692 "Making dummy/rest group for %s containing %d elements\n",
2693 title, natoms-ntot);
2695 /* Add group name "rest" */
2696 grps->nm_ind[grps->nr] = restnm;
2698 /* Assign the rest name to all atoms not currently assigned to a group */
2699 for (j = 0; (j < natoms); j++)
2701 if (cbuf[j] == NOGID)
2710 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2712 /* All atoms are part of one (or no) group, no index required */
2713 groups->ngrpnr[gtype] = 0;
2714 groups->grpnr[gtype] = nullptr;
2718 groups->ngrpnr[gtype] = natoms;
2719 snew(groups->grpnr[gtype], natoms);
2720 for (j = 0; (j < natoms); j++)
2722 groups->grpnr[gtype][j] = cbuf[j];
2728 return (bRest && grptp == egrptpPART);
2731 static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2734 const gmx_groups_t *groups;
2735 pull_params_t *pull;
2736 int natoms, ai, aj, i, j, d, g, imin, jmin;
2738 int *nrdf2, *na_vcm, na_tot;
2739 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2741 gmx_mtop_atomloop_all_t aloop;
2745 * First calc 3xnr-atoms for each group
2746 * then subtract half a degree of freedom for each constraint
2748 * Only atoms and nuclei contribute to the degrees of freedom...
2753 groups = &mtop->groups;
2754 natoms = mtop->natoms;
2756 /* Allocate one more for a possible rest group */
2757 /* We need to sum degrees of freedom into doubles,
2758 * since floats give too low nrdf's above 3 million atoms.
2760 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2761 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2762 snew(dof_vcm, groups->grps[egcVCM].nr+1);
2763 snew(na_vcm, groups->grps[egcVCM].nr+1);
2764 snew(nrdf_vcm_sub, groups->grps[egcVCM].nr+1);
2766 for (i = 0; i < groups->grps[egcTC].nr; i++)
2770 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2773 clear_ivec(dof_vcm[i]);
2775 nrdf_vcm_sub[i] = 0;
2778 snew(nrdf2, natoms);
2779 aloop = gmx_mtop_atomloop_all_init(mtop);
2781 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2784 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2786 g = ggrpnr(groups, egcFREEZE, i);
2787 for (d = 0; d < DIM; d++)
2789 if (opts->nFreeze[g][d] == 0)
2791 /* Add one DOF for particle i (counted as 2*1) */
2793 /* VCM group i has dim d as a DOF */
2794 dof_vcm[ggrpnr(groups, egcVCM, i)][d] = 1;
2797 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2798 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2803 for (const gmx_molblock_t &molb : mtop->molblock)
2805 const gmx_moltype_t &molt = mtop->moltype[molb.type];
2806 atom = molt.atoms.atom;
2807 for (mol = 0; mol < molb.nmol; mol++)
2809 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2811 ia = molt.ilist[ftype].iatoms;
2812 for (i = 0; i < molt.ilist[ftype].nr; )
2814 /* Subtract degrees of freedom for the constraints,
2815 * if the particles still have degrees of freedom left.
2816 * If one of the particles is a vsite or a shell, then all
2817 * constraint motion will go there, but since they do not
2818 * contribute to the constraints the degrees of freedom do not
2823 if (((atom[ia[1]].ptype == eptNucleus) ||
2824 (atom[ia[1]].ptype == eptAtom)) &&
2825 ((atom[ia[2]].ptype == eptNucleus) ||
2826 (atom[ia[2]].ptype == eptAtom)))
2844 imin = std::min(imin, nrdf2[ai]);
2845 jmin = std::min(jmin, nrdf2[aj]);
2848 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2849 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2850 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2851 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2853 ia += interaction_function[ftype].nratoms+1;
2854 i += interaction_function[ftype].nratoms+1;
2857 ia = molt.ilist[F_SETTLE].iatoms;
2858 for (i = 0; i < molt.ilist[F_SETTLE].nr; )
2860 /* Subtract 1 dof from every atom in the SETTLE */
2861 for (j = 0; j < 3; j++)
2864 imin = std::min(2, nrdf2[ai]);
2866 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2867 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2872 as += molt.atoms.nr;
2878 /* Correct nrdf for the COM constraints.
2879 * We correct using the TC and VCM group of the first atom
2880 * in the reference and pull group. If atoms in one pull group
2881 * belong to different TC or VCM groups it is anyhow difficult
2882 * to determine the optimal nrdf assignment.
2886 for (i = 0; i < pull->ncoord; i++)
2888 if (pull->coord[i].eType != epullCONSTRAINT)
2895 for (j = 0; j < 2; j++)
2897 const t_pull_group *pgrp;
2899 pgrp = &pull->group[pull->coord[i].group[j]];
2903 /* Subtract 1/2 dof from each group */
2905 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2906 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2907 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2909 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)]]);
2914 /* We need to subtract the whole DOF from group j=1 */
2921 if (ir->nstcomm != 0)
2925 /* We remove COM motion up to dim ndof_com() */
2926 ndim_rm_vcm = ndof_com(ir);
2928 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
2929 * the number of degrees of freedom in each vcm group when COM
2930 * translation is removed and 6 when rotation is removed as well.
2932 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2934 switch (ir->comm_mode)
2937 case ecmLINEAR_ACCELERATION_CORRECTION:
2938 nrdf_vcm_sub[j] = 0;
2939 for (d = 0; d < ndim_rm_vcm; d++)
2948 nrdf_vcm_sub[j] = 6;
2951 gmx_incons("Checking comm_mode");
2955 for (i = 0; i < groups->grps[egcTC].nr; i++)
2957 /* Count the number of atoms of TC group i for every VCM group */
2958 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2963 for (ai = 0; ai < natoms; ai++)
2965 if (ggrpnr(groups, egcTC, ai) == i)
2967 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2971 /* Correct for VCM removal according to the fraction of each VCM
2972 * group present in this TC group.
2974 nrdf_uc = nrdf_tc[i];
2976 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2978 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
2980 nrdf_tc[i] += nrdf_uc*(static_cast<double>(na_vcm[j])/static_cast<double>(na_tot))*
2981 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
2986 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2988 opts->nrdf[i] = nrdf_tc[i];
2989 if (opts->nrdf[i] < 0)
2994 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2995 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
3003 sfree(nrdf_vcm_sub);
3006 static bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
3007 const char *option, const char *val, int flag)
3009 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3010 * But since this is much larger than STRLEN, such a line can not be parsed.
3011 * The real maximum is the number of names that fit in a string: STRLEN/2.
3013 #define EGP_MAX (STRLEN/2)
3014 int nelem, i, j, k, nr;
3015 char *names[EGP_MAX];
3019 gnames = groups->grpname;
3021 nelem = str_nelem(val, EGP_MAX, names);
3024 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3026 nr = groups->grps[egcENER].nr;
3028 for (i = 0; i < nelem/2; i++)
3032 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
3038 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3039 names[2*i], option);
3043 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3049 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3050 names[2*i+1], option);
3052 if ((j < nr) && (k < nr))
3054 ir->opts.egp_flags[nr*j+k] |= flag;
3055 ir->opts.egp_flags[nr*k+j] |= flag;
3064 static void make_swap_groups(
3069 int ig = -1, i = 0, gind;
3073 /* Just a quick check here, more thorough checks are in mdrun */
3074 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3076 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3079 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3080 for (ig = 0; ig < swap->ngrp; ig++)
3082 swapg = &swap->grp[ig];
3083 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3084 swapg->nat = grps->index[gind+1] - grps->index[gind];
3088 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3089 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3090 swap->grp[ig].molname, swapg->nat);
3091 snew(swapg->ind, swapg->nat);
3092 for (i = 0; i < swapg->nat; i++)
3094 swapg->ind[i] = grps->a[grps->index[gind]+i];
3099 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3105 static void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3110 ig = search_string(IMDgname, grps->nr, gnames);
3111 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3113 if (IMDgroup->nat > 0)
3115 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3116 IMDgname, IMDgroup->nat);
3117 snew(IMDgroup->ind, IMDgroup->nat);
3118 for (i = 0; i < IMDgroup->nat; i++)
3120 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3126 void do_index(const char* mdparin, const char *ndx,
3133 gmx_groups_t *groups;
3137 char warnbuf[STRLEN], **gnames;
3138 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3141 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3142 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3143 int i, j, k, restnm;
3144 bool bExcl, bTable, bAnneal, bRest;
3145 int nQMmethod, nQMbasis, nQMg;
3146 char warn_buf[STRLEN];
3151 fprintf(stderr, "processing index file...\n");
3156 snew(grps->index, 1);
3158 atoms_all = gmx_mtop_global_atoms(mtop);
3159 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3160 done_atom(&atoms_all);
3164 grps = init_index(ndx, &gnames);
3167 groups = &mtop->groups;
3168 natoms = mtop->natoms;
3169 symtab = &mtop->symtab;
3171 snew(groups->grpname, grps->nr+1);
3173 for (i = 0; (i < grps->nr); i++)
3175 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3177 groups->grpname[i] = put_symtab(symtab, "rest");
3179 srenew(gnames, grps->nr+1);
3180 gnames[restnm] = *(groups->grpname[i]);
3181 groups->ngrpname = grps->nr+1;
3183 set_warning_line(wi, mdparin, -1);
3185 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3186 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3187 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3188 if ((ntau_t != ntcg) || (nref_t != ntcg))
3190 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3191 "%d tau-t values", ntcg, nref_t, ntau_t);
3194 const bool useReferenceTemperature = integratorHasReferenceTemperature(ir);
3195 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3196 restnm, useReferenceTemperature ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3197 nr = groups->grps[egcTC].nr;
3199 snew(ir->opts.nrdf, nr);
3200 snew(ir->opts.tau_t, nr);
3201 snew(ir->opts.ref_t, nr);
3202 if (ir->eI == eiBD && ir->bd_fric == 0)
3204 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3207 if (useReferenceTemperature)
3211 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3215 for (i = 0; (i < nr); i++)
3217 ir->opts.tau_t[i] = strtod(ptr1[i], &endptr);
3220 warning_error(wi, "Invalid value for mdp option tau-t. tau-t should only consist of real numbers separated by spaces.");
3222 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3224 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3225 warning_error(wi, warn_buf);
3228 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3230 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.");
3233 if (ir->opts.tau_t[i] >= 0)
3235 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3238 if (ir->etc != etcNO && ir->nsttcouple == -1)
3240 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3245 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3247 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");
3249 if (ir->epc == epcMTTK)
3251 if (ir->etc != etcNOSEHOOVER)
3253 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3257 if (ir->nstpcouple != ir->nsttcouple)
3259 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3260 ir->nstpcouple = ir->nsttcouple = mincouple;
3261 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);
3262 warning_note(wi, warn_buf);
3267 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3268 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3270 if (ETC_ANDERSEN(ir->etc))
3272 if (ir->nsttcouple != 1)
3275 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");
3276 warning_note(wi, warn_buf);
3279 nstcmin = tcouple_min_integration_steps(ir->etc);
3282 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3284 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3285 ETCOUPLTYPE(ir->etc),
3287 ir->nsttcouple*ir->delta_t);
3288 warning(wi, warn_buf);
3291 for (i = 0; (i < nr); i++)
3293 ir->opts.ref_t[i] = strtod(ptr2[i], &endptr);
3296 warning_error(wi, "Invalid value for mdp option ref-t. ref-t should only consist of real numbers separated by spaces.");
3298 if (ir->opts.ref_t[i] < 0)
3300 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3303 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3304 if we are in this conditional) if mc_temp is negative */
3305 if (ir->expandedvals->mc_temp < 0)
3307 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3311 /* Simulated annealing for each group. There are nr groups */
3312 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3313 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3317 if (nSA > 0 && nSA != nr)
3319 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3323 snew(ir->opts.annealing, nr);
3324 snew(ir->opts.anneal_npoints, nr);
3325 snew(ir->opts.anneal_time, nr);
3326 snew(ir->opts.anneal_temp, nr);
3327 for (i = 0; i < nr; i++)
3329 ir->opts.annealing[i] = eannNO;
3330 ir->opts.anneal_npoints[i] = 0;
3331 ir->opts.anneal_time[i] = nullptr;
3332 ir->opts.anneal_temp[i] = nullptr;
3337 for (i = 0; i < nr; i++)
3339 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3341 ir->opts.annealing[i] = eannNO;
3343 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3345 ir->opts.annealing[i] = eannSINGLE;
3348 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3350 ir->opts.annealing[i] = eannPERIODIC;
3356 /* Read the other fields too */
3357 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3358 if (nSA_points != nSA)
3360 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3362 for (k = 0, i = 0; i < nr; i++)
3364 ir->opts.anneal_npoints[i] = strtol(ptr1[i], &endptr, 10);
3367 warning_error(wi, "Invalid value for mdp option annealing-npoints. annealing should only consist of integers separated by spaces.");
3369 if (ir->opts.anneal_npoints[i] == 1)
3371 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3373 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3374 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3375 k += ir->opts.anneal_npoints[i];
3378 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3381 gmx_fatal(FARGS, "Found %d annealing-time values, wanted %d\n", nSA_time, k);
3383 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3386 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3389 for (i = 0, k = 0; i < nr; i++)
3392 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3394 ir->opts.anneal_time[i][j] = strtod(ptr1[k], &endptr);
3397 warning_error(wi, "Invalid value for mdp option anneal-time. anneal-time should only consist of real numbers separated by spaces.");
3399 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], &endptr);
3402 warning_error(wi, "Invalid value for anneal-temp. anneal-temp should only consist of real numbers separated by spaces.");
3406 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3408 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3414 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3416 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3417 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3420 if (ir->opts.anneal_temp[i][j] < 0)
3422 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3427 /* Print out some summary information, to make sure we got it right */
3428 for (i = 0, k = 0; i < nr; i++)
3430 if (ir->opts.annealing[i] != eannNO)
3432 j = groups->grps[egcTC].nm_ind[i];
3433 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3434 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3435 ir->opts.anneal_npoints[i]);
3436 fprintf(stderr, "Time (ps) Temperature (K)\n");
3437 /* All terms except the last one */
3438 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3440 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3443 /* Finally the last one */
3444 j = ir->opts.anneal_npoints[i]-1;
3445 if (ir->opts.annealing[i] == eannSINGLE)
3447 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3451 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3452 if (std::fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3454 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3465 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3467 make_pull_coords(ir->pull);
3472 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3475 if (ir->eSwapCoords != eswapNO)
3477 make_swap_groups(ir->swap, grps, gnames);
3480 /* Make indices for IMD session */
3483 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3486 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3487 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3488 if (nacg*DIM != nacc)
3490 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3493 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3494 restnm, egrptpALL_GENREST, bVerbose, wi);
3495 nr = groups->grps[egcACC].nr;
3496 snew(ir->opts.acc, nr);
3497 ir->opts.ngacc = nr;
3499 for (i = k = 0; (i < nacg); i++)
3501 for (j = 0; (j < DIM); j++, k++)
3503 ir->opts.acc[i][j] = strtod(ptr1[k], &endptr);
3506 warning_error(wi, "Invalid value for mdp option accelerate. accelerate should only consist of real numbers separated by spaces.");
3510 for (; (i < nr); i++)
3512 for (j = 0; (j < DIM); j++)
3514 ir->opts.acc[i][j] = 0;
3518 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3519 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3520 if (nfrdim != DIM*nfreeze)
3522 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3525 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3526 restnm, egrptpALL_GENREST, bVerbose, wi);
3527 nr = groups->grps[egcFREEZE].nr;
3528 ir->opts.ngfrz = nr;
3529 snew(ir->opts.nFreeze, nr);
3530 for (i = k = 0; (i < nfreeze); i++)
3532 for (j = 0; (j < DIM); j++, k++)
3534 ir->opts.nFreeze[i][j] = static_cast<int>(gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3535 if (!ir->opts.nFreeze[i][j])
3537 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3539 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3540 "(not %s)", ptr1[k]);
3541 warning(wi, warn_buf);
3546 for (; (i < nr); i++)
3548 for (j = 0; (j < DIM); j++)
3550 ir->opts.nFreeze[i][j] = 0;
3554 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3555 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3556 restnm, egrptpALL_GENREST, bVerbose, wi);
3557 add_wall_energrps(groups, ir->nwall, symtab);
3558 ir->opts.ngener = groups->grps[egcENER].nr;
3559 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3561 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3562 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3565 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3566 "This may lead to artifacts.\n"
3567 "In most cases one should use one group for the whole system.");
3570 /* Now we have filled the freeze struct, so we can calculate NRDF */
3571 calc_nrdf(mtop, ir, gnames);
3573 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3574 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3575 restnm, egrptpALL_GENREST, bVerbose, wi);
3576 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3577 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3578 restnm, egrptpALL_GENREST, bVerbose, wi);
3579 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3580 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3581 restnm, egrptpONE, bVerbose, wi);
3582 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3583 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3584 restnm, egrptpALL_GENREST, bVerbose, wi);
3586 /* QMMM input processing */
3587 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3588 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3589 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3590 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3592 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3593 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3595 /* group rest, if any, is always MM! */
3596 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3597 restnm, egrptpALL_GENREST, bVerbose, wi);
3598 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3599 ir->opts.ngQM = nQMg;
3600 snew(ir->opts.QMmethod, nr);
3601 snew(ir->opts.QMbasis, nr);
3602 for (i = 0; i < nr; i++)
3604 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3605 * converted to the corresponding enum in names.c
3607 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3609 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3613 str_nelem(is->QMmult, MAXPTR, ptr1);
3614 str_nelem(is->QMcharge, MAXPTR, ptr2);
3615 str_nelem(is->bSH, MAXPTR, ptr3);
3616 snew(ir->opts.QMmult, nr);
3617 snew(ir->opts.QMcharge, nr);
3618 snew(ir->opts.bSH, nr);
3620 for (i = 0; i < nr; i++)
3622 ir->opts.QMmult[i] = strtol(ptr1[i], &endptr, 10);
3625 warning_error(wi, "Invalid value for mdp option QMmult. QMmult should only consist of integers separated by spaces.");
3627 ir->opts.QMcharge[i] = strtol(ptr2[i], &endptr, 10);
3630 warning_error(wi, "Invalid value for mdp option QMcharge. QMcharge should only consist of integers separated by spaces.");
3632 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3635 str_nelem(is->CASelectrons, MAXPTR, ptr1);
3636 str_nelem(is->CASorbitals, MAXPTR, ptr2);
3637 snew(ir->opts.CASelectrons, nr);
3638 snew(ir->opts.CASorbitals, nr);
3639 for (i = 0; i < nr; i++)
3641 ir->opts.CASelectrons[i] = strtol(ptr1[i], &endptr, 10);
3644 warning_error(wi, "Invalid value for mdp option CASelectrons. CASelectrons should only consist of integers separated by spaces.");
3646 ir->opts.CASorbitals[i] = strtol(ptr2[i], &endptr, 10);
3649 warning_error(wi, "Invalid value for mdp option CASorbitals. CASorbitals should only consist of integers separated by spaces.");
3653 str_nelem(is->SAon, MAXPTR, ptr1);
3654 str_nelem(is->SAoff, MAXPTR, ptr2);
3655 str_nelem(is->SAsteps, MAXPTR, ptr3);
3656 snew(ir->opts.SAon, nr);
3657 snew(ir->opts.SAoff, nr);
3658 snew(ir->opts.SAsteps, nr);
3660 for (i = 0; i < nr; i++)
3662 ir->opts.SAon[i] = strtod(ptr1[i], &endptr);
3665 warning_error(wi, "Invalid value for mdp option SAon. SAon should only consist of real numbers separated by spaces.");
3667 ir->opts.SAoff[i] = strtod(ptr2[i], &endptr);
3670 warning_error(wi, "Invalid value for mdp option SAoff. SAoff should only consist of real numbers separated by spaces.");
3672 ir->opts.SAsteps[i] = strtol(ptr3[i], &endptr, 10);
3675 warning_error(wi, "Invalid value for mdp option SAsteps. SAsteps should only consist of integers separated by spaces.");
3678 /* end of QMMM input */
3682 for (i = 0; (i < egcNR); i++)
3684 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3685 for (j = 0; (j < groups->grps[i].nr); j++)
3687 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3689 fprintf(stderr, "\n");
3693 nr = groups->grps[egcENER].nr;
3694 snew(ir->opts.egp_flags, nr*nr);
3696 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3697 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3699 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3701 if (bExcl && EEL_FULL(ir->coulombtype))
3703 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3706 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3707 if (bTable && !(ir->vdwtype == evdwUSER) &&
3708 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3709 !(ir->coulombtype == eelPMEUSERSWITCH))
3711 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3714 for (i = 0; (i < grps->nr); i++)
3726 static void check_disre(gmx_mtop_t *mtop)
3728 gmx_ffparams_t *ffparams;
3729 t_functype *functype;
3731 int i, ndouble, ftype;
3732 int label, old_label;
3734 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3736 ffparams = &mtop->ffparams;
3737 functype = ffparams->functype;
3738 ip = ffparams->iparams;
3741 for (i = 0; i < ffparams->ntypes; i++)
3743 ftype = functype[i];
3744 if (ftype == F_DISRES)
3746 label = ip[i].disres.label;
3747 if (label == old_label)
3749 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3757 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3758 "probably the parameters for multiple pairs in one restraint "
3759 "are not identical\n", ndouble);
3764 static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3769 gmx_mtop_ilistloop_t iloop;
3770 const t_ilist *ilist;
3779 for (d = 0; d < DIM; d++)
3781 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3783 /* Check for freeze groups */
3784 for (g = 0; g < ir->opts.ngfrz; g++)
3786 for (d = 0; d < DIM; d++)
3788 if (ir->opts.nFreeze[g][d] != 0)
3796 /* Check for position restraints */
3797 iloop = gmx_mtop_ilistloop_init(sys);
3798 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3801 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3803 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3805 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3806 for (d = 0; d < DIM; d++)
3808 if (pr->posres.fcA[d] != 0)
3814 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3816 /* Check for flat-bottom posres */
3817 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3818 if (pr->fbposres.k != 0)
3820 switch (pr->fbposres.geom)
3822 case efbposresSPHERE:
3823 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3825 case efbposresCYLINDERX:
3826 AbsRef[YY] = AbsRef[ZZ] = 1;
3828 case efbposresCYLINDERY:
3829 AbsRef[XX] = AbsRef[ZZ] = 1;
3831 case efbposresCYLINDER:
3832 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3833 case efbposresCYLINDERZ:
3834 AbsRef[XX] = AbsRef[YY] = 1;
3836 case efbposresX: /* d=XX */
3837 case efbposresY: /* d=YY */
3838 case efbposresZ: /* d=ZZ */
3839 d = pr->fbposres.geom - efbposresX;
3843 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3844 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3852 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3856 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3857 bool *bC6ParametersWorkWithGeometricRules,
3858 bool *bC6ParametersWorkWithLBRules,
3859 bool *bLBRulesPossible)
3861 int ntypes, tpi, tpj;
3864 double c6i, c6j, c12i, c12j;
3865 double c6, c6_geometric, c6_LB;
3866 double sigmai, sigmaj, epsi, epsj;
3867 bool bCanDoLBRules, bCanDoGeometricRules;
3870 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3871 * force-field floating point parameters.
3874 ptr = getenv("GMX_LJCOMB_TOL");
3878 double gmx_unused canary;
3880 if (sscanf(ptr, "%lf%lf", &dbl, &canary) != 1)
3882 gmx_fatal(FARGS, "Could not parse a single floating-point number from GMX_LJCOMB_TOL (%s)", ptr);
3887 *bC6ParametersWorkWithLBRules = TRUE;
3888 *bC6ParametersWorkWithGeometricRules = TRUE;
3889 bCanDoLBRules = TRUE;
3890 ntypes = mtop->ffparams.atnr;
3891 snew(typecount, ntypes);
3892 gmx_mtop_count_atomtypes(mtop, state, typecount);
3893 *bLBRulesPossible = TRUE;
3894 for (tpi = 0; tpi < ntypes; ++tpi)
3896 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3897 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3898 for (tpj = tpi; tpj < ntypes; ++tpj)
3900 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3901 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3902 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3903 c6_geometric = std::sqrt(c6i * c6j);
3904 if (!gmx_numzero(c6_geometric))
3906 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3908 sigmai = gmx::sixthroot(c12i / c6i);
3909 sigmaj = gmx::sixthroot(c12j / c6j);
3910 epsi = c6i * c6i /(4.0 * c12i);
3911 epsj = c6j * c6j /(4.0 * c12j);
3912 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
3916 *bLBRulesPossible = FALSE;
3917 c6_LB = c6_geometric;
3919 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3924 *bC6ParametersWorkWithLBRules = FALSE;
3927 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3929 if (!bCanDoGeometricRules)
3931 *bC6ParametersWorkWithGeometricRules = FALSE;
3939 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3942 bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3944 check_combination_rule_differences(mtop, 0,
3945 &bC6ParametersWorkWithGeometricRules,
3946 &bC6ParametersWorkWithLBRules,
3948 if (ir->ljpme_combination_rule == eljpmeLB)
3950 if (!bC6ParametersWorkWithLBRules || !bLBRulesPossible)
3952 warning(wi, "You are using arithmetic-geometric combination rules "
3953 "in LJ-PME, but your non-bonded C6 parameters do not "
3954 "follow these rules.");
3959 if (!bC6ParametersWorkWithGeometricRules)
3961 if (ir->eDispCorr != edispcNO)
3963 warning_note(wi, "You are using geometric combination rules in "
3964 "LJ-PME, but your non-bonded C6 parameters do "
3965 "not follow these rules. "
3966 "This will introduce very small errors in the forces and energies in "
3967 "your simulations. Dispersion correction will correct total energy "
3968 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3972 warning_note(wi, "You are using geometric combination rules in "
3973 "LJ-PME, but your non-bonded C6 parameters do "
3974 "not follow these rules. "
3975 "This will introduce very small errors in the forces and energies in "
3976 "your simulations. If your system is homogeneous, consider using dispersion correction "
3977 "for the total energy and pressure.");
3983 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3986 char err_buf[STRLEN];
3991 gmx_mtop_atomloop_block_t aloopb;
3992 gmx_mtop_atomloop_all_t aloop;
3994 char warn_buf[STRLEN];
3996 set_warning_line(wi, mdparin, -1);
3998 if (ir->cutoff_scheme == ecutsVERLET &&
3999 ir->verletbuf_tol > 0 &&
4001 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
4002 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4004 /* Check if a too small Verlet buffer might potentially
4005 * cause more drift than the thermostat can couple off.
4007 /* Temperature error fraction for warning and suggestion */
4008 const real T_error_warn = 0.002;
4009 const real T_error_suggest = 0.001;
4010 /* For safety: 2 DOF per atom (typical with constraints) */
4011 const real nrdf_at = 2;
4012 real T, tau, max_T_error;
4017 for (i = 0; i < ir->opts.ngtc; i++)
4019 T = std::max(T, ir->opts.ref_t[i]);
4020 tau = std::max(tau, ir->opts.tau_t[i]);
4024 /* This is a worst case estimate of the temperature error,
4025 * assuming perfect buffer estimation and no cancelation
4026 * of errors. The factor 0.5 is because energy distributes
4027 * equally over Ekin and Epot.
4029 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4030 if (max_T_error > T_error_warn)
4032 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.",
4033 ir->verletbuf_tol, T, tau,
4035 100*T_error_suggest,
4036 ir->verletbuf_tol*T_error_suggest/max_T_error);
4037 warning(wi, warn_buf);
4042 if (ETC_ANDERSEN(ir->etc))
4046 for (i = 0; i < ir->opts.ngtc; i++)
4048 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4049 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4050 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
4051 i, ir->opts.tau_t[i]);
4052 CHECK(ir->opts.tau_t[i] < 0);
4055 if (ir->etc == etcANDERSENMASSIVE && ir->comm_mode != ecmNO)
4057 for (i = 0; i < ir->opts.ngtc; i++)
4059 int nsteps = static_cast<int>(ir->opts.tau_t[i]/ir->delta_t + 0.5);
4060 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);
4061 CHECK(nsteps % ir->nstcomm != 0);
4066 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4067 ir->comm_mode == ecmNO &&
4068 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4069 !ETC_ANDERSEN(ir->etc))
4071 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");
4074 /* Check for pressure coupling with absolute position restraints */
4075 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4077 absolute_reference(ir, sys, TRUE, AbsRef);
4079 for (m = 0; m < DIM; m++)
4081 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4083 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4091 aloopb = gmx_mtop_atomloop_block_init(sys);
4093 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4095 if (atom->q != 0 || atom->qB != 0)
4103 if (EEL_FULL(ir->coulombtype))
4106 "You are using full electrostatics treatment %s for a system without charges.\n"
4107 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4108 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4109 warning(wi, err_buf);
4114 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0)
4117 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4118 "You might want to consider using %s electrostatics.\n",
4120 warning_note(wi, err_buf);
4124 /* Check if combination rules used in LJ-PME are the same as in the force field */
4125 if (EVDW_PME(ir->vdwtype))
4127 check_combination_rules(ir, sys, wi);
4130 /* Generalized reaction field */
4131 if (ir->opts.ngtc == 0)
4133 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4135 CHECK(ir->coulombtype == eelGRF);
4139 sprintf(err_buf, "When using coulombtype = %s"
4140 " ref-t for temperature coupling should be > 0",
4142 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4146 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4148 for (m = 0; (m < DIM); m++)
4150 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4159 snew(mgrp, sys->groups.grps[egcACC].nr);
4160 aloop = gmx_mtop_atomloop_all_init(sys);
4162 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4164 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4167 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4169 for (m = 0; (m < DIM); m++)
4171 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4175 for (m = 0; (m < DIM); m++)
4177 if (fabs(acc[m]) > 1e-6)
4179 const char *dim[DIM] = { "X", "Y", "Z" };
4181 "Net Acceleration in %s direction, will %s be corrected\n",
4182 dim[m], ir->nstcomm != 0 ? "" : "not");
4183 if (ir->nstcomm != 0 && m < ndof_com(ir))
4186 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4188 ir->opts.acc[i][m] -= acc[m];
4196 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4197 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4199 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4207 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4209 if (ir->pull->coord[i].group[0] == 0 ||
4210 ir->pull->coord[i].group[1] == 0)
4212 absolute_reference(ir, sys, FALSE, AbsRef);
4213 for (m = 0; m < DIM; m++)
4215 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4217 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.");
4225 for (i = 0; i < 3; i++)
4227 for (m = 0; m <= i; m++)
4229 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4230 ir->deform[i][m] != 0)
4232 for (c = 0; c < ir->pull->ncoord; c++)
4234 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4235 ir->pull->coord[c].vec[m] != 0)
4237 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4248 void double_check(t_inputrec *ir, matrix box,
4249 bool bHasNormalConstraints,
4250 bool bHasAnyConstraints,
4254 char warn_buf[STRLEN];
4257 ptr = check_box(ir->ePBC, box);
4260 warning_error(wi, ptr);
4263 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4265 if (ir->shake_tol <= 0.0)
4267 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4269 warning_error(wi, warn_buf);
4273 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4275 /* If we have Lincs constraints: */
4276 if (ir->eI == eiMD && ir->etc == etcNO &&
4277 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4279 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4280 warning_note(wi, warn_buf);
4283 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4285 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4286 warning_note(wi, warn_buf);
4288 if (ir->epc == epcMTTK)
4290 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4294 if (bHasAnyConstraints && ir->epc == epcMTTK)
4296 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4299 if (ir->LincsWarnAngle > 90.0)
4301 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4302 warning(wi, warn_buf);
4303 ir->LincsWarnAngle = 90.0;
4306 if (ir->ePBC != epbcNONE)
4308 if (ir->nstlist == 0)
4310 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4312 if (ir->ns_type == ensGRID)
4314 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4316 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");
4317 warning_error(wi, warn_buf);
4322 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4323 if (2*ir->rlist >= min_size)
4325 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4326 warning_error(wi, warn_buf);
4329 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4336 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4340 real rvdw1, rvdw2, rcoul1, rcoul2;
4341 char warn_buf[STRLEN];
4343 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4347 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4352 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4358 if (rvdw1 + rvdw2 > ir->rlist ||
4359 rcoul1 + rcoul2 > ir->rlist)
4362 "The sum of the two largest charge group radii (%f) "
4363 "is larger than rlist (%f)\n",
4364 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4365 warning(wi, warn_buf);
4369 /* Here we do not use the zero at cut-off macro,
4370 * since user defined interactions might purposely
4371 * not be zero at the cut-off.
4373 if (ir_vdw_is_zero_at_cutoff(ir) &&
4374 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4376 sprintf(warn_buf, "The sum of the two largest charge group "
4377 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4378 "With exact cut-offs, better performance can be "
4379 "obtained with cutoff-scheme = %s, because it "
4380 "does not use charge groups at all.",
4382 ir->rlist, ir->rvdw,
4383 ecutscheme_names[ecutsVERLET]);
4386 warning(wi, warn_buf);
4390 warning_note(wi, warn_buf);
4393 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4394 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4396 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4397 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4399 ir->rlist, ir->rcoulomb,
4400 ecutscheme_names[ecutsVERLET]);
4403 warning(wi, warn_buf);
4407 warning_note(wi, warn_buf);