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, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
49 #include "gromacs/fileio/readinp.h"
50 #include "gromacs/fileio/warninp.h"
51 #include "gromacs/gmxlib/chargegroup.h"
52 #include "gromacs/gmxlib/network.h"
53 #include "gromacs/gmxpreprocess/toputil.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdlib/calc_verletbuf.h"
58 #include "gromacs/mdtypes/inputrec.h"
59 #include "gromacs/mdtypes/md_enums.h"
60 #include "gromacs/mdtypes/pull-params.h"
61 #include "gromacs/pbcutil/pbc.h"
62 #include "gromacs/topology/block.h"
63 #include "gromacs/topology/ifunc.h"
64 #include "gromacs/topology/index.h"
65 #include "gromacs/topology/mtop_util.h"
66 #include "gromacs/topology/symtab.h"
67 #include "gromacs/topology/topology.h"
68 #include "gromacs/utility/cstringutil.h"
69 #include "gromacs/utility/fatalerror.h"
70 #include "gromacs/utility/smalloc.h"
75 /* Resource parameters
76 * Do not change any of these until you read the instruction
77 * in readinp.h. Some cpp's do not take spaces after the backslash
78 * (like the c-shell), which will give you a very weird compiler
82 typedef struct t_inputrec_strings
84 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
85 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
86 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
87 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
88 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
90 char fep_lambda[efptNR][STRLEN];
91 char lambda_weights[STRLEN];
94 char anneal[STRLEN], anneal_npoints[STRLEN],
95 anneal_time[STRLEN], anneal_temp[STRLEN];
96 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
97 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
98 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
99 char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
100 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
102 } gmx_inputrec_strings;
104 static gmx_inputrec_strings *is = NULL;
106 void init_inputrec_strings()
110 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.");
115 void done_inputrec_strings()
123 egrptpALL, /* All particles have to be a member of a group. */
124 egrptpALL_GENREST, /* A rest group with name is generated for particles *
125 * that are not part of any group. */
126 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
127 * for the rest group. */
128 egrptpONE /* Merge all selected groups into one group, *
129 * make a rest group for the remaining particles. */
132 static const char *constraints[eshNR+1] = {
133 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
136 static const char *couple_lam[ecouplamNR+1] = {
137 "vdw-q", "vdw", "q", "none", NULL
140 void init_ir(t_inputrec *ir, t_gromppopts *opts)
142 snew(opts->include, STRLEN);
143 snew(opts->define, STRLEN);
144 snew(ir->fepvals, 1);
145 snew(ir->expandedvals, 1);
146 snew(ir->simtempvals, 1);
149 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
154 for (i = 0; i < ntemps; i++)
156 /* simple linear scaling -- allows more control */
157 if (simtemp->eSimTempScale == esimtempLINEAR)
159 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
161 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
163 simtemp->temperatures[i] = simtemp->simtemp_low * std::pow(simtemp->simtemp_high/simtemp->simtemp_low, static_cast<real>((1.0*i)/(ntemps-1)));
165 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
167 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*(std::expm1(temperature_lambdas[i])/std::expm1(1.0));
172 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
173 gmx_fatal(FARGS, errorstr);
180 static void _low_check(gmx_bool b, const char *s, warninp_t wi)
184 warning_error(wi, s);
188 static void check_nst(const char *desc_nst, int nst,
189 const char *desc_p, int *p,
194 if (*p > 0 && *p % nst != 0)
196 /* Round up to the next multiple of nst */
197 *p = ((*p)/nst + 1)*nst;
198 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
199 desc_p, desc_nst, desc_p, *p);
204 static gmx_bool ir_NVE(const t_inputrec *ir)
206 return (EI_MD(ir->eI) && ir->etc == etcNO);
209 static int lcd(int n1, int n2)
214 for (i = 2; (i <= n1 && i <= n2); i++)
216 if (n1 % i == 0 && n2 % i == 0)
225 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
227 if (*eintmod == eintmodPOTSHIFT_VERLET)
229 if (ir->cutoff_scheme == ecutsVERLET)
231 *eintmod = eintmodPOTSHIFT;
235 *eintmod = eintmodNONE;
240 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
242 /* Check internal consistency.
243 * NOTE: index groups are not set here yet, don't check things
244 * like temperature coupling group options here, but in triple_check
247 /* Strange macro: first one fills the err_buf, and then one can check
248 * the condition, which will print the message and increase the error
251 #define CHECK(b) _low_check(b, err_buf, wi)
252 char err_buf[256], warn_buf[STRLEN];
255 t_lambda *fep = ir->fepvals;
256 t_expanded *expand = ir->expandedvals;
258 set_warning_line(wi, mdparin, -1);
260 if (ir->coulombtype == eelRF_NEC_UNSUPPORTED)
262 sprintf(warn_buf, "%s electrostatics is no longer supported",
263 eel_names[eelRF_NEC_UNSUPPORTED]);
264 warning_error(wi, warn_buf);
267 /* BASIC CUT-OFF STUFF */
268 if (ir->rcoulomb < 0)
270 warning_error(wi, "rcoulomb should be >= 0");
274 warning_error(wi, "rvdw should be >= 0");
277 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
279 warning_error(wi, "rlist should be >= 0");
281 sprintf(err_buf, "nstlist can not be smaller than 0. (If you were trying to use the heuristic neighbour-list update scheme for efficient buffering for improved energy conservation, please use the Verlet cut-off scheme instead.)");
282 CHECK(ir->nstlist < 0);
284 process_interaction_modifier(ir, &ir->coulomb_modifier);
285 process_interaction_modifier(ir, &ir->vdw_modifier);
287 if (ir->cutoff_scheme == ecutsGROUP)
290 "The group cutoff scheme is deprecated since GROMACS 5.0 and will be removed in a future "
291 "release when all interaction forms are supported for the verlet scheme. The verlet "
292 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
294 if (ir->rlist > 0 && ir->rlist < ir->rcoulomb)
296 gmx_fatal(FARGS, "rcoulomb must not be greater than rlist (twin-range schemes are not supported)");
298 if (ir->rlist > 0 && ir->rlist < ir->rvdw)
300 gmx_fatal(FARGS, "rvdw must not be greater than rlist (twin-range schemes are not supported)");
303 if (ir->rlist == 0 && ir->ePBC != epbcNONE)
305 warning_error(wi, "Can not have an infinite cut-off with PBC");
309 if (ir->cutoff_scheme == ecutsVERLET)
313 /* Normal Verlet type neighbor-list, currently only limited feature support */
314 if (inputrec2nboundeddim(ir) < 3)
316 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
319 // We don't (yet) have general Verlet kernels for rcoulomb!=rvdw
320 if (ir->rcoulomb != ir->rvdw)
322 // Since we have PME coulomb + LJ cut-off kernels with rcoulomb>rvdw
323 // for PME load balancing, we can support this exception.
324 bool bUsesPmeTwinRangeKernel = (EEL_PME_EWALD(ir->coulombtype) &&
325 ir->vdwtype == evdwCUT &&
326 ir->rcoulomb > ir->rvdw);
327 if (!bUsesPmeTwinRangeKernel)
329 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported (except for rcoulomb>rvdw with PME electrostatics)");
333 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
335 if (ir->vdw_modifier == eintmodNONE ||
336 ir->vdw_modifier == eintmodPOTSHIFT)
338 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
340 sprintf(warn_buf, "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
341 warning_note(wi, warn_buf);
343 ir->vdwtype = evdwCUT;
347 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
348 warning_error(wi, warn_buf);
352 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
354 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
356 if (!(ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype) ||
357 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
359 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
361 if (!(ir->coulomb_modifier == eintmodNONE ||
362 ir->coulomb_modifier == eintmodPOTSHIFT))
364 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
365 warning_error(wi, warn_buf);
368 if (ir->implicit_solvent != eisNO)
370 warning_error(wi, "Implicit solvent is not (yet) supported with the with Verlet lists.");
373 if (ir->nstlist <= 0)
375 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
378 if (ir->nstlist < 10)
380 warning_note(wi, "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
383 rc_max = std::max(ir->rvdw, ir->rcoulomb);
385 if (ir->verletbuf_tol <= 0)
387 if (ir->verletbuf_tol == 0)
389 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
392 if (ir->rlist < rc_max)
394 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
397 if (ir->rlist == rc_max && ir->nstlist > 1)
399 warning_note(wi, "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
404 if (ir->rlist > rc_max)
406 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
409 if (ir->nstlist == 1)
411 /* No buffer required */
416 if (EI_DYNAMICS(ir->eI))
418 if (inputrec2nboundeddim(ir) < 3)
420 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
422 /* Set rlist temporarily so we can continue processing */
427 /* Set the buffer to 5% of the cut-off */
428 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
434 /* GENERAL INTEGRATOR STUFF */
437 if (ir->etc != etcNO)
439 if (EI_RANDOM(ir->eI))
441 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. %s handles temperature coupling implicitly. See the documentation for more information on which parameters affect temperature for %s.", etcoupl_names[ir->etc], ei_names[ir->eI], ei_names[ir->eI]);
445 sprintf(warn_buf, "Setting tcoupl from '%s' to 'no'. Temperature coupling does not apply to %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
447 warning_note(wi, warn_buf);
451 if (ir->eI == eiVVAK)
453 sprintf(warn_buf, "Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s", ei_names[eiVVAK], ei_names[eiMD], ei_names[eiVV]);
454 warning_note(wi, warn_buf);
456 if (!EI_DYNAMICS(ir->eI))
458 if (ir->epc != epcNO)
460 sprintf(warn_buf, "Setting pcoupl from '%s' to 'no'. Pressure coupling does not apply to %s.", epcoupl_names[ir->epc], ei_names[ir->eI]);
461 warning_note(wi, warn_buf);
465 if (EI_DYNAMICS(ir->eI))
467 if (ir->nstcalcenergy < 0)
469 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
470 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
472 /* nstcalcenergy larger than nstener does not make sense.
473 * We ideally want nstcalcenergy=nstener.
477 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
481 ir->nstcalcenergy = ir->nstenergy;
485 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
486 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
487 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
490 const char *nsten = "nstenergy";
491 const char *nstdh = "nstdhdl";
492 const char *min_name = nsten;
493 int min_nst = ir->nstenergy;
495 /* find the smallest of ( nstenergy, nstdhdl ) */
496 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
497 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
499 min_nst = ir->fepvals->nstdhdl;
502 /* If the user sets nstenergy small, we should respect that */
504 "Setting nstcalcenergy (%d) equal to %s (%d)",
505 ir->nstcalcenergy, min_name, min_nst);
506 warning_note(wi, warn_buf);
507 ir->nstcalcenergy = min_nst;
510 if (ir->epc != epcNO)
512 if (ir->nstpcouple < 0)
514 ir->nstpcouple = ir_optimal_nstpcouple(ir);
518 if (ir->nstcalcenergy > 0)
520 if (ir->efep != efepNO)
522 /* nstdhdl should be a multiple of nstcalcenergy */
523 check_nst("nstcalcenergy", ir->nstcalcenergy,
524 "nstdhdl", &ir->fepvals->nstdhdl, wi);
525 /* nstexpanded should be a multiple of nstcalcenergy */
526 check_nst("nstcalcenergy", ir->nstcalcenergy,
527 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
529 /* for storing exact averages nstenergy should be
530 * a multiple of nstcalcenergy
532 check_nst("nstcalcenergy", ir->nstcalcenergy,
533 "nstenergy", &ir->nstenergy, wi);
537 if (ir->nsteps == 0 && !ir->bContinuation)
539 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
543 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
544 ir->bContinuation && ir->ld_seed != -1)
546 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
552 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
553 CHECK(ir->ePBC != epbcXYZ);
554 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
555 CHECK(ir->ns_type != ensGRID);
556 sprintf(err_buf, "with TPI nstlist should be larger than zero");
557 CHECK(ir->nstlist <= 0);
558 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
559 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
560 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
561 CHECK(ir->cutoff_scheme == ecutsVERLET);
565 if ( (opts->nshake > 0) && (opts->bMorse) )
568 "Using morse bond-potentials while constraining bonds is useless");
569 warning(wi, warn_buf);
572 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
573 ir->bContinuation && ir->ld_seed != -1)
575 warning_note(wi, "You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
577 /* verify simulated tempering options */
581 gmx_bool bAllTempZero = TRUE;
582 for (i = 0; i < fep->n_lambda; i++)
584 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[efptTEMPERATURE], fep->all_lambda[efptTEMPERATURE][i]);
585 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
586 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
588 bAllTempZero = FALSE;
591 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
592 CHECK(bAllTempZero == TRUE);
594 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
595 CHECK(ir->eI != eiVV);
597 /* check compatability of the temperature coupling with simulated tempering */
599 if (ir->etc == etcNOSEHOOVER)
601 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
602 warning_note(wi, warn_buf);
605 /* check that the temperatures make sense */
607 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= than the simulated tempering lower temperature (%g)", ir->simtempvals->simtemp_high, ir->simtempvals->simtemp_low);
608 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
610 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
611 CHECK(ir->simtempvals->simtemp_high <= 0);
613 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
614 CHECK(ir->simtempvals->simtemp_low <= 0);
617 /* verify free energy options */
619 if (ir->efep != efepNO)
622 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
624 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
626 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
627 (int)fep->sc_r_power);
628 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
630 sprintf(err_buf, "Can't use positive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
631 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
633 sprintf(err_buf, "Can't use positive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
634 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
636 sprintf(err_buf, "Can only use expanded ensemble with md-vv (for now)");
637 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
639 sprintf(err_buf, "Free-energy not implemented for Ewald");
640 CHECK(ir->coulombtype == eelEWALD);
642 /* check validty of lambda inputs */
643 if (fep->n_lambda == 0)
645 /* Clear output in case of no states:*/
646 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
647 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
651 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
652 CHECK((fep->init_fep_state >= fep->n_lambda));
655 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
656 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
658 sprintf(err_buf, "init-lambda=%g while init-lambda-state=%d. Lambda state must be set either with init-lambda-state or with init-lambda, but not both",
659 fep->init_lambda, fep->init_fep_state);
660 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
664 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
668 for (i = 0; i < efptNR; i++)
670 if (fep->separate_dvdl[i])
675 if (n_lambda_terms > 1)
677 sprintf(warn_buf, "If lambda vector states (fep-lambdas, coul-lambdas etc.) are set, don't use init-lambda to set lambda state (except for slow growth). Use init-lambda-state instead.");
678 warning(wi, warn_buf);
681 if (n_lambda_terms < 2 && fep->n_lambda > 0)
684 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
688 for (j = 0; j < efptNR; j++)
690 for (i = 0; i < fep->n_lambda; i++)
692 sprintf(err_buf, "Entry %d for %s must be between 0 and 1, instead is %g", i, efpt_names[j], fep->all_lambda[j][i]);
693 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
697 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
699 for (i = 0; i < fep->n_lambda; i++)
701 sprintf(err_buf, "For state %d, vdw-lambdas (%f) is changing with vdw softcore, while coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to crashes, and is not supported.", i, fep->all_lambda[efptVDW][i],
702 fep->all_lambda[efptCOUL][i]);
703 CHECK((fep->sc_alpha > 0) &&
704 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
705 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
706 ((fep->all_lambda[efptVDW][i] > 0.0) &&
707 (fep->all_lambda[efptVDW][i] < 1.0))));
711 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
713 real sigma, lambda, r_sc;
716 /* Maximum estimate for A and B charges equal with lambda power 1 */
718 r_sc = std::pow(lambda*fep->sc_alpha*std::pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
719 sprintf(warn_buf, "With PME there is a minor soft core effect present at the cut-off, proportional to (LJsigma/rcoulomb)^%g. This could have a minor effect on energy conservation, but usually other effects dominate. With a common sigma value of %g nm the fraction of the particle-particle potential at the cut-off at lambda=%g is around %.1e, while ewald-rtol is %.1e.",
721 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
722 warning_note(wi, warn_buf);
725 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
726 be treated differently, but that's the next step */
728 for (i = 0; i < efptNR; i++)
730 for (j = 0; j < fep->n_lambda; j++)
732 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
733 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
738 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
742 /* checking equilibration of weights inputs for validity */
744 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
745 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
746 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
748 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
749 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
750 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
752 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
753 expand->equil_steps, elmceq_names[elmceqSTEPS]);
754 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
756 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
757 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
758 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
760 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
761 expand->equil_ratio, elmceq_names[elmceqRATIO]);
762 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
764 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
765 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
766 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
768 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
769 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
770 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
772 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
773 expand->equil_steps, elmceq_names[elmceqSTEPS]);
774 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
776 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
777 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
778 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
780 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
781 expand->equil_ratio, elmceq_names[elmceqRATIO]);
782 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
784 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
785 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
786 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
788 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
789 CHECK((expand->lmc_repeats <= 0));
790 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
791 CHECK((expand->minvarmin <= 0));
792 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
793 CHECK((expand->c_range < 0));
794 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
795 fep->init_fep_state, expand->lmc_forced_nstart);
796 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
797 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
798 CHECK((expand->lmc_forced_nstart < 0));
799 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
800 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
802 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
803 CHECK((expand->init_wl_delta < 0));
804 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
805 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
806 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
807 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
809 /* if there is no temperature control, we need to specify an MC temperature */
810 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
811 if (expand->nstTij > 0)
813 sprintf(err_buf, "nstlog must be non-zero");
814 CHECK(ir->nstlog != 0);
815 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
816 expand->nstTij, ir->nstlog);
817 CHECK((expand->nstTij % ir->nstlog) != 0);
822 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
823 CHECK(ir->nwall && ir->ePBC != epbcXY);
826 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
828 if (ir->ePBC == epbcNONE)
830 if (ir->epc != epcNO)
832 warning(wi, "Turning off pressure coupling for vacuum system");
838 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
839 epbc_names[ir->ePBC]);
840 CHECK(ir->epc != epcNO);
842 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
843 CHECK(EEL_FULL(ir->coulombtype));
845 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
846 epbc_names[ir->ePBC]);
847 CHECK(ir->eDispCorr != edispcNO);
850 if (ir->rlist == 0.0)
852 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
853 "with coulombtype = %s or coulombtype = %s\n"
854 "without periodic boundary conditions (pbc = %s) and\n"
855 "rcoulomb and rvdw set to zero",
856 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
857 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
858 (ir->ePBC != epbcNONE) ||
859 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
863 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
868 if (ir->nstcomm == 0)
870 ir->comm_mode = ecmNO;
872 if (ir->comm_mode != ecmNO)
876 warning(wi, "If you want to remove the rotation around the center of mass, you should set comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to its absolute value");
877 ir->nstcomm = abs(ir->nstcomm);
880 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
882 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
883 ir->nstcomm = ir->nstcalcenergy;
886 if (ir->comm_mode == ecmANGULAR)
888 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
889 CHECK(ir->bPeriodicMols);
890 if (ir->ePBC != epbcNONE)
892 warning(wi, "Removing the rotation around the center of mass in a periodic system, this can lead to artifacts. Only use this on a single (cluster of) molecules. This cluster should not cross periodic boundaries.");
897 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
899 warning_note(wi, "Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
902 /* TEMPERATURE COUPLING */
903 if (ir->etc == etcYES)
905 ir->etc = etcBERENDSEN;
906 warning_note(wi, "Old option for temperature coupling given: "
907 "changing \"yes\" to \"Berendsen\"\n");
910 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
912 if (ir->opts.nhchainlength < 1)
914 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
915 ir->opts.nhchainlength = 1;
916 warning(wi, warn_buf);
919 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
921 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
922 ir->opts.nhchainlength = 1;
927 ir->opts.nhchainlength = 0;
930 if (ir->eI == eiVVAK)
932 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
934 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
937 if (ETC_ANDERSEN(ir->etc))
939 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
940 CHECK(!(EI_VV(ir->eI)));
942 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
944 sprintf(warn_buf, "Center of mass removal not necessary for %s. All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.", etcoupl_names[ir->etc]);
945 warning_note(wi, warn_buf);
948 sprintf(err_buf, "nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step", ir->nstcomm, etcoupl_names[ir->etc]);
949 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
952 if (ir->etc == etcBERENDSEN)
954 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
955 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
956 warning_note(wi, warn_buf);
959 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
960 && ir->epc == epcBERENDSEN)
962 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
963 "true ensemble for the thermostat");
964 warning(wi, warn_buf);
967 /* PRESSURE COUPLING */
968 if (ir->epc == epcISOTROPIC)
970 ir->epc = epcBERENDSEN;
971 warning_note(wi, "Old option for pressure coupling given: "
972 "changing \"Isotropic\" to \"Berendsen\"\n");
975 if (ir->epc != epcNO)
977 dt_pcoupl = ir->nstpcouple*ir->delta_t;
979 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
980 CHECK(ir->tau_p <= 0);
982 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
984 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
985 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
986 warning(wi, warn_buf);
989 sprintf(err_buf, "compressibility must be > 0 when using pressure"
990 " coupling %s\n", EPCOUPLTYPE(ir->epc));
991 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
992 ir->compress[ZZ][ZZ] < 0 ||
993 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
994 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
996 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
999 "You are generating velocities so I am assuming you "
1000 "are equilibrating a system. You are using "
1001 "%s pressure coupling, but this can be "
1002 "unstable for equilibration. If your system crashes, try "
1003 "equilibrating first with Berendsen pressure coupling. If "
1004 "you are not equilibrating the system, you can probably "
1005 "ignore this warning.",
1006 epcoupl_names[ir->epc]);
1007 warning(wi, warn_buf);
1013 if (ir->epc > epcNO)
1015 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1017 warning_error(wi, "for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
1023 if (ir->epc == epcMTTK)
1025 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1029 /* ELECTROSTATICS */
1030 /* More checks are in triple check (grompp.c) */
1032 if (ir->coulombtype == eelSWITCH)
1034 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1035 "artifacts, advice: use coulombtype = %s",
1036 eel_names[ir->coulombtype],
1037 eel_names[eelRF_ZERO]);
1038 warning(wi, warn_buf);
1041 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1043 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1044 warning_note(wi, warn_buf);
1047 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1049 sprintf(warn_buf, "epsilon-r = %g and epsilon-rf = 1 with reaction field, proceeding assuming old format and exchanging epsilon-r and epsilon-rf", ir->epsilon_r);
1050 warning(wi, warn_buf);
1051 ir->epsilon_rf = ir->epsilon_r;
1052 ir->epsilon_r = 1.0;
1055 if (ir->epsilon_r == 0)
1058 "It is pointless to use long-range or Generalized Born electrostatics with infinite relative permittivity."
1059 "Since you are effectively turning of electrostatics, a plain cutoff will be much faster.");
1060 CHECK(EEL_FULL(ir->coulombtype) || ir->implicit_solvent == eisGBSA);
1063 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == NULL)
1065 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1066 CHECK(ir->epsilon_r < 0);
1069 if (EEL_RF(ir->coulombtype))
1071 /* reaction field (at the cut-off) */
1073 if (ir->coulombtype == eelRF_ZERO && ir->epsilon_rf != 0)
1075 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1076 eel_names[ir->coulombtype]);
1077 warning(wi, warn_buf);
1078 ir->epsilon_rf = 0.0;
1081 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1082 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1083 (ir->epsilon_r == 0));
1084 if (ir->epsilon_rf == ir->epsilon_r)
1086 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1087 eel_names[ir->coulombtype]);
1088 warning(wi, warn_buf);
1091 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1092 * means the interaction is zero outside rcoulomb, but it helps to
1093 * provide accurate energy conservation.
1095 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1097 if (ir_coulomb_switched(ir))
1100 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1101 eel_names[ir->coulombtype]);
1102 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1105 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1107 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1109 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1110 eel_names[ir->coulombtype]);
1111 CHECK(ir->rlist > ir->rcoulomb);
1115 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1118 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1119 CHECK( ir->coulomb_modifier != eintmodNONE);
1121 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1124 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1125 CHECK( ir->vdw_modifier != eintmodNONE);
1128 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1129 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1132 "The switch/shift interaction settings are just for compatibility; you will get better "
1133 "performance from applying potential modifiers to your interactions!\n");
1134 warning_note(wi, warn_buf);
1137 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1139 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1141 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1142 sprintf(warn_buf, "The switching range should be 5%% or less (currently %.2f%% using a switching range of %4f-%4f) for accurate electrostatic energies, energy conservation will be good regardless, since ewald_rtol = %g.",
1143 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1144 warning(wi, warn_buf);
1148 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1150 if (ir->rvdw_switch == 0)
1152 sprintf(warn_buf, "rvdw-switch is equal 0 even though you are using a switched Lennard-Jones potential. This suggests it was not set in the mdp, which can lead to large energy errors. In GROMACS, 0.05 to 0.1 nm is often a reasonable vdw switching range.");
1153 warning(wi, warn_buf);
1157 if (EEL_FULL(ir->coulombtype))
1159 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1160 ir->coulombtype == eelPMEUSERSWITCH)
1162 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1163 eel_names[ir->coulombtype]);
1164 CHECK(ir->rcoulomb > ir->rlist);
1166 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1168 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1171 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist.\n"
1172 "For optimal energy conservation,consider using\n"
1173 "a potential modifier.", eel_names[ir->coulombtype]);
1174 CHECK(ir->rcoulomb != ir->rlist);
1179 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1181 if (ir->pme_order < 3)
1183 warning_error(wi, "pme-order can not be smaller than 3");
1187 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1189 if (ir->ewald_geometry == eewg3D)
1191 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1192 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1193 warning(wi, warn_buf);
1195 /* This check avoids extra pbc coding for exclusion corrections */
1196 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1197 CHECK(ir->wall_ewald_zfac < 2);
1199 if ((ir->ewald_geometry == eewg3DC) && (ir->ePBC != epbcXY) &&
1200 EEL_FULL(ir->coulombtype))
1202 sprintf(warn_buf, "With %s and ewald_geometry = %s you should use pbc = %s",
1203 eel_names[ir->coulombtype], eewg_names[eewg3DC], epbc_names[epbcXY]);
1204 warning(wi, warn_buf);
1206 if ((ir->epsilon_surface != 0) && EEL_FULL(ir->coulombtype))
1208 if (ir->cutoff_scheme == ecutsVERLET)
1210 sprintf(warn_buf, "Since molecules/charge groups are broken using the Verlet scheme, you can not use a dipole correction to the %s electrostatics.",
1211 eel_names[ir->coulombtype]);
1212 warning(wi, warn_buf);
1216 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",
1217 eel_names[ir->coulombtype]);
1218 warning_note(wi, warn_buf);
1222 if (ir_vdw_switched(ir))
1224 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1225 CHECK(ir->rvdw_switch >= ir->rvdw);
1227 if (ir->rvdw_switch < 0.5*ir->rvdw)
1229 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.",
1230 ir->rvdw_switch, ir->rvdw);
1231 warning_note(wi, warn_buf);
1234 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1236 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1238 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1239 CHECK(ir->rlist > ir->rvdw);
1243 if (ir->vdwtype == evdwPME)
1245 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1247 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s and %s",
1248 evdw_names[ir->vdwtype],
1249 eintmod_names[eintmodPOTSHIFT],
1250 eintmod_names[eintmodNONE]);
1254 if (ir->cutoff_scheme == ecutsGROUP)
1256 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1257 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)))
1259 warning_note(wi, "With exact cut-offs, rlist should be "
1260 "larger than rcoulomb and rvdw, so that there "
1261 "is a buffer region for particle motion "
1262 "between neighborsearch steps");
1265 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlist <= ir->rcoulomb)
1267 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rcoulomb.");
1268 warning_note(wi, warn_buf);
1270 if (ir_vdw_switched(ir) && (ir->rlist <= ir->rvdw))
1272 sprintf(warn_buf, "For energy conservation with switch/shift potentials, rlist should be 0.1 to 0.3 nm larger than rvdw.");
1273 warning_note(wi, warn_buf);
1277 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1279 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.");
1282 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1285 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1288 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1290 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1293 /* ENERGY CONSERVATION */
1294 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1296 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1298 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1299 evdw_names[evdwSHIFT]);
1300 warning_note(wi, warn_buf);
1302 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1304 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1305 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1306 warning_note(wi, warn_buf);
1310 /* IMPLICIT SOLVENT */
1311 if (ir->coulombtype == eelGB_NOTUSED)
1313 sprintf(warn_buf, "Invalid option %s for coulombtype",
1314 eel_names[ir->coulombtype]);
1315 warning_error(wi, warn_buf);
1318 if (ir->sa_algorithm == esaSTILL)
1320 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1321 CHECK(ir->sa_algorithm == esaSTILL);
1324 if (ir->implicit_solvent == eisGBSA)
1326 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1327 CHECK(ir->rgbradii != ir->rlist);
1329 if (ir->coulombtype != eelCUT)
1331 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1332 CHECK(ir->coulombtype != eelCUT);
1334 if (ir->vdwtype != evdwCUT)
1336 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1337 CHECK(ir->vdwtype != evdwCUT);
1339 if (ir->nstgbradii < 1)
1341 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1342 warning_note(wi, warn_buf);
1345 if (ir->sa_algorithm == esaNO)
1347 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1348 warning_note(wi, warn_buf);
1350 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1352 sprintf(warn_buf, "Value of sa_surface_tension is < 0. Changing it to 2.05016 or 2.25936 kJ/nm^2/mol for Still and HCT/OBC respectively\n");
1353 warning_note(wi, warn_buf);
1355 if (ir->gb_algorithm == egbSTILL)
1357 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1361 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1364 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1366 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1367 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1374 gmx_fatal(FARGS, "AdResS simulations are no longer supported");
1378 /* count the number of text elemets separated by whitespace in a string.
1379 str = the input string
1380 maxptr = the maximum number of allowed elements
1381 ptr = the output array of pointers to the first character of each element
1382 returns: the number of elements. */
1383 int str_nelem(const char *str, int maxptr, char *ptr[])
1388 copy0 = gmx_strdup(str);
1391 while (*copy != '\0')
1395 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1403 while ((*copy != '\0') && !isspace(*copy))
1422 /* interpret a number of doubles from a string and put them in an array,
1423 after allocating space for them.
1424 str = the input string
1425 n = the (pre-allocated) number of doubles read
1426 r = the output array of doubles. */
1427 static void parse_n_real(char *str, int *n, real **r, warninp_t wi)
1432 char warn_buf[STRLEN];
1434 *n = str_nelem(str, MAXPTR, ptr);
1437 for (i = 0; i < *n; i++)
1439 (*r)[i] = strtod(ptr[i], &endptr);
1442 sprintf(warn_buf, "Invalid value %s in string in mdp file. Expected a real number.", ptr[i]);
1443 warning_error(wi, warn_buf);
1448 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN], warninp_t wi)
1451 int i, j, max_n_lambda, nweights, nfep[efptNR];
1452 t_lambda *fep = ir->fepvals;
1453 t_expanded *expand = ir->expandedvals;
1454 real **count_fep_lambdas;
1455 gmx_bool bOneLambda = TRUE;
1457 snew(count_fep_lambdas, efptNR);
1459 /* FEP input processing */
1460 /* first, identify the number of lambda values for each type.
1461 All that are nonzero must have the same number */
1463 for (i = 0; i < efptNR; i++)
1465 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]), wi);
1468 /* now, determine the number of components. All must be either zero, or equal. */
1471 for (i = 0; i < efptNR; i++)
1473 if (nfep[i] > max_n_lambda)
1475 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1476 must have the same number if its not zero.*/
1481 for (i = 0; i < efptNR; i++)
1485 ir->fepvals->separate_dvdl[i] = FALSE;
1487 else if (nfep[i] == max_n_lambda)
1489 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1490 respect to the temperature currently */
1492 ir->fepvals->separate_dvdl[i] = TRUE;
1497 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1498 nfep[i], efpt_names[i], max_n_lambda);
1501 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1502 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1504 /* the number of lambdas is the number we've read in, which is either zero
1505 or the same for all */
1506 fep->n_lambda = max_n_lambda;
1508 /* allocate space for the array of lambda values */
1509 snew(fep->all_lambda, efptNR);
1510 /* if init_lambda is defined, we need to set lambda */
1511 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1513 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1515 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1516 for (i = 0; i < efptNR; i++)
1518 snew(fep->all_lambda[i], fep->n_lambda);
1519 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1522 for (j = 0; j < fep->n_lambda; j++)
1524 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1526 sfree(count_fep_lambdas[i]);
1529 sfree(count_fep_lambdas);
1531 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1532 bookkeeping -- for now, init_lambda */
1534 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1536 for (i = 0; i < fep->n_lambda; i++)
1538 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1542 /* check to see if only a single component lambda is defined, and soft core is defined.
1543 In this case, turn on coulomb soft core */
1545 if (max_n_lambda == 0)
1551 for (i = 0; i < efptNR; i++)
1553 if ((nfep[i] != 0) && (i != efptFEP))
1559 if ((bOneLambda) && (fep->sc_alpha > 0))
1561 fep->bScCoul = TRUE;
1564 /* Fill in the others with the efptFEP if they are not explicitly
1565 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1566 they are all zero. */
1568 for (i = 0; i < efptNR; i++)
1570 if ((nfep[i] == 0) && (i != efptFEP))
1572 for (j = 0; j < fep->n_lambda; j++)
1574 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1580 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1581 if (fep->sc_r_power == 48)
1583 if (fep->sc_alpha > 0.1)
1585 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1589 /* now read in the weights */
1590 parse_n_real(weights, &nweights, &(expand->init_lambda_weights), wi);
1593 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1595 else if (nweights != fep->n_lambda)
1597 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1598 nweights, fep->n_lambda);
1600 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1602 expand->nstexpanded = fep->nstdhdl;
1603 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1605 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1607 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1608 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1609 2*tau_t just to be careful so it's not to frequent */
1614 static void do_simtemp_params(t_inputrec *ir)
1617 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1618 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1623 static void do_wall_params(t_inputrec *ir,
1624 char *wall_atomtype, char *wall_density,
1628 char *names[MAXPTR];
1631 opts->wall_atomtype[0] = NULL;
1632 opts->wall_atomtype[1] = NULL;
1634 ir->wall_atomtype[0] = -1;
1635 ir->wall_atomtype[1] = -1;
1636 ir->wall_density[0] = 0;
1637 ir->wall_density[1] = 0;
1641 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1642 if (nstr != ir->nwall)
1644 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1647 for (i = 0; i < ir->nwall; i++)
1649 opts->wall_atomtype[i] = gmx_strdup(names[i]);
1652 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1654 nstr = str_nelem(wall_density, MAXPTR, names);
1655 if (nstr != ir->nwall)
1657 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1659 for (i = 0; i < ir->nwall; i++)
1661 sscanf(names[i], "%lf", &dbl);
1664 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1666 ir->wall_density[i] = dbl;
1672 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1680 srenew(groups->grpname, groups->ngrpname+nwall);
1681 grps = &(groups->grps[egcENER]);
1682 srenew(grps->nm_ind, grps->nr+nwall);
1683 for (i = 0; i < nwall; i++)
1685 sprintf(str, "wall%d", i);
1686 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1687 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1692 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1693 t_expanded *expand, warninp_t wi)
1701 /* read expanded ensemble parameters */
1702 CCTYPE ("expanded ensemble variables");
1703 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1704 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1705 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1706 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1707 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1708 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1709 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1710 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1711 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1712 CCTYPE("Seed for Monte Carlo in lambda space");
1713 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1714 RTYPE ("mc-temperature", expand->mc_temp, -1);
1715 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1716 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1717 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1718 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1719 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1720 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1721 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1722 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1723 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1724 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1725 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1733 /*! \brief Return whether an end state with the given coupling-lambda
1734 * value describes fully-interacting VDW.
1736 * \param[in] couple_lambda_value Enumeration ecouplam value describing the end state
1737 * \return Whether VDW is on (i.e. the user chose vdw or vdw-q in the .mdp file)
1739 static gmx_bool couple_lambda_has_vdw_on(int couple_lambda_value)
1741 return (couple_lambda_value == ecouplamVDW ||
1742 couple_lambda_value == ecouplamVDWQ);
1745 void get_ir(const char *mdparin, const char *mdparout,
1746 t_inputrec *ir, t_gromppopts *opts,
1750 double dumdub[2][6];
1754 char warn_buf[STRLEN];
1755 t_lambda *fep = ir->fepvals;
1756 t_expanded *expand = ir->expandedvals;
1758 init_inputrec_strings();
1759 inp = read_inpfile(mdparin, &ninp, wi);
1761 snew(dumstr[0], STRLEN);
1762 snew(dumstr[1], STRLEN);
1764 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1767 "%s did not specify a value for the .mdp option "
1768 "\"cutoff-scheme\". Probably it was first intended for use "
1769 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1770 "introduced, but the group scheme was still the default. "
1771 "The default is now the Verlet scheme, so you will observe "
1772 "different behaviour.", mdparin);
1773 warning_note(wi, warn_buf);
1776 /* ignore the following deprecated commands */
1779 REM_TYPE("domain-decomposition");
1780 REM_TYPE("andersen-seed");
1782 REM_TYPE("dihre-fc");
1783 REM_TYPE("dihre-tau");
1784 REM_TYPE("nstdihreout");
1785 REM_TYPE("nstcheckpoint");
1786 REM_TYPE("optimize-fft");
1787 REM_TYPE("adress_type");
1788 REM_TYPE("adress_const_wf");
1789 REM_TYPE("adress_ex_width");
1790 REM_TYPE("adress_hy_width");
1791 REM_TYPE("adress_ex_forcecap");
1792 REM_TYPE("adress_interface_correction");
1793 REM_TYPE("adress_site");
1794 REM_TYPE("adress_reference_coords");
1795 REM_TYPE("adress_tf_grp_names");
1796 REM_TYPE("adress_cg_grp_names");
1797 REM_TYPE("adress_do_hybridpairs");
1798 REM_TYPE("rlistlong");
1799 REM_TYPE("nstcalclr");
1800 REM_TYPE("pull-print-com2");
1802 /* replace the following commands with the clearer new versions*/
1803 REPL_TYPE("unconstrained-start", "continuation");
1804 REPL_TYPE("foreign-lambda", "fep-lambdas");
1805 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1806 REPL_TYPE("nstxtcout", "nstxout-compressed");
1807 REPL_TYPE("xtc-grps", "compressed-x-grps");
1808 REPL_TYPE("xtc-precision", "compressed-x-precision");
1809 REPL_TYPE("pull-print-com1", "pull-print-com");
1811 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1812 CTYPE ("Preprocessor information: use cpp syntax.");
1813 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1814 STYPE ("include", opts->include, NULL);
1815 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1816 STYPE ("define", opts->define, NULL);
1818 CCTYPE ("RUN CONTROL PARAMETERS");
1819 EETYPE("integrator", ir->eI, ei_names);
1820 CTYPE ("Start time and timestep in ps");
1821 RTYPE ("tinit", ir->init_t, 0.0);
1822 RTYPE ("dt", ir->delta_t, 0.001);
1823 STEPTYPE ("nsteps", ir->nsteps, 0);
1824 CTYPE ("For exact run continuation or redoing part of a run");
1825 STEPTYPE ("init-step", ir->init_step, 0);
1826 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1827 ITYPE ("simulation-part", ir->simulation_part, 1);
1828 CTYPE ("mode for center of mass motion removal");
1829 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1830 CTYPE ("number of steps for center of mass motion removal");
1831 ITYPE ("nstcomm", ir->nstcomm, 100);
1832 CTYPE ("group(s) for center of mass motion removal");
1833 STYPE ("comm-grps", is->vcm, NULL);
1835 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1836 CTYPE ("Friction coefficient (amu/ps) and random seed");
1837 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1838 STEPTYPE ("ld-seed", ir->ld_seed, -1);
1841 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1842 CTYPE ("Force tolerance and initial step-size");
1843 RTYPE ("emtol", ir->em_tol, 10.0);
1844 RTYPE ("emstep", ir->em_stepsize, 0.01);
1845 CTYPE ("Max number of iterations in relax-shells");
1846 ITYPE ("niter", ir->niter, 20);
1847 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1848 RTYPE ("fcstep", ir->fc_stepsize, 0);
1849 CTYPE ("Frequency of steepest descents steps when doing CG");
1850 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1851 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1853 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1854 RTYPE ("rtpi", ir->rtpi, 0.05);
1856 /* Output options */
1857 CCTYPE ("OUTPUT CONTROL OPTIONS");
1858 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1859 ITYPE ("nstxout", ir->nstxout, 0);
1860 ITYPE ("nstvout", ir->nstvout, 0);
1861 ITYPE ("nstfout", ir->nstfout, 0);
1862 CTYPE ("Output frequency for energies to log file and energy file");
1863 ITYPE ("nstlog", ir->nstlog, 1000);
1864 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1865 ITYPE ("nstenergy", ir->nstenergy, 1000);
1866 CTYPE ("Output frequency and precision for .xtc file");
1867 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
1868 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1869 CTYPE ("This selects the subset of atoms for the compressed");
1870 CTYPE ("trajectory file. You can select multiple groups. By");
1871 CTYPE ("default, all atoms will be written.");
1872 STYPE ("compressed-x-grps", is->x_compressed_groups, NULL);
1873 CTYPE ("Selection of energy groups");
1874 STYPE ("energygrps", is->energy, NULL);
1876 /* Neighbor searching */
1877 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1878 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1879 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1880 CTYPE ("nblist update frequency");
1881 ITYPE ("nstlist", ir->nstlist, 10);
1882 CTYPE ("ns algorithm (simple or grid)");
1883 EETYPE("ns-type", ir->ns_type, ens_names);
1884 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1885 EETYPE("pbc", ir->ePBC, epbc_names);
1886 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1887 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1888 CTYPE ("a value of -1 means: use rlist");
1889 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1890 CTYPE ("nblist cut-off");
1891 RTYPE ("rlist", ir->rlist, 1.0);
1892 CTYPE ("long-range cut-off for switched potentials");
1894 /* Electrostatics */
1895 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1896 CTYPE ("Method for doing electrostatics");
1897 EETYPE("coulombtype", ir->coulombtype, eel_names);
1898 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1899 CTYPE ("cut-off lengths");
1900 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1901 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1902 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1903 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1904 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1905 CTYPE ("Method for doing Van der Waals");
1906 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1907 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1908 CTYPE ("cut-off lengths");
1909 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1910 RTYPE ("rvdw", ir->rvdw, 1.0);
1911 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1912 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1913 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1914 RTYPE ("table-extension", ir->tabext, 1.0);
1915 CTYPE ("Separate tables between energy group pairs");
1916 STYPE ("energygrp-table", is->egptable, NULL);
1917 CTYPE ("Spacing for the PME/PPPM FFT grid");
1918 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1919 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1920 ITYPE ("fourier-nx", ir->nkx, 0);
1921 ITYPE ("fourier-ny", ir->nky, 0);
1922 ITYPE ("fourier-nz", ir->nkz, 0);
1923 CTYPE ("EWALD/PME/PPPM parameters");
1924 ITYPE ("pme-order", ir->pme_order, 4);
1925 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1926 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1927 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1928 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1929 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1931 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1932 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1934 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1935 CTYPE ("Algorithm for calculating Born radii");
1936 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1937 CTYPE ("Frequency of calculating the Born radii inside rlist");
1938 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1939 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1940 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1941 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1942 CTYPE ("Dielectric coefficient of the implicit solvent");
1943 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1944 CTYPE ("Salt concentration in M for Generalized Born models");
1945 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1946 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1947 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1948 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1949 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1950 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1951 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1952 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1953 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1954 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1956 /* Coupling stuff */
1957 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1958 CTYPE ("Temperature coupling");
1959 EETYPE("tcoupl", ir->etc, etcoupl_names);
1960 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1961 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
1962 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1963 CTYPE ("Groups to couple separately");
1964 STYPE ("tc-grps", is->tcgrps, NULL);
1965 CTYPE ("Time constant (ps) and reference temperature (K)");
1966 STYPE ("tau-t", is->tau_t, NULL);
1967 STYPE ("ref-t", is->ref_t, NULL);
1968 CTYPE ("pressure coupling");
1969 EETYPE("pcoupl", ir->epc, epcoupl_names);
1970 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1971 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1972 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1973 RTYPE ("tau-p", ir->tau_p, 1.0);
1974 STYPE ("compressibility", dumstr[0], NULL);
1975 STYPE ("ref-p", dumstr[1], NULL);
1976 CTYPE ("Scaling of reference coordinates, No, All or COM");
1977 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1980 CCTYPE ("OPTIONS FOR QMMM calculations");
1981 EETYPE("QMMM", ir->bQMMM, yesno_names);
1982 CTYPE ("Groups treated Quantum Mechanically");
1983 STYPE ("QMMM-grps", is->QMMM, NULL);
1984 CTYPE ("QM method");
1985 STYPE("QMmethod", is->QMmethod, NULL);
1986 CTYPE ("QMMM scheme");
1987 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1988 CTYPE ("QM basisset");
1989 STYPE("QMbasis", is->QMbasis, NULL);
1990 CTYPE ("QM charge");
1991 STYPE ("QMcharge", is->QMcharge, NULL);
1992 CTYPE ("QM multiplicity");
1993 STYPE ("QMmult", is->QMmult, NULL);
1994 CTYPE ("Surface Hopping");
1995 STYPE ("SH", is->bSH, NULL);
1996 CTYPE ("CAS space options");
1997 STYPE ("CASorbitals", is->CASorbitals, NULL);
1998 STYPE ("CASelectrons", is->CASelectrons, NULL);
1999 STYPE ("SAon", is->SAon, NULL);
2000 STYPE ("SAoff", is->SAoff, NULL);
2001 STYPE ("SAsteps", is->SAsteps, NULL);
2002 CTYPE ("Scale factor for MM charges");
2003 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
2004 CTYPE ("Optimization of QM subsystem");
2005 STYPE ("bOPT", is->bOPT, NULL);
2006 STYPE ("bTS", is->bTS, NULL);
2008 /* Simulated annealing */
2009 CCTYPE("SIMULATED ANNEALING");
2010 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
2011 STYPE ("annealing", is->anneal, NULL);
2012 CTYPE ("Number of time points to use for specifying annealing in each group");
2013 STYPE ("annealing-npoints", is->anneal_npoints, NULL);
2014 CTYPE ("List of times at the annealing points for each group");
2015 STYPE ("annealing-time", is->anneal_time, NULL);
2016 CTYPE ("Temp. at each annealing point, for each group.");
2017 STYPE ("annealing-temp", is->anneal_temp, NULL);
2020 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
2021 EETYPE("gen-vel", opts->bGenVel, yesno_names);
2022 RTYPE ("gen-temp", opts->tempi, 300.0);
2023 ITYPE ("gen-seed", opts->seed, -1);
2026 CCTYPE ("OPTIONS FOR BONDS");
2027 EETYPE("constraints", opts->nshake, constraints);
2028 CTYPE ("Type of constraint algorithm");
2029 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
2030 CTYPE ("Do not constrain the start configuration");
2031 EETYPE("continuation", ir->bContinuation, yesno_names);
2032 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2033 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2034 CTYPE ("Relative tolerance of shake");
2035 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2036 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2037 ITYPE ("lincs-order", ir->nProjOrder, 4);
2038 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2039 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2040 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2041 ITYPE ("lincs-iter", ir->nLincsIter, 1);
2042 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2043 CTYPE ("rotates over more degrees than");
2044 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2045 CTYPE ("Convert harmonic bonds to morse potentials");
2046 EETYPE("morse", opts->bMorse, yesno_names);
2048 /* Energy group exclusions */
2049 CCTYPE ("ENERGY GROUP EXCLUSIONS");
2050 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2051 STYPE ("energygrp-excl", is->egpexcl, NULL);
2055 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2056 ITYPE ("nwall", ir->nwall, 0);
2057 EETYPE("wall-type", ir->wall_type, ewt_names);
2058 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2059 STYPE ("wall-atomtype", is->wall_atomtype, NULL);
2060 STYPE ("wall-density", is->wall_density, NULL);
2061 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2064 CCTYPE("COM PULLING");
2065 EETYPE("pull", ir->bPull, yesno_names);
2069 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, wi);
2072 /* Enforced rotation */
2073 CCTYPE("ENFORCED ROTATION");
2074 CTYPE("Enforced rotation: No or Yes");
2075 EETYPE("rotation", ir->bRot, yesno_names);
2079 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2082 /* Interactive MD */
2084 CCTYPE("Group to display and/or manipulate in interactive MD session");
2085 STYPE ("IMD-group", is->imd_grp, NULL);
2086 if (is->imd_grp[0] != '\0')
2093 CCTYPE("NMR refinement stuff");
2094 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2095 EETYPE("disre", ir->eDisre, edisre_names);
2096 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2097 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2098 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2099 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2100 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2101 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2102 CTYPE ("Output frequency for pair distances to energy file");
2103 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2104 CTYPE ("Orientation restraints: No or Yes");
2105 EETYPE("orire", opts->bOrire, yesno_names);
2106 CTYPE ("Orientation restraints force constant and tau for time averaging");
2107 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2108 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2109 STYPE ("orire-fitgrp", is->orirefitgrp, NULL);
2110 CTYPE ("Output frequency for trace(SD) and S to energy file");
2111 ITYPE ("nstorireout", ir->nstorireout, 100);
2113 /* free energy variables */
2114 CCTYPE ("Free energy variables");
2115 EETYPE("free-energy", ir->efep, efep_names);
2116 STYPE ("couple-moltype", is->couple_moltype, NULL);
2117 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2118 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2119 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2121 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2123 it was not entered */
2124 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2125 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2126 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2127 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
2128 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
2129 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
2130 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
2131 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
2132 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
2133 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
2134 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2135 STYPE ("init-lambda-weights", is->lambda_weights, NULL);
2136 EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
2137 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2138 ITYPE ("sc-power", fep->sc_power, 1);
2139 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2140 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2141 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2142 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2143 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2144 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2145 separate_dhdl_file_names);
2146 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2147 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2148 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2150 /* Non-equilibrium MD stuff */
2151 CCTYPE("Non-equilibrium MD stuff");
2152 STYPE ("acc-grps", is->accgrps, NULL);
2153 STYPE ("accelerate", is->acc, NULL);
2154 STYPE ("freezegrps", is->freeze, NULL);
2155 STYPE ("freezedim", is->frdim, NULL);
2156 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2157 STYPE ("deform", is->deform, NULL);
2159 /* simulated tempering variables */
2160 CCTYPE("simulated tempering variables");
2161 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2162 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2163 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2164 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2166 /* expanded ensemble variables */
2167 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2169 read_expandedparams(&ninp, &inp, expand, wi);
2172 /* Electric fields */
2173 CCTYPE("Electric fields");
2174 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2175 CTYPE ("and a phase angle (real)");
2176 STYPE ("E-x", is->efield_x, NULL);
2177 CTYPE ("Time dependent (pulsed) electric field. Format is omega, time for pulse");
2178 CTYPE ("peak, and sigma (width) for pulse. Sigma = 0 removes pulse, leaving");
2179 CTYPE ("the field to be a cosine function.");
2180 STYPE ("E-xt", is->efield_xt, NULL);
2181 STYPE ("E-y", is->efield_y, NULL);
2182 STYPE ("E-yt", is->efield_yt, NULL);
2183 STYPE ("E-z", is->efield_z, NULL);
2184 STYPE ("E-zt", is->efield_zt, NULL);
2186 /* Ion/water position swapping ("computational electrophysiology") */
2187 CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2188 CTYPE("Swap positions along direction: no, X, Y, Z");
2189 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2190 if (ir->eSwapCoords != eswapNO)
2197 CTYPE("Swap attempt frequency");
2198 ITYPE("swap-frequency", ir->swap->nstswap, 1);
2199 CTYPE("Number of ion types to be controlled");
2200 ITYPE("iontypes", nIonTypes, 1);
2203 warning_error(wi, "You need to provide at least one ion type for position exchanges.");
2205 ir->swap->ngrp = nIonTypes + eSwapFixedGrpNR;
2206 snew(ir->swap->grp, ir->swap->ngrp);
2207 for (i = 0; i < ir->swap->ngrp; i++)
2209 snew(ir->swap->grp[i].molname, STRLEN);
2211 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2212 STYPE("split-group0", ir->swap->grp[eGrpSplit0].molname, NULL);
2213 STYPE("split-group1", ir->swap->grp[eGrpSplit1].molname, NULL);
2214 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2215 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2216 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2218 CTYPE("Name of solvent molecules");
2219 STYPE("solvent-group", ir->swap->grp[eGrpSolvent].molname, NULL);
2221 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2222 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2223 CTYPE("however, if correctly defined, the permeation events are recorded per channel");
2224 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2225 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2226 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2227 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2228 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2229 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2231 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2232 ITYPE("coupl-steps", ir->swap->nAverage, 10);
2234 CTYPE("Names of the ion types that can be exchanged with solvent molecules,");
2235 CTYPE("and the requested number of ions of this type in compartments A and B");
2236 CTYPE("-1 means fix the numbers as found in step 0");
2237 for (i = 0; i < nIonTypes; i++)
2239 int ig = eSwapFixedGrpNR + i;
2241 sprintf(buf, "iontype%d-name", i);
2242 STYPE(buf, ir->swap->grp[ig].molname, NULL);
2243 sprintf(buf, "iontype%d-in-A", i);
2244 ITYPE(buf, ir->swap->grp[ig].nmolReq[0], -1);
2245 sprintf(buf, "iontype%d-in-B", i);
2246 ITYPE(buf, ir->swap->grp[ig].nmolReq[1], -1);
2249 CTYPE("By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers");
2250 CTYPE("at maximum distance (= bulk concentration) to the split group layers. However,");
2251 CTYPE("an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0");
2252 CTYPE("towards one of the compartment-partitioning layers (at +/- 1.0).");
2253 RTYPE("bulk-offsetA", ir->swap->bulkOffset[0], 0.0);
2254 RTYPE("bulk-offsetB", ir->swap->bulkOffset[1], 0.0);
2255 if (!(ir->swap->bulkOffset[0] > -1.0 && ir->swap->bulkOffset[0] < 1.0)
2256 || !(ir->swap->bulkOffset[1] > -1.0 && ir->swap->bulkOffset[1] < 1.0) )
2258 warning_error(wi, "Bulk layer offsets must be > -1.0 and < 1.0 !");
2261 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2262 RTYPE("threshold", ir->swap->threshold, 1.0);
2265 /* AdResS is no longer supported, but we need mdrun to be able to refuse to run old AdResS .tpr files */
2266 EETYPE("adress", ir->bAdress, yesno_names);
2268 /* User defined thingies */
2269 CCTYPE ("User defined thingies");
2270 STYPE ("user1-grps", is->user1, NULL);
2271 STYPE ("user2-grps", is->user2, NULL);
2272 ITYPE ("userint1", ir->userint1, 0);
2273 ITYPE ("userint2", ir->userint2, 0);
2274 ITYPE ("userint3", ir->userint3, 0);
2275 ITYPE ("userint4", ir->userint4, 0);
2276 RTYPE ("userreal1", ir->userreal1, 0);
2277 RTYPE ("userreal2", ir->userreal2, 0);
2278 RTYPE ("userreal3", ir->userreal3, 0);
2279 RTYPE ("userreal4", ir->userreal4, 0);
2282 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2283 for (i = 0; (i < ninp); i++)
2286 sfree(inp[i].value);
2290 /* Process options if necessary */
2291 for (m = 0; m < 2; m++)
2293 for (i = 0; i < 2*DIM; i++)
2302 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2304 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2306 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2308 case epctSEMIISOTROPIC:
2309 case epctSURFACETENSION:
2310 if (sscanf(dumstr[m], "%lf%lf",
2311 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2313 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2315 dumdub[m][YY] = dumdub[m][XX];
2317 case epctANISOTROPIC:
2318 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2319 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2320 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2322 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2326 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2327 epcoupltype_names[ir->epct]);
2331 clear_mat(ir->ref_p);
2332 clear_mat(ir->compress);
2333 for (i = 0; i < DIM; i++)
2335 ir->ref_p[i][i] = dumdub[1][i];
2336 ir->compress[i][i] = dumdub[0][i];
2338 if (ir->epct == epctANISOTROPIC)
2340 ir->ref_p[XX][YY] = dumdub[1][3];
2341 ir->ref_p[XX][ZZ] = dumdub[1][4];
2342 ir->ref_p[YY][ZZ] = dumdub[1][5];
2343 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2345 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2347 ir->compress[XX][YY] = dumdub[0][3];
2348 ir->compress[XX][ZZ] = dumdub[0][4];
2349 ir->compress[YY][ZZ] = dumdub[0][5];
2350 for (i = 0; i < DIM; i++)
2352 for (m = 0; m < i; m++)
2354 ir->ref_p[i][m] = ir->ref_p[m][i];
2355 ir->compress[i][m] = ir->compress[m][i];
2360 if (ir->comm_mode == ecmNO)
2365 opts->couple_moltype = NULL;
2366 if (strlen(is->couple_moltype) > 0)
2368 if (ir->efep != efepNO)
2370 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2371 if (opts->couple_lam0 == opts->couple_lam1)
2373 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2375 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2376 opts->couple_lam1 == ecouplamNONE))
2378 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2383 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2386 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2387 if (ir->efep != efepNO)
2389 if (fep->delta_lambda > 0)
2391 ir->efep = efepSLOWGROWTH;
2395 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2397 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2398 warning_note(wi, "Old option for dhdl-print-energy given: "
2399 "changing \"yes\" to \"total\"\n");
2402 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2404 /* always print out the energy to dhdl if we are doing
2405 expanded ensemble, since we need the total energy for
2406 analysis if the temperature is changing. In some
2407 conditions one may only want the potential energy, so
2408 we will allow that if the appropriate mdp setting has
2409 been enabled. Otherwise, total it is:
2411 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2414 if ((ir->efep != efepNO) || ir->bSimTemp)
2416 ir->bExpanded = FALSE;
2417 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2419 ir->bExpanded = TRUE;
2421 do_fep_params(ir, is->fep_lambda, is->lambda_weights, wi);
2422 if (ir->bSimTemp) /* done after fep params */
2424 do_simtemp_params(ir);
2427 /* Because sc-coul (=FALSE by default) only acts on the lambda state
2428 * setup and not on the old way of specifying the free-energy setup,
2429 * we should check for using soft-core when not needed, since that
2430 * can complicate the sampling significantly.
2431 * Note that we only check for the automated coupling setup.
2432 * If the (advanced) user does FEP through manual topology changes,
2433 * this check will not be triggered.
2435 if (ir->efep != efepNO && ir->fepvals->n_lambda == 0 &&
2436 ir->fepvals->sc_alpha != 0 &&
2437 (couple_lambda_has_vdw_on(opts->couple_lam0) &&
2438 couple_lambda_has_vdw_on(opts->couple_lam1)))
2440 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.");
2445 ir->fepvals->n_lambda = 0;
2448 /* WALL PARAMETERS */
2450 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2452 /* ORIENTATION RESTRAINT PARAMETERS */
2454 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
2456 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2459 /* DEFORMATION PARAMETERS */
2461 clear_mat(ir->deform);
2462 for (i = 0; i < 6; i++)
2466 sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2467 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2468 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2469 for (i = 0; i < 3; i++)
2471 ir->deform[i][i] = dumdub[0][i];
2473 ir->deform[YY][XX] = dumdub[0][3];
2474 ir->deform[ZZ][XX] = dumdub[0][4];
2475 ir->deform[ZZ][YY] = dumdub[0][5];
2476 if (ir->epc != epcNO)
2478 for (i = 0; i < 3; i++)
2480 for (j = 0; j <= i; j++)
2482 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2484 warning_error(wi, "A box element has deform set and compressibility > 0");
2488 for (i = 0; i < 3; i++)
2490 for (j = 0; j < i; j++)
2492 if (ir->deform[i][j] != 0)
2494 for (m = j; m < DIM; m++)
2496 if (ir->compress[m][j] != 0)
2498 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.");
2499 warning(wi, warn_buf);
2507 /* Ion/water position swapping checks */
2508 if (ir->eSwapCoords != eswapNO)
2510 if (ir->swap->nstswap < 1)
2512 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2514 if (ir->swap->nAverage < 1)
2516 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2518 if (ir->swap->threshold < 1.0)
2520 warning_error(wi, "Ion count threshold must be at least 1.\n");
2528 static int search_QMstring(const char *s, int ng, const char *gn[])
2530 /* same as normal search_string, but this one searches QM strings */
2533 for (i = 0; (i < ng); i++)
2535 if (gmx_strcasecmp(s, gn[i]) == 0)
2541 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2545 } /* search_QMstring */
2547 /* We would like gn to be const as well, but C doesn't allow this */
2548 /* TODO this is utility functionality (search for the index of a
2549 string in a collection), so should be refactored and located more
2551 int search_string(const char *s, int ng, char *gn[])
2555 for (i = 0; (i < ng); i++)
2557 if (gmx_strcasecmp(s, gn[i]) == 0)
2564 "Group %s referenced in the .mdp file was not found in the index file.\n"
2565 "Group names must match either [moleculetype] names or custom index group\n"
2566 "names, in which case you must supply an index file to the '-n' option\n"
2573 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2574 t_blocka *block, char *gnames[],
2575 int gtype, int restnm,
2576 int grptp, gmx_bool bVerbose,
2579 unsigned short *cbuf;
2580 t_grps *grps = &(groups->grps[gtype]);
2581 int i, j, gid, aj, ognr, ntot = 0;
2584 char warn_buf[STRLEN];
2588 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2591 title = gtypes[gtype];
2594 /* Mark all id's as not set */
2595 for (i = 0; (i < natoms); i++)
2600 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2601 for (i = 0; (i < ng); i++)
2603 /* Lookup the group name in the block structure */
2604 gid = search_string(ptrs[i], block->nr, gnames);
2605 if ((grptp != egrptpONE) || (i == 0))
2607 grps->nm_ind[grps->nr++] = gid;
2611 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2614 /* Now go over the atoms in the group */
2615 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2620 /* Range checking */
2621 if ((aj < 0) || (aj >= natoms))
2623 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2625 /* Lookup up the old group number */
2629 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2630 aj+1, title, ognr+1, i+1);
2634 /* Store the group number in buffer */
2635 if (grptp == egrptpONE)
2648 /* Now check whether we have done all atoms */
2652 if (grptp == egrptpALL)
2654 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2655 natoms-ntot, title);
2657 else if (grptp == egrptpPART)
2659 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2660 natoms-ntot, title);
2661 warning_note(wi, warn_buf);
2663 /* Assign all atoms currently unassigned to a rest group */
2664 for (j = 0; (j < natoms); j++)
2666 if (cbuf[j] == NOGID)
2672 if (grptp != egrptpPART)
2677 "Making dummy/rest group for %s containing %d elements\n",
2678 title, natoms-ntot);
2680 /* Add group name "rest" */
2681 grps->nm_ind[grps->nr] = restnm;
2683 /* Assign the rest name to all atoms not currently assigned to a group */
2684 for (j = 0; (j < natoms); j++)
2686 if (cbuf[j] == NOGID)
2695 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2697 /* All atoms are part of one (or no) group, no index required */
2698 groups->ngrpnr[gtype] = 0;
2699 groups->grpnr[gtype] = NULL;
2703 groups->ngrpnr[gtype] = natoms;
2704 snew(groups->grpnr[gtype], natoms);
2705 for (j = 0; (j < natoms); j++)
2707 groups->grpnr[gtype][j] = cbuf[j];
2713 return (bRest && grptp == egrptpPART);
2716 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2719 gmx_groups_t *groups;
2720 pull_params_t *pull;
2721 int natoms, ai, aj, i, j, d, g, imin, jmin;
2723 int *nrdf2, *na_vcm, na_tot;
2724 double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
2726 gmx_mtop_atomloop_all_t aloop;
2728 int mb, mol, ftype, as;
2729 gmx_molblock_t *molb;
2730 gmx_moltype_t *molt;
2733 * First calc 3xnr-atoms for each group
2734 * then subtract half a degree of freedom for each constraint
2736 * Only atoms and nuclei contribute to the degrees of freedom...
2741 groups = &mtop->groups;
2742 natoms = mtop->natoms;
2744 /* Allocate one more for a possible rest group */
2745 /* We need to sum degrees of freedom into doubles,
2746 * since floats give too low nrdf's above 3 million atoms.
2748 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2749 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2750 snew(dof_vcm, groups->grps[egcVCM].nr+1);
2751 snew(na_vcm, groups->grps[egcVCM].nr+1);
2752 snew(nrdf_vcm_sub, groups->grps[egcVCM].nr+1);
2754 for (i = 0; i < groups->grps[egcTC].nr; i++)
2758 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2761 clear_ivec(dof_vcm[i]);
2763 nrdf_vcm_sub[i] = 0;
2766 snew(nrdf2, natoms);
2767 aloop = gmx_mtop_atomloop_all_init(mtop);
2768 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2771 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2773 g = ggrpnr(groups, egcFREEZE, i);
2774 for (d = 0; d < DIM; d++)
2776 if (opts->nFreeze[g][d] == 0)
2778 /* Add one DOF for particle i (counted as 2*1) */
2780 /* VCM group i has dim d as a DOF */
2781 dof_vcm[ggrpnr(groups, egcVCM, i)][d] = 1;
2784 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2785 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2790 for (mb = 0; mb < mtop->nmolblock; mb++)
2792 molb = &mtop->molblock[mb];
2793 molt = &mtop->moltype[molb->type];
2794 atom = molt->atoms.atom;
2795 for (mol = 0; mol < molb->nmol; mol++)
2797 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2799 ia = molt->ilist[ftype].iatoms;
2800 for (i = 0; i < molt->ilist[ftype].nr; )
2802 /* Subtract degrees of freedom for the constraints,
2803 * if the particles still have degrees of freedom left.
2804 * If one of the particles is a vsite or a shell, then all
2805 * constraint motion will go there, but since they do not
2806 * contribute to the constraints the degrees of freedom do not
2811 if (((atom[ia[1]].ptype == eptNucleus) ||
2812 (atom[ia[1]].ptype == eptAtom)) &&
2813 ((atom[ia[2]].ptype == eptNucleus) ||
2814 (atom[ia[2]].ptype == eptAtom)))
2832 imin = std::min(imin, nrdf2[ai]);
2833 jmin = std::min(jmin, nrdf2[aj]);
2836 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2837 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2838 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2839 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2841 ia += interaction_function[ftype].nratoms+1;
2842 i += interaction_function[ftype].nratoms+1;
2845 ia = molt->ilist[F_SETTLE].iatoms;
2846 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2848 /* Subtract 1 dof from every atom in the SETTLE */
2849 for (j = 0; j < 3; j++)
2852 imin = std::min(2, nrdf2[ai]);
2854 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2855 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2860 as += molt->atoms.nr;
2866 /* Correct nrdf for the COM constraints.
2867 * We correct using the TC and VCM group of the first atom
2868 * in the reference and pull group. If atoms in one pull group
2869 * belong to different TC or VCM groups it is anyhow difficult
2870 * to determine the optimal nrdf assignment.
2874 for (i = 0; i < pull->ncoord; i++)
2876 if (pull->coord[i].eType != epullCONSTRAINT)
2883 for (j = 0; j < 2; j++)
2885 const t_pull_group *pgrp;
2887 pgrp = &pull->group[pull->coord[i].group[j]];
2891 /* Subtract 1/2 dof from each group */
2893 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2894 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2895 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2897 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)]]);
2902 /* We need to subtract the whole DOF from group j=1 */
2909 if (ir->nstcomm != 0)
2913 /* We remove COM motion up to dim ndof_com() */
2914 ndim_rm_vcm = ndof_com(ir);
2916 /* Subtract ndim_rm_vcm (or less with frozen dimensions) from
2917 * the number of degrees of freedom in each vcm group when COM
2918 * translation is removed and 6 when rotation is removed as well.
2920 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2922 switch (ir->comm_mode)
2925 nrdf_vcm_sub[j] = 0;
2926 for (d = 0; d < ndim_rm_vcm; d++)
2935 nrdf_vcm_sub[j] = 6;
2938 gmx_incons("Checking comm_mode");
2942 for (i = 0; i < groups->grps[egcTC].nr; i++)
2944 /* Count the number of atoms of TC group i for every VCM group */
2945 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2950 for (ai = 0; ai < natoms; ai++)
2952 if (ggrpnr(groups, egcTC, ai) == i)
2954 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2958 /* Correct for VCM removal according to the fraction of each VCM
2959 * group present in this TC group.
2961 nrdf_uc = nrdf_tc[i];
2964 fprintf(debug, "T-group[%d] nrdf_uc = %g\n", i, nrdf_uc);
2967 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2969 if (nrdf_vcm[j] > nrdf_vcm_sub[j])
2971 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2972 (nrdf_vcm[j] - nrdf_vcm_sub[j])/nrdf_vcm[j];
2976 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2977 j, nrdf_vcm[j], nrdf_tc[i]);
2982 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2984 opts->nrdf[i] = nrdf_tc[i];
2985 if (opts->nrdf[i] < 0)
2990 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2991 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2999 sfree(nrdf_vcm_sub);
3002 static void decode_cos(char *s, t_cosines *cosine)
3005 char format[STRLEN], f1[STRLEN];
3017 sscanf(t, "%d", &(cosine->n));
3024 snew(cosine->a, cosine->n);
3025 snew(cosine->phi, cosine->n);
3027 sprintf(format, "%%*d");
3028 for (i = 0; (i < cosine->n); i++)
3031 strcat(f1, "%lf%lf");
3032 if (sscanf(t, f1, &a, &phi) < 2)
3034 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
3037 cosine->phi[i] = phi;
3038 strcat(format, "%*lf%*lf");
3045 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
3046 const char *option, const char *val, int flag)
3048 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
3049 * But since this is much larger than STRLEN, such a line can not be parsed.
3050 * The real maximum is the number of names that fit in a string: STRLEN/2.
3052 #define EGP_MAX (STRLEN/2)
3053 int nelem, i, j, k, nr;
3054 char *names[EGP_MAX];
3058 gnames = groups->grpname;
3060 nelem = str_nelem(val, EGP_MAX, names);
3063 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
3065 nr = groups->grps[egcENER].nr;
3067 for (i = 0; i < nelem/2; i++)
3071 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
3077 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3078 names[2*i], option);
3082 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3088 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3089 names[2*i+1], option);
3091 if ((j < nr) && (k < nr))
3093 ir->opts.egp_flags[nr*j+k] |= flag;
3094 ir->opts.egp_flags[nr*k+j] |= flag;
3103 static void make_swap_groups(
3108 int ig = -1, i = 0, gind;
3112 /* Just a quick check here, more thorough checks are in mdrun */
3113 if (strcmp(swap->grp[eGrpSplit0].molname, swap->grp[eGrpSplit1].molname) == 0)
3115 gmx_fatal(FARGS, "The split groups can not both be '%s'.", swap->grp[eGrpSplit0].molname);
3118 /* Get the index atoms of the split0, split1, solvent, and swap groups */
3119 for (ig = 0; ig < swap->ngrp; ig++)
3121 swapg = &swap->grp[ig];
3122 gind = search_string(swap->grp[ig].molname, grps->nr, gnames);
3123 swapg->nat = grps->index[gind+1] - grps->index[gind];
3127 fprintf(stderr, "%s group '%s' contains %d atoms.\n",
3128 ig < 3 ? eSwapFixedGrp_names[ig] : "Swap",
3129 swap->grp[ig].molname, swapg->nat);
3130 snew(swapg->ind, swapg->nat);
3131 for (i = 0; i < swapg->nat; i++)
3133 swapg->ind[i] = grps->a[grps->index[gind]+i];
3138 gmx_fatal(FARGS, "Swap group %s does not contain any atoms.", swap->grp[ig].molname);
3144 void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3149 ig = search_string(IMDgname, grps->nr, gnames);
3150 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3152 if (IMDgroup->nat > 0)
3154 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3155 IMDgname, IMDgroup->nat);
3156 snew(IMDgroup->ind, IMDgroup->nat);
3157 for (i = 0; i < IMDgroup->nat; i++)
3159 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3165 void do_index(const char* mdparin, const char *ndx,
3172 gmx_groups_t *groups;
3176 char warnbuf[STRLEN], **gnames;
3177 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3180 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3181 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3182 int i, j, k, restnm;
3183 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3184 int nQMmethod, nQMbasis, nQMg;
3185 char warn_buf[STRLEN];
3190 fprintf(stderr, "processing index file...\n");
3195 snew(grps->index, 1);
3197 atoms_all = gmx_mtop_global_atoms(mtop);
3198 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3199 done_atom(&atoms_all);
3203 grps = init_index(ndx, &gnames);
3206 groups = &mtop->groups;
3207 natoms = mtop->natoms;
3208 symtab = &mtop->symtab;
3210 snew(groups->grpname, grps->nr+1);
3212 for (i = 0; (i < grps->nr); i++)
3214 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3216 groups->grpname[i] = put_symtab(symtab, "rest");
3218 srenew(gnames, grps->nr+1);
3219 gnames[restnm] = *(groups->grpname[i]);
3220 groups->ngrpname = grps->nr+1;
3222 set_warning_line(wi, mdparin, -1);
3224 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3225 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3226 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3227 if ((ntau_t != ntcg) || (nref_t != ntcg))
3229 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3230 "%d tau-t values", ntcg, nref_t, ntau_t);
3233 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3234 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3235 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3236 nr = groups->grps[egcTC].nr;
3238 snew(ir->opts.nrdf, nr);
3239 snew(ir->opts.tau_t, nr);
3240 snew(ir->opts.ref_t, nr);
3241 if (ir->eI == eiBD && ir->bd_fric == 0)
3243 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3250 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3254 for (i = 0; (i < nr); i++)
3256 ir->opts.tau_t[i] = strtod(ptr1[i], &endptr);
3259 warning_error(wi, "Invalid value for mdp option tau-t. tau-t should only consist of real numbers separated by spaces.");
3261 if ((ir->eI == eiBD) && ir->opts.tau_t[i] <= 0)
3263 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3264 warning_error(wi, warn_buf);
3267 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3269 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.");
3272 if (ir->opts.tau_t[i] >= 0)
3274 tau_min = std::min(tau_min, ir->opts.tau_t[i]);
3277 if (ir->etc != etcNO && ir->nsttcouple == -1)
3279 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3284 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3286 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");
3288 if (ir->epc == epcMTTK)
3290 if (ir->etc != etcNOSEHOOVER)
3292 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
3296 if (ir->nstpcouple != ir->nsttcouple)
3298 int mincouple = std::min(ir->nstpcouple, ir->nsttcouple);
3299 ir->nstpcouple = ir->nsttcouple = mincouple;
3300 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);
3301 warning_note(wi, warn_buf);
3306 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3307 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3309 if (ETC_ANDERSEN(ir->etc))
3311 if (ir->nsttcouple != 1)
3314 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");
3315 warning_note(wi, warn_buf);
3318 nstcmin = tcouple_min_integration_steps(ir->etc);
3321 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3323 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3324 ETCOUPLTYPE(ir->etc),
3326 ir->nsttcouple*ir->delta_t);
3327 warning(wi, warn_buf);
3330 for (i = 0; (i < nr); i++)
3332 ir->opts.ref_t[i] = strtod(ptr2[i], &endptr);
3335 warning_error(wi, "Invalid value for mdp option ref-t. ref-t should only consist of real numbers separated by spaces.");
3337 if (ir->opts.ref_t[i] < 0)
3339 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3342 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3343 if we are in this conditional) if mc_temp is negative */
3344 if (ir->expandedvals->mc_temp < 0)
3346 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3350 /* Simulated annealing for each group. There are nr groups */
3351 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3352 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3356 if (nSA > 0 && nSA != nr)
3358 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3362 snew(ir->opts.annealing, nr);
3363 snew(ir->opts.anneal_npoints, nr);
3364 snew(ir->opts.anneal_time, nr);
3365 snew(ir->opts.anneal_temp, nr);
3366 for (i = 0; i < nr; i++)
3368 ir->opts.annealing[i] = eannNO;
3369 ir->opts.anneal_npoints[i] = 0;
3370 ir->opts.anneal_time[i] = NULL;
3371 ir->opts.anneal_temp[i] = NULL;
3376 for (i = 0; i < nr; i++)
3378 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3380 ir->opts.annealing[i] = eannNO;
3382 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3384 ir->opts.annealing[i] = eannSINGLE;
3387 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3389 ir->opts.annealing[i] = eannPERIODIC;
3395 /* Read the other fields too */
3396 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3397 if (nSA_points != nSA)
3399 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3401 for (k = 0, i = 0; i < nr; i++)
3403 ir->opts.anneal_npoints[i] = strtol(ptr1[i], &endptr, 10);
3406 warning_error(wi, "Invalid value for mdp option annealing-npoints. annealing should only consist of integers separated by spaces.");
3408 if (ir->opts.anneal_npoints[i] == 1)
3410 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3412 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3413 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3414 k += ir->opts.anneal_npoints[i];
3417 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3420 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3422 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3425 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3428 for (i = 0, k = 0; i < nr; i++)
3431 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3433 ir->opts.anneal_time[i][j] = strtod(ptr1[k], &endptr);
3436 warning_error(wi, "Invalid value for mdp option anneal-time. anneal-time should only consist of real numbers separated by spaces.");
3438 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], &endptr);
3441 warning_error(wi, "Invalid value for anneal-temp. anneal-temp should only consist of real numbers separated by spaces.");
3445 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3447 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3453 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3455 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3456 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3459 if (ir->opts.anneal_temp[i][j] < 0)
3461 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3466 /* Print out some summary information, to make sure we got it right */
3467 for (i = 0, k = 0; i < nr; i++)
3469 if (ir->opts.annealing[i] != eannNO)
3471 j = groups->grps[egcTC].nm_ind[i];
3472 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3473 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3474 ir->opts.anneal_npoints[i]);
3475 fprintf(stderr, "Time (ps) Temperature (K)\n");
3476 /* All terms except the last one */
3477 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3479 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3482 /* Finally the last one */
3483 j = ir->opts.anneal_npoints[i]-1;
3484 if (ir->opts.annealing[i] == eannSINGLE)
3486 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3490 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3491 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3493 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3504 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3506 make_pull_coords(ir->pull);
3511 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3514 if (ir->eSwapCoords != eswapNO)
3516 make_swap_groups(ir->swap, grps, gnames);
3519 /* Make indices for IMD session */
3522 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3525 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3526 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3527 if (nacg*DIM != nacc)
3529 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3532 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3533 restnm, egrptpALL_GENREST, bVerbose, wi);
3534 nr = groups->grps[egcACC].nr;
3535 snew(ir->opts.acc, nr);
3536 ir->opts.ngacc = nr;
3538 for (i = k = 0; (i < nacg); i++)
3540 for (j = 0; (j < DIM); j++, k++)
3542 ir->opts.acc[i][j] = strtod(ptr1[k], &endptr);
3545 warning_error(wi, "Invalid value for mdp option accelerate. accelerate should only consist of real numbers separated by spaces.");
3549 for (; (i < nr); i++)
3551 for (j = 0; (j < DIM); j++)
3553 ir->opts.acc[i][j] = 0;
3557 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3558 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3559 if (nfrdim != DIM*nfreeze)
3561 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3564 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3565 restnm, egrptpALL_GENREST, bVerbose, wi);
3566 nr = groups->grps[egcFREEZE].nr;
3567 ir->opts.ngfrz = nr;
3568 snew(ir->opts.nFreeze, nr);
3569 for (i = k = 0; (i < nfreeze); i++)
3571 for (j = 0; (j < DIM); j++, k++)
3573 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3574 if (!ir->opts.nFreeze[i][j])
3576 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3578 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3579 "(not %s)", ptr1[k]);
3580 warning(wi, warn_buf);
3585 for (; (i < nr); i++)
3587 for (j = 0; (j < DIM); j++)
3589 ir->opts.nFreeze[i][j] = 0;
3593 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3594 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3595 restnm, egrptpALL_GENREST, bVerbose, wi);
3596 add_wall_energrps(groups, ir->nwall, symtab);
3597 ir->opts.ngener = groups->grps[egcENER].nr;
3598 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3600 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3601 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3604 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3605 "This may lead to artifacts.\n"
3606 "In most cases one should use one group for the whole system.");
3609 /* Now we have filled the freeze struct, so we can calculate NRDF */
3610 calc_nrdf(mtop, ir, gnames);
3612 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3613 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3614 restnm, egrptpALL_GENREST, bVerbose, wi);
3615 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3616 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3617 restnm, egrptpALL_GENREST, bVerbose, wi);
3618 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3619 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3620 restnm, egrptpONE, bVerbose, wi);
3621 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3622 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3623 restnm, egrptpALL_GENREST, bVerbose, wi);
3625 /* QMMM input processing */
3626 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3627 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3628 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3629 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3631 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3632 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3634 /* group rest, if any, is always MM! */
3635 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3636 restnm, egrptpALL_GENREST, bVerbose, wi);
3637 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3638 ir->opts.ngQM = nQMg;
3639 snew(ir->opts.QMmethod, nr);
3640 snew(ir->opts.QMbasis, nr);
3641 for (i = 0; i < nr; i++)
3643 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3644 * converted to the corresponding enum in names.c
3646 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3648 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3652 str_nelem(is->QMmult, MAXPTR, ptr1);
3653 str_nelem(is->QMcharge, MAXPTR, ptr2);
3654 str_nelem(is->bSH, MAXPTR, ptr3);
3655 snew(ir->opts.QMmult, nr);
3656 snew(ir->opts.QMcharge, nr);
3657 snew(ir->opts.bSH, nr);
3659 for (i = 0; i < nr; i++)
3661 ir->opts.QMmult[i] = strtol(ptr1[i], &endptr, 10);
3664 warning_error(wi, "Invalid value for mdp option QMmult. QMmult should only consist of integers separated by spaces.");
3666 ir->opts.QMcharge[i] = strtol(ptr2[i], &endptr, 10);
3669 warning_error(wi, "Invalid value for mdp option QMcharge. QMcharge should only consist of integers separated by spaces.");
3671 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3674 str_nelem(is->CASelectrons, MAXPTR, ptr1);
3675 str_nelem(is->CASorbitals, MAXPTR, ptr2);
3676 snew(ir->opts.CASelectrons, nr);
3677 snew(ir->opts.CASorbitals, nr);
3678 for (i = 0; i < nr; i++)
3680 ir->opts.CASelectrons[i] = strtol(ptr1[i], &endptr, 10);
3683 warning_error(wi, "Invalid value for mdp option CASelectrons. CASelectrons should only consist of integers separated by spaces.");
3685 ir->opts.CASorbitals[i] = strtol(ptr2[i], &endptr, 10);
3688 warning_error(wi, "Invalid value for mdp option CASorbitals. CASorbitals should only consist of integers separated by spaces.");
3691 /* special optimization options */
3693 str_nelem(is->bOPT, MAXPTR, ptr1);
3694 str_nelem(is->bTS, MAXPTR, ptr2);
3695 snew(ir->opts.bOPT, nr);
3696 snew(ir->opts.bTS, nr);
3697 for (i = 0; i < nr; i++)
3699 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3700 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3702 str_nelem(is->SAon, MAXPTR, ptr1);
3703 str_nelem(is->SAoff, MAXPTR, ptr2);
3704 str_nelem(is->SAsteps, MAXPTR, ptr3);
3705 snew(ir->opts.SAon, nr);
3706 snew(ir->opts.SAoff, nr);
3707 snew(ir->opts.SAsteps, nr);
3709 for (i = 0; i < nr; i++)
3711 ir->opts.SAon[i] = strtod(ptr1[i], &endptr);
3714 warning_error(wi, "Invalid value for mdp option SAon. SAon should only consist of real numbers separated by spaces.");
3716 ir->opts.SAoff[i] = strtod(ptr2[i], &endptr);
3719 warning_error(wi, "Invalid value for mdp option SAoff. SAoff should only consist of real numbers separated by spaces.");
3721 ir->opts.SAsteps[i] = strtol(ptr3[i], &endptr, 10);
3724 warning_error(wi, "Invalid value for mdp option SAsteps. SAsteps should only consist of integers separated by spaces.");
3727 /* end of QMMM input */
3731 for (i = 0; (i < egcNR); i++)
3733 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3734 for (j = 0; (j < groups->grps[i].nr); j++)
3736 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3738 fprintf(stderr, "\n");
3742 nr = groups->grps[egcENER].nr;
3743 snew(ir->opts.egp_flags, nr*nr);
3745 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3746 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3748 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3750 if (bExcl && EEL_FULL(ir->coulombtype))
3752 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3755 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3756 if (bTable && !(ir->vdwtype == evdwUSER) &&
3757 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3758 !(ir->coulombtype == eelPMEUSERSWITCH))
3760 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3763 decode_cos(is->efield_x, &(ir->ex[XX]));
3764 decode_cos(is->efield_xt, &(ir->et[XX]));
3765 decode_cos(is->efield_y, &(ir->ex[YY]));
3766 decode_cos(is->efield_yt, &(ir->et[YY]));
3767 decode_cos(is->efield_z, &(ir->ex[ZZ]));
3768 decode_cos(is->efield_zt, &(ir->et[ZZ]));
3770 for (i = 0; (i < grps->nr); i++)
3782 static void check_disre(gmx_mtop_t *mtop)
3784 gmx_ffparams_t *ffparams;
3785 t_functype *functype;
3787 int i, ndouble, ftype;
3788 int label, old_label;
3790 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3792 ffparams = &mtop->ffparams;
3793 functype = ffparams->functype;
3794 ip = ffparams->iparams;
3797 for (i = 0; i < ffparams->ntypes; i++)
3799 ftype = functype[i];
3800 if (ftype == F_DISRES)
3802 label = ip[i].disres.label;
3803 if (label == old_label)
3805 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3813 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3814 "probably the parameters for multiple pairs in one restraint "
3815 "are not identical\n", ndouble);
3820 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3821 gmx_bool posres_only,
3825 gmx_mtop_ilistloop_t iloop;
3835 for (d = 0; d < DIM; d++)
3837 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3839 /* Check for freeze groups */
3840 for (g = 0; g < ir->opts.ngfrz; g++)
3842 for (d = 0; d < DIM; d++)
3844 if (ir->opts.nFreeze[g][d] != 0)
3852 /* Check for position restraints */
3853 iloop = gmx_mtop_ilistloop_init(sys);
3854 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3857 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3859 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3861 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3862 for (d = 0; d < DIM; d++)
3864 if (pr->posres.fcA[d] != 0)
3870 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3872 /* Check for flat-bottom posres */
3873 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3874 if (pr->fbposres.k != 0)
3876 switch (pr->fbposres.geom)
3878 case efbposresSPHERE:
3879 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3881 case efbposresCYLINDERX:
3882 AbsRef[YY] = AbsRef[ZZ] = 1;
3884 case efbposresCYLINDERY:
3885 AbsRef[XX] = AbsRef[ZZ] = 1;
3887 case efbposresCYLINDER:
3888 /* efbposres is a synonym for efbposresCYLINDERZ for backwards compatibility */
3889 case efbposresCYLINDERZ:
3890 AbsRef[XX] = AbsRef[YY] = 1;
3892 case efbposresX: /* d=XX */
3893 case efbposresY: /* d=YY */
3894 case efbposresZ: /* d=ZZ */
3895 d = pr->fbposres.geom - efbposresX;
3899 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3900 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3908 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3912 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3913 gmx_bool *bC6ParametersWorkWithGeometricRules,
3914 gmx_bool *bC6ParametersWorkWithLBRules,
3915 gmx_bool *bLBRulesPossible)
3917 int ntypes, tpi, tpj;
3920 double c6i, c6j, c12i, c12j;
3921 double c6, c6_geometric, c6_LB;
3922 double sigmai, sigmaj, epsi, epsj;
3923 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3926 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3927 * force-field floating point parameters.
3930 ptr = getenv("GMX_LJCOMB_TOL");
3935 sscanf(ptr, "%lf", &dbl);
3939 *bC6ParametersWorkWithLBRules = TRUE;
3940 *bC6ParametersWorkWithGeometricRules = TRUE;
3941 bCanDoLBRules = TRUE;
3942 ntypes = mtop->ffparams.atnr;
3943 snew(typecount, ntypes);
3944 gmx_mtop_count_atomtypes(mtop, state, typecount);
3945 *bLBRulesPossible = TRUE;
3946 for (tpi = 0; tpi < ntypes; ++tpi)
3948 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3949 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3950 for (tpj = tpi; tpj < ntypes; ++tpj)
3952 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3953 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3954 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3955 c6_geometric = std::sqrt(c6i * c6j);
3956 if (!gmx_numzero(c6_geometric))
3958 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3960 sigmai = gmx::sixthroot(c12i / c6i);
3961 sigmaj = gmx::sixthroot(c12j / c6j);
3962 epsi = c6i * c6i /(4.0 * c12i);
3963 epsj = c6j * c6j /(4.0 * c12j);
3964 c6_LB = 4.0 * std::sqrt(epsi * epsj) * gmx::power6(0.5 * (sigmai + sigmaj));
3968 *bLBRulesPossible = FALSE;
3969 c6_LB = c6_geometric;
3971 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3974 if (FALSE == bCanDoLBRules)
3976 *bC6ParametersWorkWithLBRules = FALSE;
3979 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3981 if (FALSE == bCanDoGeometricRules)
3983 *bC6ParametersWorkWithGeometricRules = FALSE;
3991 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3994 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3996 check_combination_rule_differences(mtop, 0,
3997 &bC6ParametersWorkWithGeometricRules,
3998 &bC6ParametersWorkWithLBRules,
4000 if (ir->ljpme_combination_rule == eljpmeLB)
4002 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
4004 warning(wi, "You are using arithmetic-geometric combination rules "
4005 "in LJ-PME, but your non-bonded C6 parameters do not "
4006 "follow these rules.");
4011 if (FALSE == bC6ParametersWorkWithGeometricRules)
4013 if (ir->eDispCorr != edispcNO)
4015 warning_note(wi, "You are using geometric combination rules in "
4016 "LJ-PME, but your non-bonded C6 parameters do "
4017 "not follow these rules. "
4018 "This will introduce very small errors in the forces and energies in "
4019 "your simulations. Dispersion correction will correct total energy "
4020 "and/or pressure for isotropic systems, but not forces or surface tensions.");
4024 warning_note(wi, "You are using geometric combination rules in "
4025 "LJ-PME, but your non-bonded C6 parameters do "
4026 "not follow these rules. "
4027 "This will introduce very small errors in the forces and energies in "
4028 "your simulations. If your system is homogeneous, consider using dispersion correction "
4029 "for the total energy and pressure.");
4035 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
4038 char err_buf[STRLEN];
4040 gmx_bool bCharge, bAcc;
4043 gmx_mtop_atomloop_block_t aloopb;
4044 gmx_mtop_atomloop_all_t aloop;
4047 char warn_buf[STRLEN];
4049 set_warning_line(wi, mdparin, -1);
4051 if (ir->cutoff_scheme == ecutsVERLET &&
4052 ir->verletbuf_tol > 0 &&
4054 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
4055 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
4057 /* Check if a too small Verlet buffer might potentially
4058 * cause more drift than the thermostat can couple off.
4060 /* Temperature error fraction for warning and suggestion */
4061 const real T_error_warn = 0.002;
4062 const real T_error_suggest = 0.001;
4063 /* For safety: 2 DOF per atom (typical with constraints) */
4064 const real nrdf_at = 2;
4065 real T, tau, max_T_error;
4070 for (i = 0; i < ir->opts.ngtc; i++)
4072 T = std::max(T, ir->opts.ref_t[i]);
4073 tau = std::max(tau, ir->opts.tau_t[i]);
4077 /* This is a worst case estimate of the temperature error,
4078 * assuming perfect buffer estimation and no cancelation
4079 * of errors. The factor 0.5 is because energy distributes
4080 * equally over Ekin and Epot.
4082 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4083 if (max_T_error > T_error_warn)
4085 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.",
4086 ir->verletbuf_tol, T, tau,
4088 100*T_error_suggest,
4089 ir->verletbuf_tol*T_error_suggest/max_T_error);
4090 warning(wi, warn_buf);
4095 if (ETC_ANDERSEN(ir->etc))
4099 for (i = 0; i < ir->opts.ngtc; i++)
4101 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4102 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4103 sprintf(err_buf, "all tau_t must be positive using Andersen temperature control, tau_t[%d]=%10.6f",
4104 i, ir->opts.tau_t[i]);
4105 CHECK(ir->opts.tau_t[i] < 0);
4108 for (i = 0; i < ir->opts.ngtc; i++)
4110 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
4111 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);
4112 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
4116 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4117 ir->comm_mode == ecmNO &&
4118 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4119 !ETC_ANDERSEN(ir->etc))
4121 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");
4124 /* Check for pressure coupling with absolute position restraints */
4125 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4127 absolute_reference(ir, sys, TRUE, AbsRef);
4129 for (m = 0; m < DIM; m++)
4131 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4133 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4141 aloopb = gmx_mtop_atomloop_block_init(sys);
4142 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4144 if (atom->q != 0 || atom->qB != 0)
4152 if (EEL_FULL(ir->coulombtype))
4155 "You are using full electrostatics treatment %s for a system without charges.\n"
4156 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4157 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4158 warning(wi, err_buf);
4163 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4166 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4167 "You might want to consider using %s electrostatics.\n",
4169 warning_note(wi, err_buf);
4173 /* Check if combination rules used in LJ-PME are the same as in the force field */
4174 if (EVDW_PME(ir->vdwtype))
4176 check_combination_rules(ir, sys, wi);
4179 /* Generalized reaction field */
4180 if (ir->opts.ngtc == 0)
4182 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4184 CHECK(ir->coulombtype == eelGRF);
4188 sprintf(err_buf, "When using coulombtype = %s"
4189 " ref-t for temperature coupling should be > 0",
4191 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4195 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4197 for (m = 0; (m < DIM); m++)
4199 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4208 snew(mgrp, sys->groups.grps[egcACC].nr);
4209 aloop = gmx_mtop_atomloop_all_init(sys);
4210 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4212 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4215 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4217 for (m = 0; (m < DIM); m++)
4219 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4223 for (m = 0; (m < DIM); m++)
4225 if (fabs(acc[m]) > 1e-6)
4227 const char *dim[DIM] = { "X", "Y", "Z" };
4229 "Net Acceleration in %s direction, will %s be corrected\n",
4230 dim[m], ir->nstcomm != 0 ? "" : "not");
4231 if (ir->nstcomm != 0 && m < ndof_com(ir))
4234 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4236 ir->opts.acc[i][m] -= acc[m];
4244 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4245 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4247 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4255 for (i = 0; i < ir->pull->ncoord && !bWarned; i++)
4257 if (ir->pull->coord[i].group[0] == 0 ||
4258 ir->pull->coord[i].group[1] == 0)
4260 absolute_reference(ir, sys, FALSE, AbsRef);
4261 for (m = 0; m < DIM; m++)
4263 if (ir->pull->coord[i].dim[m] && !AbsRef[m])
4265 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.");
4273 for (i = 0; i < 3; i++)
4275 for (m = 0; m <= i; m++)
4277 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4278 ir->deform[i][m] != 0)
4280 for (c = 0; c < ir->pull->ncoord; c++)
4282 if (ir->pull->coord[c].eGeom == epullgDIRPBC &&
4283 ir->pull->coord[c].vec[m] != 0)
4285 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->coord[c].eGeom), 'x'+m);
4296 void double_check(t_inputrec *ir, matrix box,
4297 gmx_bool bHasNormalConstraints,
4298 gmx_bool bHasAnyConstraints,
4302 char warn_buf[STRLEN];
4305 ptr = check_box(ir->ePBC, box);
4308 warning_error(wi, ptr);
4311 if (bHasNormalConstraints && ir->eConstrAlg == econtSHAKE)
4313 if (ir->shake_tol <= 0.0)
4315 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4317 warning_error(wi, warn_buf);
4321 if ( (ir->eConstrAlg == econtLINCS) && bHasNormalConstraints)
4323 /* If we have Lincs constraints: */
4324 if (ir->eI == eiMD && ir->etc == etcNO &&
4325 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4327 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4328 warning_note(wi, warn_buf);
4331 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4333 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4334 warning_note(wi, warn_buf);
4336 if (ir->epc == epcMTTK)
4338 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4342 if (bHasAnyConstraints && ir->epc == epcMTTK)
4344 warning_error(wi, "Constraints are not implemented with MTTK pressure control.");
4347 if (ir->LincsWarnAngle > 90.0)
4349 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4350 warning(wi, warn_buf);
4351 ir->LincsWarnAngle = 90.0;
4354 if (ir->ePBC != epbcNONE)
4356 if (ir->nstlist == 0)
4358 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4360 if (ir->ns_type == ensGRID)
4362 if (gmx::square(ir->rlist) >= max_cutoff2(ir->ePBC, box))
4364 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");
4365 warning_error(wi, warn_buf);
4370 min_size = std::min(box[XX][XX], std::min(box[YY][YY], box[ZZ][ZZ]));
4371 if (2*ir->rlist >= min_size)
4373 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4374 warning_error(wi, warn_buf);
4377 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4384 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4388 real rvdw1, rvdw2, rcoul1, rcoul2;
4389 char warn_buf[STRLEN];
4391 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4395 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4400 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4406 if (rvdw1 + rvdw2 > ir->rlist ||
4407 rcoul1 + rcoul2 > ir->rlist)
4410 "The sum of the two largest charge group radii (%f) "
4411 "is larger than rlist (%f)\n",
4412 std::max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4413 warning(wi, warn_buf);
4417 /* Here we do not use the zero at cut-off macro,
4418 * since user defined interactions might purposely
4419 * not be zero at the cut-off.
4421 if (ir_vdw_is_zero_at_cutoff(ir) &&
4422 rvdw1 + rvdw2 > ir->rlist - ir->rvdw)
4424 sprintf(warn_buf, "The sum of the two largest charge group "
4425 "radii (%f) is larger than rlist (%f) - rvdw (%f).\n"
4426 "With exact cut-offs, better performance can be "
4427 "obtained with cutoff-scheme = %s, because it "
4428 "does not use charge groups at all.",
4430 ir->rlist, ir->rvdw,
4431 ecutscheme_names[ecutsVERLET]);
4434 warning(wi, warn_buf);
4438 warning_note(wi, warn_buf);
4441 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4442 rcoul1 + rcoul2 > ir->rlist - ir->rcoulomb)
4444 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than rlist (%f) - rcoulomb (%f).\n"
4445 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4447 ir->rlist, ir->rcoulomb,
4448 ecutscheme_names[ecutsVERLET]);
4451 warning(wi, warn_buf);
4455 warning_note(wi, warn_buf);