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, 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.
43 #include "gromacs/legacyheaders/typedefs.h"
44 #include "gromacs/math/units.h"
45 #include "gromacs/legacyheaders/names.h"
46 #include "gromacs/legacyheaders/macros.h"
47 #include "gromacs/topology/index.h"
48 #include "gromacs/utility/cstringutil.h"
49 #include "gromacs/legacyheaders/readinp.h"
50 #include "gromacs/legacyheaders/warninp.h"
53 #include "gromacs/legacyheaders/network.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/topology/mtop_util.h"
57 #include "gromacs/legacyheaders/chargegroup.h"
58 #include "gromacs/legacyheaders/inputrec.h"
59 #include "calc_verletbuf.h"
61 #include "gromacs/topology/block.h"
62 #include "gromacs/topology/symtab.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/smalloc.h"
69 /* Resource parameters
70 * Do not change any of these until you read the instruction
71 * in readinp.h. Some cpp's do not take spaces after the backslash
72 * (like the c-shell), which will give you a very weird compiler
76 typedef struct t_inputrec_strings
78 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
79 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
80 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
81 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
82 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
84 char fep_lambda[efptNR][STRLEN];
85 char lambda_weights[STRLEN];
88 char anneal[STRLEN], anneal_npoints[STRLEN],
89 anneal_time[STRLEN], anneal_temp[STRLEN];
90 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
91 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
92 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
93 char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
94 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
96 } gmx_inputrec_strings;
98 static gmx_inputrec_strings *is = NULL;
100 void init_inputrec_strings()
104 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.");
109 void done_inputrec_strings()
115 static char swapgrp[STRLEN], splitgrp0[STRLEN], splitgrp1[STRLEN], solgrp[STRLEN];
118 egrptpALL, /* All particles have to be a member of a group. */
119 egrptpALL_GENREST, /* A rest group with name is generated for particles *
120 * that are not part of any group. */
121 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
122 * for the rest group. */
123 egrptpONE /* Merge all selected groups into one group, *
124 * make a rest group for the remaining particles. */
127 static const char *constraints[eshNR+1] = {
128 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
131 static const char *couple_lam[ecouplamNR+1] = {
132 "vdw-q", "vdw", "q", "none", NULL
135 void init_ir(t_inputrec *ir, t_gromppopts *opts)
137 snew(opts->include, STRLEN);
138 snew(opts->define, STRLEN);
139 snew(ir->fepvals, 1);
140 snew(ir->expandedvals, 1);
141 snew(ir->simtempvals, 1);
144 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
149 for (i = 0; i < ntemps; i++)
151 /* simple linear scaling -- allows more control */
152 if (simtemp->eSimTempScale == esimtempLINEAR)
154 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
156 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
158 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
160 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
162 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
167 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
168 gmx_fatal(FARGS, errorstr);
175 static void _low_check(gmx_bool b, char *s, warninp_t wi)
179 warning_error(wi, s);
183 static void check_nst(const char *desc_nst, int nst,
184 const char *desc_p, int *p,
189 if (*p > 0 && *p % nst != 0)
191 /* Round up to the next multiple of nst */
192 *p = ((*p)/nst + 1)*nst;
193 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
194 desc_p, desc_nst, desc_p, *p);
199 static gmx_bool ir_NVE(const t_inputrec *ir)
201 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
204 static int lcd(int n1, int n2)
209 for (i = 2; (i <= n1 && i <= n2); i++)
211 if (n1 % i == 0 && n2 % i == 0)
220 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
222 if (*eintmod == eintmodPOTSHIFT_VERLET)
224 if (ir->cutoff_scheme == ecutsVERLET)
226 *eintmod = eintmodPOTSHIFT;
230 *eintmod = eintmodNONE;
235 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
237 /* Check internal consistency.
238 * NOTE: index groups are not set here yet, don't check things
239 * like temperature coupling group options here, but in triple_check
242 /* Strange macro: first one fills the err_buf, and then one can check
243 * the condition, which will print the message and increase the error
246 #define CHECK(b) _low_check(b, err_buf, wi)
247 char err_buf[256], warn_buf[STRLEN];
253 t_lambda *fep = ir->fepvals;
254 t_expanded *expand = ir->expandedvals;
256 set_warning_line(wi, mdparin, -1);
258 /* BASIC CUT-OFF STUFF */
259 if (ir->rcoulomb < 0)
261 warning_error(wi, "rcoulomb should be >= 0");
265 warning_error(wi, "rvdw should be >= 0");
268 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
270 warning_error(wi, "rlist should be >= 0");
273 process_interaction_modifier(ir, &ir->coulomb_modifier);
274 process_interaction_modifier(ir, &ir->vdw_modifier);
276 if (ir->cutoff_scheme == ecutsGROUP)
279 "The group cutoff scheme is deprecated in Gromacs 5.0 and will be removed in a future "
280 "release when all interaction forms are supported for the verlet scheme. The verlet "
281 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
283 /* BASIC CUT-OFF STUFF */
284 if (ir->rlist == 0 ||
285 !((ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > ir->rlist) ||
286 (ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > ir->rlist)))
288 /* No switched potential and/or no twin-range:
289 * we can set the long-range cut-off to the maximum of the other cut-offs.
291 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
293 else if (ir->rlistlong < 0)
295 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
296 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
298 warning(wi, warn_buf);
300 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
302 warning_error(wi, "Can not have an infinite cut-off with PBC");
304 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
306 warning_error(wi, "rlistlong can not be shorter than rlist");
308 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
310 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
314 if (ir->rlistlong == ir->rlist)
318 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
320 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
323 if (ir->cutoff_scheme == ecutsVERLET)
327 /* Normal Verlet type neighbor-list, currently only limited feature support */
328 if (inputrec2nboundeddim(ir) < 3)
330 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
332 if (ir->rcoulomb != ir->rvdw)
334 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
336 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
338 if (ir->vdw_modifier == eintmodNONE ||
339 ir->vdw_modifier == eintmodPOTSHIFT)
341 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
343 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]);
344 warning_note(wi, warn_buf);
346 ir->vdwtype = evdwCUT;
350 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
351 warning_error(wi, warn_buf);
355 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
357 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
359 if (!(ir->coulombtype == eelCUT ||
360 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
361 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
363 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
365 if (!(ir->coulomb_modifier == eintmodNONE ||
366 ir->coulomb_modifier == eintmodPOTSHIFT))
368 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
369 warning_error(wi, warn_buf);
372 if (ir->implicit_solvent != eisNO)
374 warning_error(wi, "Implicit solvent is not (yet) supported with the with Verlet lists.");
377 if (ir->nstlist <= 0)
379 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
382 if (ir->nstlist < 10)
384 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.");
387 rc_max = max(ir->rvdw, ir->rcoulomb);
389 if (ir->verletbuf_tol <= 0)
391 if (ir->verletbuf_tol == 0)
393 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
396 if (ir->rlist < rc_max)
398 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
401 if (ir->rlist == rc_max && ir->nstlist > 1)
403 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.");
408 if (ir->rlist > rc_max)
410 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.");
413 if (ir->nstlist == 1)
415 /* No buffer required */
420 if (EI_DYNAMICS(ir->eI))
422 if (inputrec2nboundeddim(ir) < 3)
424 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.");
426 /* Set rlist temporarily so we can continue processing */
431 /* Set the buffer to 5% of the cut-off */
432 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
437 /* No twin-range calculations with Verlet lists */
438 ir->rlistlong = ir->rlist;
441 if (ir->nstcalclr == -1)
443 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
444 ir->nstcalclr = ir->nstlist;
446 else if (ir->nstcalclr > 0)
448 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
450 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
453 else if (ir->nstcalclr < -1)
455 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
458 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
460 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
463 /* GENERAL INTEGRATOR STUFF */
464 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
468 if (ir->eI == eiVVAK)
470 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]);
471 warning_note(wi, warn_buf);
473 if (!EI_DYNAMICS(ir->eI))
477 if (EI_DYNAMICS(ir->eI))
479 if (ir->nstcalcenergy < 0)
481 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
482 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
484 /* nstcalcenergy larger than nstener does not make sense.
485 * We ideally want nstcalcenergy=nstener.
489 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
493 ir->nstcalcenergy = ir->nstenergy;
497 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
498 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
499 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
502 const char *nsten = "nstenergy";
503 const char *nstdh = "nstdhdl";
504 const char *min_name = nsten;
505 int min_nst = ir->nstenergy;
507 /* find the smallest of ( nstenergy, nstdhdl ) */
508 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
509 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
511 min_nst = ir->fepvals->nstdhdl;
514 /* If the user sets nstenergy small, we should respect that */
516 "Setting nstcalcenergy (%d) equal to %s (%d)",
517 ir->nstcalcenergy, min_name, min_nst);
518 warning_note(wi, warn_buf);
519 ir->nstcalcenergy = min_nst;
522 if (ir->epc != epcNO)
524 if (ir->nstpcouple < 0)
526 ir->nstpcouple = ir_optimal_nstpcouple(ir);
529 if (IR_TWINRANGE(*ir))
531 check_nst("nstlist", ir->nstlist,
532 "nstcalcenergy", &ir->nstcalcenergy, wi);
533 if (ir->epc != epcNO)
535 check_nst("nstlist", ir->nstlist,
536 "nstpcouple", &ir->nstpcouple, wi);
540 if (ir->nstcalcenergy > 0)
542 if (ir->efep != efepNO)
544 /* nstdhdl should be a multiple of nstcalcenergy */
545 check_nst("nstcalcenergy", ir->nstcalcenergy,
546 "nstdhdl", &ir->fepvals->nstdhdl, wi);
547 /* nstexpanded should be a multiple of nstcalcenergy */
548 check_nst("nstcalcenergy", ir->nstcalcenergy,
549 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
551 /* for storing exact averages nstenergy should be
552 * a multiple of nstcalcenergy
554 check_nst("nstcalcenergy", ir->nstcalcenergy,
555 "nstenergy", &ir->nstenergy, wi);
559 if (ir->nsteps == 0 && !ir->bContinuation)
561 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
565 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
566 ir->bContinuation && ir->ld_seed != -1)
568 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)");
574 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
575 CHECK(ir->ePBC != epbcXYZ);
576 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
577 CHECK(ir->ns_type != ensGRID);
578 sprintf(err_buf, "with TPI nstlist should be larger than zero");
579 CHECK(ir->nstlist <= 0);
580 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
581 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
582 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
583 CHECK(ir->cutoff_scheme == ecutsVERLET);
587 if ( (opts->nshake > 0) && (opts->bMorse) )
590 "Using morse bond-potentials while constraining bonds is useless");
591 warning(wi, warn_buf);
594 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
595 ir->bContinuation && ir->ld_seed != -1)
597 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)");
599 /* verify simulated tempering options */
603 gmx_bool bAllTempZero = TRUE;
604 for (i = 0; i < fep->n_lambda; i++)
606 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]);
607 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
608 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
610 bAllTempZero = FALSE;
613 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
614 CHECK(bAllTempZero == TRUE);
616 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
617 CHECK(ir->eI != eiVV);
619 /* check compatability of the temperature coupling with simulated tempering */
621 if (ir->etc == etcNOSEHOOVER)
623 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
624 warning_note(wi, warn_buf);
627 /* check that the temperatures make sense */
629 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);
630 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
632 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
633 CHECK(ir->simtempvals->simtemp_high <= 0);
635 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
636 CHECK(ir->simtempvals->simtemp_low <= 0);
639 /* verify free energy options */
641 if (ir->efep != efepNO)
644 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
646 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
648 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
649 (int)fep->sc_r_power);
650 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
652 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
653 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
655 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
656 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
658 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
659 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
661 sprintf(err_buf, "Free-energy not implemented for Ewald");
662 CHECK(ir->coulombtype == eelEWALD);
664 /* check validty of lambda inputs */
665 if (fep->n_lambda == 0)
667 /* Clear output in case of no states:*/
668 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
669 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
673 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
674 CHECK((fep->init_fep_state >= fep->n_lambda));
677 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
678 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
680 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",
681 fep->init_lambda, fep->init_fep_state);
682 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
686 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
690 for (i = 0; i < efptNR; i++)
692 if (fep->separate_dvdl[i])
697 if (n_lambda_terms > 1)
699 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.");
700 warning(wi, warn_buf);
703 if (n_lambda_terms < 2 && fep->n_lambda > 0)
706 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
710 for (j = 0; j < efptNR; j++)
712 for (i = 0; i < fep->n_lambda; i++)
714 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]);
715 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
719 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
721 for (i = 0; i < fep->n_lambda; i++)
723 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],
724 fep->all_lambda[efptCOUL][i]);
725 CHECK((fep->sc_alpha > 0) &&
726 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
727 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
728 ((fep->all_lambda[efptVDW][i] > 0.0) &&
729 (fep->all_lambda[efptVDW][i] < 1.0))));
733 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
735 real sigma, lambda, r_sc;
738 /* Maximum estimate for A and B charges equal with lambda power 1 */
740 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
741 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.",
743 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
744 warning_note(wi, warn_buf);
747 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
748 be treated differently, but that's the next step */
750 for (i = 0; i < efptNR; i++)
752 for (j = 0; j < fep->n_lambda; j++)
754 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
755 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
760 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
763 expand = ir->expandedvals;
765 /* checking equilibration of weights inputs for validity */
767 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
768 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
769 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
771 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
772 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
773 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
775 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
776 expand->equil_steps, elmceq_names[elmceqSTEPS]);
777 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
779 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
780 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
781 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
783 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
784 expand->equil_ratio, elmceq_names[elmceqRATIO]);
785 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
787 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
788 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
789 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
791 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
792 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
793 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
795 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
796 expand->equil_steps, elmceq_names[elmceqSTEPS]);
797 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
799 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
800 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
801 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
803 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
804 expand->equil_ratio, elmceq_names[elmceqRATIO]);
805 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
807 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
808 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
809 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
811 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
812 CHECK((expand->lmc_repeats <= 0));
813 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
814 CHECK((expand->minvarmin <= 0));
815 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
816 CHECK((expand->c_range < 0));
817 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
818 fep->init_fep_state, expand->lmc_forced_nstart);
819 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
820 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
821 CHECK((expand->lmc_forced_nstart < 0));
822 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
823 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
825 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
826 CHECK((expand->init_wl_delta < 0));
827 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
828 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
829 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
830 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
832 /* if there is no temperature control, we need to specify an MC temperature */
833 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
834 if (expand->nstTij > 0)
836 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
837 expand->nstTij, ir->nstlog);
838 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
843 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
844 CHECK(ir->nwall && ir->ePBC != epbcXY);
847 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
849 if (ir->ePBC == epbcNONE)
851 if (ir->epc != epcNO)
853 warning(wi, "Turning off pressure coupling for vacuum system");
859 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
860 epbc_names[ir->ePBC]);
861 CHECK(ir->epc != epcNO);
863 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
864 CHECK(EEL_FULL(ir->coulombtype));
866 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
867 epbc_names[ir->ePBC]);
868 CHECK(ir->eDispCorr != edispcNO);
871 if (ir->rlist == 0.0)
873 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
874 "with coulombtype = %s or coulombtype = %s\n"
875 "without periodic boundary conditions (pbc = %s) and\n"
876 "rcoulomb and rvdw set to zero",
877 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
878 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
879 (ir->ePBC != epbcNONE) ||
880 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
884 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
888 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
893 if (ir->nstcomm == 0)
895 ir->comm_mode = ecmNO;
897 if (ir->comm_mode != ecmNO)
901 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");
902 ir->nstcomm = abs(ir->nstcomm);
905 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
907 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
908 ir->nstcomm = ir->nstcalcenergy;
911 if (ir->comm_mode == ecmANGULAR)
913 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
914 CHECK(ir->bPeriodicMols);
915 if (ir->ePBC != epbcNONE)
917 warning(wi, "Removing the rotation around the center of mass in a periodic system (this is not a problem when you have only one molecule).");
922 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
924 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.");
927 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
928 " algorithm not implemented");
929 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
930 && (ir->ns_type == ensSIMPLE));
932 /* TEMPERATURE COUPLING */
933 if (ir->etc == etcYES)
935 ir->etc = etcBERENDSEN;
936 warning_note(wi, "Old option for temperature coupling given: "
937 "changing \"yes\" to \"Berendsen\"\n");
940 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
942 if (ir->opts.nhchainlength < 1)
944 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
945 ir->opts.nhchainlength = 1;
946 warning(wi, warn_buf);
949 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
951 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
952 ir->opts.nhchainlength = 1;
957 ir->opts.nhchainlength = 0;
960 if (ir->eI == eiVVAK)
962 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
964 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
967 if (ETC_ANDERSEN(ir->etc))
969 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
970 CHECK(!(EI_VV(ir->eI)));
972 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
974 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]);
975 warning_note(wi, warn_buf);
978 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]);
979 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
982 if (ir->etc == etcBERENDSEN)
984 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
985 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
986 warning_note(wi, warn_buf);
989 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
990 && ir->epc == epcBERENDSEN)
992 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
993 "true ensemble for the thermostat");
994 warning(wi, warn_buf);
997 /* PRESSURE COUPLING */
998 if (ir->epc == epcISOTROPIC)
1000 ir->epc = epcBERENDSEN;
1001 warning_note(wi, "Old option for pressure coupling given: "
1002 "changing \"Isotropic\" to \"Berendsen\"\n");
1005 if (ir->epc != epcNO)
1007 dt_pcoupl = ir->nstpcouple*ir->delta_t;
1009 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1010 CHECK(ir->tau_p <= 0);
1012 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
1014 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
1015 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1016 warning(wi, warn_buf);
1019 sprintf(err_buf, "compressibility must be > 0 when using pressure"
1020 " coupling %s\n", EPCOUPLTYPE(ir->epc));
1021 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
1022 ir->compress[ZZ][ZZ] < 0 ||
1023 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
1024 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1026 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1029 "You are generating velocities so I am assuming you "
1030 "are equilibrating a system. You are using "
1031 "%s pressure coupling, but this can be "
1032 "unstable for equilibration. If your system crashes, try "
1033 "equilibrating first with Berendsen pressure coupling. If "
1034 "you are not equilibrating the system, you can probably "
1035 "ignore this warning.",
1036 epcoupl_names[ir->epc]);
1037 warning(wi, warn_buf);
1043 if (ir->epc > epcNO)
1045 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1047 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.");
1053 if (ir->epc == epcMTTK)
1055 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1059 /* ELECTROSTATICS */
1060 /* More checks are in triple check (grompp.c) */
1062 if (ir->coulombtype == eelSWITCH)
1064 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1065 "artifacts, advice: use coulombtype = %s",
1066 eel_names[ir->coulombtype],
1067 eel_names[eelRF_ZERO]);
1068 warning(wi, warn_buf);
1071 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1073 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1074 warning_note(wi, warn_buf);
1077 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1079 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);
1080 warning(wi, warn_buf);
1081 ir->epsilon_rf = ir->epsilon_r;
1082 ir->epsilon_r = 1.0;
1085 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == NULL)
1087 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1088 CHECK(ir->epsilon_r < 0);
1091 if (EEL_RF(ir->coulombtype))
1093 /* reaction field (at the cut-off) */
1095 if (ir->coulombtype == eelRF_ZERO)
1097 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1098 eel_names[ir->coulombtype]);
1099 CHECK(ir->epsilon_rf != 0);
1100 ir->epsilon_rf = 0.0;
1103 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1104 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1105 (ir->epsilon_r == 0));
1106 if (ir->epsilon_rf == ir->epsilon_r)
1108 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1109 eel_names[ir->coulombtype]);
1110 warning(wi, warn_buf);
1113 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1114 * means the interaction is zero outside rcoulomb, but it helps to
1115 * provide accurate energy conservation.
1117 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1119 if (ir_coulomb_switched(ir))
1122 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1123 eel_names[ir->coulombtype]);
1124 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1127 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1129 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1131 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1132 eel_names[ir->coulombtype]);
1133 CHECK(ir->rlist > ir->rcoulomb);
1137 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1140 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1141 CHECK( ir->coulomb_modifier != eintmodNONE);
1143 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1146 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1147 CHECK( ir->vdw_modifier != eintmodNONE);
1150 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1151 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1154 "The switch/shift interaction settings are just for compatibility; you will get better "
1155 "performance from applying potential modifiers to your interactions!\n");
1156 warning_note(wi, warn_buf);
1159 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1161 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1163 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1164 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.",
1165 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1166 warning(wi, warn_buf);
1170 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1172 if (ir->rvdw_switch == 0)
1174 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.");
1175 warning(wi, warn_buf);
1179 if (EEL_FULL(ir->coulombtype))
1181 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1182 ir->coulombtype == eelPMEUSERSWITCH)
1184 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1185 eel_names[ir->coulombtype]);
1186 CHECK(ir->rcoulomb > ir->rlist);
1188 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1190 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1193 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1194 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1195 "a potential modifier.", eel_names[ir->coulombtype]);
1196 if (ir->nstcalclr == 1)
1198 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1202 CHECK(ir->rcoulomb != ir->rlist);
1208 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1210 if (ir->pme_order < 3)
1212 warning_error(wi, "pme-order can not be smaller than 3");
1216 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1218 if (ir->ewald_geometry == eewg3D)
1220 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1221 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1222 warning(wi, warn_buf);
1224 /* This check avoids extra pbc coding for exclusion corrections */
1225 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1226 CHECK(ir->wall_ewald_zfac < 2);
1229 if (ir_vdw_switched(ir))
1231 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1232 CHECK(ir->rvdw_switch >= ir->rvdw);
1234 if (ir->rvdw_switch < 0.5*ir->rvdw)
1236 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.",
1237 ir->rvdw_switch, ir->rvdw);
1238 warning_note(wi, warn_buf);
1241 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1243 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1245 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1246 CHECK(ir->rlist > ir->rvdw);
1250 if (ir->vdwtype == evdwPME)
1252 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1254 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s a\
1256 evdw_names[ir->vdwtype],
1257 eintmod_names[eintmodPOTSHIFT],
1258 eintmod_names[eintmodNONE]);
1262 if (ir->cutoff_scheme == ecutsGROUP)
1264 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1265 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1268 warning_note(wi, "With exact cut-offs, rlist should be "
1269 "larger than rcoulomb and rvdw, so that there "
1270 "is a buffer region for particle motion "
1271 "between neighborsearch steps");
1274 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlistlong <= ir->rcoulomb)
1276 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1277 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1278 warning_note(wi, warn_buf);
1280 if (ir_vdw_switched(ir) && (ir->rlistlong <= ir->rvdw))
1282 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1283 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1284 warning_note(wi, warn_buf);
1288 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1290 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.");
1293 if (ir->nstlist == -1)
1295 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1296 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1298 sprintf(err_buf, "nstlist can not be smaller than -1");
1299 CHECK(ir->nstlist < -1);
1301 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1304 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1307 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1309 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1312 /* ENERGY CONSERVATION */
1313 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1315 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1317 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1318 evdw_names[evdwSHIFT]);
1319 warning_note(wi, warn_buf);
1321 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1323 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1324 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1325 warning_note(wi, warn_buf);
1329 if (EI_VV(ir->eI) && IR_TWINRANGE(*ir) && ir->nstlist > 1)
1331 sprintf(warn_buf, "Twin-range multiple time stepping does not work with integrator %s.", ei_names[ir->eI]);
1332 warning_error(wi, warn_buf);
1335 /* IMPLICIT SOLVENT */
1336 if (ir->coulombtype == eelGB_NOTUSED)
1338 sprintf(warn_buf, "Invalid option %s for coulombtype",
1339 eel_names[ir->coulombtype]);
1340 warning_error(wi, warn_buf);
1343 if (ir->sa_algorithm == esaSTILL)
1345 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1346 CHECK(ir->sa_algorithm == esaSTILL);
1349 if (ir->implicit_solvent == eisGBSA)
1351 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1352 CHECK(ir->rgbradii != ir->rlist);
1354 if (ir->coulombtype != eelCUT)
1356 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1357 CHECK(ir->coulombtype != eelCUT);
1359 if (ir->vdwtype != evdwCUT)
1361 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1362 CHECK(ir->vdwtype != evdwCUT);
1364 if (ir->nstgbradii < 1)
1366 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1367 warning_note(wi, warn_buf);
1370 if (ir->sa_algorithm == esaNO)
1372 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1373 warning_note(wi, warn_buf);
1375 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1377 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");
1378 warning_note(wi, warn_buf);
1380 if (ir->gb_algorithm == egbSTILL)
1382 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1386 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1389 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1391 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1392 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1399 if (ir->cutoff_scheme != ecutsGROUP)
1401 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1405 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1407 if (ir->epc != epcNO)
1409 warning_error(wi, "AdresS simulation does not support pressure coupling");
1411 if (EEL_FULL(ir->coulombtype))
1413 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1418 /* count the number of text elemets separated by whitespace in a string.
1419 str = the input string
1420 maxptr = the maximum number of allowed elements
1421 ptr = the output array of pointers to the first character of each element
1422 returns: the number of elements. */
1423 int str_nelem(const char *str, int maxptr, char *ptr[])
1428 copy0 = gmx_strdup(str);
1431 while (*copy != '\0')
1435 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1443 while ((*copy != '\0') && !isspace(*copy))
1462 /* interpret a number of doubles from a string and put them in an array,
1463 after allocating space for them.
1464 str = the input string
1465 n = the (pre-allocated) number of doubles read
1466 r = the output array of doubles. */
1467 static void parse_n_real(char *str, int *n, real **r)
1472 *n = str_nelem(str, MAXPTR, ptr);
1475 for (i = 0; i < *n; i++)
1477 (*r)[i] = strtod(ptr[i], NULL);
1481 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1484 int i, j, max_n_lambda, nweights, nfep[efptNR];
1485 t_lambda *fep = ir->fepvals;
1486 t_expanded *expand = ir->expandedvals;
1487 real **count_fep_lambdas;
1488 gmx_bool bOneLambda = TRUE;
1490 snew(count_fep_lambdas, efptNR);
1492 /* FEP input processing */
1493 /* first, identify the number of lambda values for each type.
1494 All that are nonzero must have the same number */
1496 for (i = 0; i < efptNR; i++)
1498 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1501 /* now, determine the number of components. All must be either zero, or equal. */
1504 for (i = 0; i < efptNR; i++)
1506 if (nfep[i] > max_n_lambda)
1508 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1509 must have the same number if its not zero.*/
1514 for (i = 0; i < efptNR; i++)
1518 ir->fepvals->separate_dvdl[i] = FALSE;
1520 else if (nfep[i] == max_n_lambda)
1522 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1523 respect to the temperature currently */
1525 ir->fepvals->separate_dvdl[i] = TRUE;
1530 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1531 nfep[i], efpt_names[i], max_n_lambda);
1534 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1535 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1537 /* the number of lambdas is the number we've read in, which is either zero
1538 or the same for all */
1539 fep->n_lambda = max_n_lambda;
1541 /* allocate space for the array of lambda values */
1542 snew(fep->all_lambda, efptNR);
1543 /* if init_lambda is defined, we need to set lambda */
1544 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1546 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1548 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1549 for (i = 0; i < efptNR; i++)
1551 snew(fep->all_lambda[i], fep->n_lambda);
1552 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1555 for (j = 0; j < fep->n_lambda; j++)
1557 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1559 sfree(count_fep_lambdas[i]);
1562 sfree(count_fep_lambdas);
1564 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1565 bookkeeping -- for now, init_lambda */
1567 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1569 for (i = 0; i < fep->n_lambda; i++)
1571 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1575 /* check to see if only a single component lambda is defined, and soft core is defined.
1576 In this case, turn on coulomb soft core */
1578 if (max_n_lambda == 0)
1584 for (i = 0; i < efptNR; i++)
1586 if ((nfep[i] != 0) && (i != efptFEP))
1592 if ((bOneLambda) && (fep->sc_alpha > 0))
1594 fep->bScCoul = TRUE;
1597 /* Fill in the others with the efptFEP if they are not explicitly
1598 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1599 they are all zero. */
1601 for (i = 0; i < efptNR; i++)
1603 if ((nfep[i] == 0) && (i != efptFEP))
1605 for (j = 0; j < fep->n_lambda; j++)
1607 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1613 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1614 if (fep->sc_r_power == 48)
1616 if (fep->sc_alpha > 0.1)
1618 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1622 expand = ir->expandedvals;
1623 /* now read in the weights */
1624 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1627 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1629 else if (nweights != fep->n_lambda)
1631 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1632 nweights, fep->n_lambda);
1634 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1636 expand->nstexpanded = fep->nstdhdl;
1637 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1639 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1641 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1642 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1643 2*tau_t just to be careful so it's not to frequent */
1648 static void do_simtemp_params(t_inputrec *ir)
1651 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1652 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1657 static void do_wall_params(t_inputrec *ir,
1658 char *wall_atomtype, char *wall_density,
1662 char *names[MAXPTR];
1665 opts->wall_atomtype[0] = NULL;
1666 opts->wall_atomtype[1] = NULL;
1668 ir->wall_atomtype[0] = -1;
1669 ir->wall_atomtype[1] = -1;
1670 ir->wall_density[0] = 0;
1671 ir->wall_density[1] = 0;
1675 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1676 if (nstr != ir->nwall)
1678 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1681 for (i = 0; i < ir->nwall; i++)
1683 opts->wall_atomtype[i] = gmx_strdup(names[i]);
1686 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1688 nstr = str_nelem(wall_density, MAXPTR, names);
1689 if (nstr != ir->nwall)
1691 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1693 for (i = 0; i < ir->nwall; i++)
1695 sscanf(names[i], "%lf", &dbl);
1698 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1700 ir->wall_density[i] = dbl;
1706 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1714 srenew(groups->grpname, groups->ngrpname+nwall);
1715 grps = &(groups->grps[egcENER]);
1716 srenew(grps->nm_ind, grps->nr+nwall);
1717 for (i = 0; i < nwall; i++)
1719 sprintf(str, "wall%d", i);
1720 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1721 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1726 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1727 t_expanded *expand, warninp_t wi)
1729 int ninp, nerror = 0;
1735 /* read expanded ensemble parameters */
1736 CCTYPE ("expanded ensemble variables");
1737 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1738 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1739 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1740 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1741 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1742 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1743 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1744 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1745 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1746 CCTYPE("Seed for Monte Carlo in lambda space");
1747 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1748 RTYPE ("mc-temperature", expand->mc_temp, -1);
1749 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1750 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1751 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1752 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1753 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1754 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1755 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1756 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1757 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1758 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1759 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1767 void get_ir(const char *mdparin, const char *mdparout,
1768 t_inputrec *ir, t_gromppopts *opts,
1772 double dumdub[2][6];
1776 char warn_buf[STRLEN];
1777 t_lambda *fep = ir->fepvals;
1778 t_expanded *expand = ir->expandedvals;
1780 init_inputrec_strings();
1781 inp = read_inpfile(mdparin, &ninp, wi);
1783 snew(dumstr[0], STRLEN);
1784 snew(dumstr[1], STRLEN);
1786 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1789 "%s did not specify a value for the .mdp option "
1790 "\"cutoff-scheme\". Probably it was first intended for use "
1791 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1792 "introduced, but the group scheme was still the default. "
1793 "The default is now the Verlet scheme, so you will observe "
1794 "different behaviour.", mdparin);
1795 warning_note(wi, warn_buf);
1798 /* ignore the following deprecated commands */
1801 REM_TYPE("domain-decomposition");
1802 REM_TYPE("andersen-seed");
1804 REM_TYPE("dihre-fc");
1805 REM_TYPE("dihre-tau");
1806 REM_TYPE("nstdihreout");
1807 REM_TYPE("nstcheckpoint");
1808 REM_TYPE("optimize-fft");
1810 /* replace the following commands with the clearer new versions*/
1811 REPL_TYPE("unconstrained-start", "continuation");
1812 REPL_TYPE("foreign-lambda", "fep-lambdas");
1813 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1814 REPL_TYPE("nstxtcout", "nstxout-compressed");
1815 REPL_TYPE("xtc-grps", "compressed-x-grps");
1816 REPL_TYPE("xtc-precision", "compressed-x-precision");
1818 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1819 CTYPE ("Preprocessor information: use cpp syntax.");
1820 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1821 STYPE ("include", opts->include, NULL);
1822 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1823 STYPE ("define", opts->define, NULL);
1825 CCTYPE ("RUN CONTROL PARAMETERS");
1826 EETYPE("integrator", ir->eI, ei_names);
1827 CTYPE ("Start time and timestep in ps");
1828 RTYPE ("tinit", ir->init_t, 0.0);
1829 RTYPE ("dt", ir->delta_t, 0.001);
1830 STEPTYPE ("nsteps", ir->nsteps, 0);
1831 CTYPE ("For exact run continuation or redoing part of a run");
1832 STEPTYPE ("init-step", ir->init_step, 0);
1833 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1834 ITYPE ("simulation-part", ir->simulation_part, 1);
1835 CTYPE ("mode for center of mass motion removal");
1836 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1837 CTYPE ("number of steps for center of mass motion removal");
1838 ITYPE ("nstcomm", ir->nstcomm, 100);
1839 CTYPE ("group(s) for center of mass motion removal");
1840 STYPE ("comm-grps", is->vcm, NULL);
1842 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1843 CTYPE ("Friction coefficient (amu/ps) and random seed");
1844 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1845 STEPTYPE ("ld-seed", ir->ld_seed, -1);
1848 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1849 CTYPE ("Force tolerance and initial step-size");
1850 RTYPE ("emtol", ir->em_tol, 10.0);
1851 RTYPE ("emstep", ir->em_stepsize, 0.01);
1852 CTYPE ("Max number of iterations in relax-shells");
1853 ITYPE ("niter", ir->niter, 20);
1854 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1855 RTYPE ("fcstep", ir->fc_stepsize, 0);
1856 CTYPE ("Frequency of steepest descents steps when doing CG");
1857 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1858 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1860 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1861 RTYPE ("rtpi", ir->rtpi, 0.05);
1863 /* Output options */
1864 CCTYPE ("OUTPUT CONTROL OPTIONS");
1865 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1866 ITYPE ("nstxout", ir->nstxout, 0);
1867 ITYPE ("nstvout", ir->nstvout, 0);
1868 ITYPE ("nstfout", ir->nstfout, 0);
1869 CTYPE ("Output frequency for energies to log file and energy file");
1870 ITYPE ("nstlog", ir->nstlog, 1000);
1871 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1872 ITYPE ("nstenergy", ir->nstenergy, 1000);
1873 CTYPE ("Output frequency and precision for .xtc file");
1874 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
1875 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1876 CTYPE ("This selects the subset of atoms for the compressed");
1877 CTYPE ("trajectory file. You can select multiple groups. By");
1878 CTYPE ("default, all atoms will be written.");
1879 STYPE ("compressed-x-grps", is->x_compressed_groups, NULL);
1880 CTYPE ("Selection of energy groups");
1881 STYPE ("energygrps", is->energy, NULL);
1883 /* Neighbor searching */
1884 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1885 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1886 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1887 CTYPE ("nblist update frequency");
1888 ITYPE ("nstlist", ir->nstlist, 10);
1889 CTYPE ("ns algorithm (simple or grid)");
1890 EETYPE("ns-type", ir->ns_type, ens_names);
1891 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1892 EETYPE("pbc", ir->ePBC, epbc_names);
1893 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1894 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1895 CTYPE ("a value of -1 means: use rlist");
1896 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1897 CTYPE ("nblist cut-off");
1898 RTYPE ("rlist", ir->rlist, 1.0);
1899 CTYPE ("long-range cut-off for switched potentials");
1900 RTYPE ("rlistlong", ir->rlistlong, -1);
1901 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1903 /* Electrostatics */
1904 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1905 CTYPE ("Method for doing electrostatics");
1906 EETYPE("coulombtype", ir->coulombtype, eel_names);
1907 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1908 CTYPE ("cut-off lengths");
1909 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1910 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1911 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1912 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1913 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1914 CTYPE ("Method for doing Van der Waals");
1915 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1916 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1917 CTYPE ("cut-off lengths");
1918 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1919 RTYPE ("rvdw", ir->rvdw, 1.0);
1920 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1921 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1922 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1923 RTYPE ("table-extension", ir->tabext, 1.0);
1924 CTYPE ("Separate tables between energy group pairs");
1925 STYPE ("energygrp-table", is->egptable, NULL);
1926 CTYPE ("Spacing for the PME/PPPM FFT grid");
1927 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1928 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1929 ITYPE ("fourier-nx", ir->nkx, 0);
1930 ITYPE ("fourier-ny", ir->nky, 0);
1931 ITYPE ("fourier-nz", ir->nkz, 0);
1932 CTYPE ("EWALD/PME/PPPM parameters");
1933 ITYPE ("pme-order", ir->pme_order, 4);
1934 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1935 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1936 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1937 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1938 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1940 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1941 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1943 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1944 CTYPE ("Algorithm for calculating Born radii");
1945 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1946 CTYPE ("Frequency of calculating the Born radii inside rlist");
1947 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1948 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1949 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1950 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1951 CTYPE ("Dielectric coefficient of the implicit solvent");
1952 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1953 CTYPE ("Salt concentration in M for Generalized Born models");
1954 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1955 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1956 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1957 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1958 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1959 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1960 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1961 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1962 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1963 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1965 /* Coupling stuff */
1966 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1967 CTYPE ("Temperature coupling");
1968 EETYPE("tcoupl", ir->etc, etcoupl_names);
1969 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1970 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
1971 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1972 CTYPE ("Groups to couple separately");
1973 STYPE ("tc-grps", is->tcgrps, NULL);
1974 CTYPE ("Time constant (ps) and reference temperature (K)");
1975 STYPE ("tau-t", is->tau_t, NULL);
1976 STYPE ("ref-t", is->ref_t, NULL);
1977 CTYPE ("pressure coupling");
1978 EETYPE("pcoupl", ir->epc, epcoupl_names);
1979 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1980 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1981 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1982 RTYPE ("tau-p", ir->tau_p, 1.0);
1983 STYPE ("compressibility", dumstr[0], NULL);
1984 STYPE ("ref-p", dumstr[1], NULL);
1985 CTYPE ("Scaling of reference coordinates, No, All or COM");
1986 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1989 CCTYPE ("OPTIONS FOR QMMM calculations");
1990 EETYPE("QMMM", ir->bQMMM, yesno_names);
1991 CTYPE ("Groups treated Quantum Mechanically");
1992 STYPE ("QMMM-grps", is->QMMM, NULL);
1993 CTYPE ("QM method");
1994 STYPE("QMmethod", is->QMmethod, NULL);
1995 CTYPE ("QMMM scheme");
1996 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1997 CTYPE ("QM basisset");
1998 STYPE("QMbasis", is->QMbasis, NULL);
1999 CTYPE ("QM charge");
2000 STYPE ("QMcharge", is->QMcharge, NULL);
2001 CTYPE ("QM multiplicity");
2002 STYPE ("QMmult", is->QMmult, NULL);
2003 CTYPE ("Surface Hopping");
2004 STYPE ("SH", is->bSH, NULL);
2005 CTYPE ("CAS space options");
2006 STYPE ("CASorbitals", is->CASorbitals, NULL);
2007 STYPE ("CASelectrons", is->CASelectrons, NULL);
2008 STYPE ("SAon", is->SAon, NULL);
2009 STYPE ("SAoff", is->SAoff, NULL);
2010 STYPE ("SAsteps", is->SAsteps, NULL);
2011 CTYPE ("Scale factor for MM charges");
2012 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
2013 CTYPE ("Optimization of QM subsystem");
2014 STYPE ("bOPT", is->bOPT, NULL);
2015 STYPE ("bTS", is->bTS, NULL);
2017 /* Simulated annealing */
2018 CCTYPE("SIMULATED ANNEALING");
2019 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
2020 STYPE ("annealing", is->anneal, NULL);
2021 CTYPE ("Number of time points to use for specifying annealing in each group");
2022 STYPE ("annealing-npoints", is->anneal_npoints, NULL);
2023 CTYPE ("List of times at the annealing points for each group");
2024 STYPE ("annealing-time", is->anneal_time, NULL);
2025 CTYPE ("Temp. at each annealing point, for each group.");
2026 STYPE ("annealing-temp", is->anneal_temp, NULL);
2029 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
2030 EETYPE("gen-vel", opts->bGenVel, yesno_names);
2031 RTYPE ("gen-temp", opts->tempi, 300.0);
2032 ITYPE ("gen-seed", opts->seed, -1);
2035 CCTYPE ("OPTIONS FOR BONDS");
2036 EETYPE("constraints", opts->nshake, constraints);
2037 CTYPE ("Type of constraint algorithm");
2038 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
2039 CTYPE ("Do not constrain the start configuration");
2040 EETYPE("continuation", ir->bContinuation, yesno_names);
2041 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2042 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2043 CTYPE ("Relative tolerance of shake");
2044 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2045 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2046 ITYPE ("lincs-order", ir->nProjOrder, 4);
2047 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2048 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2049 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2050 ITYPE ("lincs-iter", ir->nLincsIter, 1);
2051 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2052 CTYPE ("rotates over more degrees than");
2053 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2054 CTYPE ("Convert harmonic bonds to morse potentials");
2055 EETYPE("morse", opts->bMorse, yesno_names);
2057 /* Energy group exclusions */
2058 CCTYPE ("ENERGY GROUP EXCLUSIONS");
2059 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2060 STYPE ("energygrp-excl", is->egpexcl, NULL);
2064 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2065 ITYPE ("nwall", ir->nwall, 0);
2066 EETYPE("wall-type", ir->wall_type, ewt_names);
2067 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2068 STYPE ("wall-atomtype", is->wall_atomtype, NULL);
2069 STYPE ("wall-density", is->wall_density, NULL);
2070 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2073 CCTYPE("COM PULLING");
2074 CTYPE("Pull type: no, umbrella, constraint or constant-force");
2075 EETYPE("pull", ir->ePull, epull_names);
2076 if (ir->ePull != epullNO)
2079 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
2082 /* Enforced rotation */
2083 CCTYPE("ENFORCED ROTATION");
2084 CTYPE("Enforced rotation: No or Yes");
2085 EETYPE("rotation", ir->bRot, yesno_names);
2089 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2092 /* Interactive MD */
2094 CCTYPE("Group to display and/or manipulate in interactive MD session");
2095 STYPE ("IMD-group", is->imd_grp, NULL);
2096 if (is->imd_grp[0] != '\0')
2103 CCTYPE("NMR refinement stuff");
2104 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2105 EETYPE("disre", ir->eDisre, edisre_names);
2106 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2107 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2108 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2109 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2110 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2111 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2112 CTYPE ("Output frequency for pair distances to energy file");
2113 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2114 CTYPE ("Orientation restraints: No or Yes");
2115 EETYPE("orire", opts->bOrire, yesno_names);
2116 CTYPE ("Orientation restraints force constant and tau for time averaging");
2117 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2118 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2119 STYPE ("orire-fitgrp", is->orirefitgrp, NULL);
2120 CTYPE ("Output frequency for trace(SD) and S to energy file");
2121 ITYPE ("nstorireout", ir->nstorireout, 100);
2123 /* free energy variables */
2124 CCTYPE ("Free energy variables");
2125 EETYPE("free-energy", ir->efep, efep_names);
2126 STYPE ("couple-moltype", is->couple_moltype, NULL);
2127 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2128 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2129 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2131 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2133 it was not entered */
2134 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2135 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2136 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2137 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
2138 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
2139 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
2140 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
2141 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
2142 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
2143 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
2144 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2145 STYPE ("init-lambda-weights", is->lambda_weights, NULL);
2146 EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
2147 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2148 ITYPE ("sc-power", fep->sc_power, 1);
2149 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2150 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2151 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2152 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2153 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2154 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2155 separate_dhdl_file_names);
2156 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2157 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2158 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2160 /* Non-equilibrium MD stuff */
2161 CCTYPE("Non-equilibrium MD stuff");
2162 STYPE ("acc-grps", is->accgrps, NULL);
2163 STYPE ("accelerate", is->acc, NULL);
2164 STYPE ("freezegrps", is->freeze, NULL);
2165 STYPE ("freezedim", is->frdim, NULL);
2166 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2167 STYPE ("deform", is->deform, NULL);
2169 /* simulated tempering variables */
2170 CCTYPE("simulated tempering variables");
2171 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2172 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2173 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2174 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2176 /* expanded ensemble variables */
2177 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2179 read_expandedparams(&ninp, &inp, expand, wi);
2182 /* Electric fields */
2183 CCTYPE("Electric fields");
2184 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2185 CTYPE ("and a phase angle (real)");
2186 STYPE ("E-x", is->efield_x, NULL);
2187 STYPE ("E-xt", is->efield_xt, NULL);
2188 STYPE ("E-y", is->efield_y, NULL);
2189 STYPE ("E-yt", is->efield_yt, NULL);
2190 STYPE ("E-z", is->efield_z, NULL);
2191 STYPE ("E-zt", is->efield_zt, NULL);
2193 CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2194 CTYPE("Swap positions along direction: no, X, Y, Z");
2195 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2196 if (ir->eSwapCoords != eswapNO)
2199 CTYPE("Swap attempt frequency");
2200 ITYPE("swap-frequency", ir->swap->nstswap, 1);
2201 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2202 STYPE("split-group0", splitgrp0, NULL);
2203 STYPE("split-group1", splitgrp1, NULL);
2204 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2205 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2206 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2208 CTYPE("Group name of ions that can be exchanged with solvent molecules");
2209 STYPE("swap-group", swapgrp, NULL);
2210 CTYPE("Group name of solvent molecules");
2211 STYPE("solvent-group", solgrp, NULL);
2213 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2214 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2215 CTYPE("however, if correctly defined, the ion permeation events are counted per channel");
2216 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2217 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2218 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2219 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2220 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2221 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2223 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2224 ITYPE("coupl-steps", ir->swap->nAverage, 10);
2225 CTYPE("Requested number of anions and cations for each of the two compartments");
2226 CTYPE("-1 means fix the numbers as found in time step 0");
2227 ITYPE("anionsA", ir->swap->nanions[0], -1);
2228 ITYPE("cationsA", ir->swap->ncations[0], -1);
2229 ITYPE("anionsB", ir->swap->nanions[1], -1);
2230 ITYPE("cationsB", ir->swap->ncations[1], -1);
2231 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2232 RTYPE("threshold", ir->swap->threshold, 1.0);
2235 /* AdResS defined thingies */
2236 CCTYPE ("AdResS parameters");
2237 EETYPE("adress", ir->bAdress, yesno_names);
2240 snew(ir->adress, 1);
2241 read_adressparams(&ninp, &inp, ir->adress, wi);
2244 /* User defined thingies */
2245 CCTYPE ("User defined thingies");
2246 STYPE ("user1-grps", is->user1, NULL);
2247 STYPE ("user2-grps", is->user2, NULL);
2248 ITYPE ("userint1", ir->userint1, 0);
2249 ITYPE ("userint2", ir->userint2, 0);
2250 ITYPE ("userint3", ir->userint3, 0);
2251 ITYPE ("userint4", ir->userint4, 0);
2252 RTYPE ("userreal1", ir->userreal1, 0);
2253 RTYPE ("userreal2", ir->userreal2, 0);
2254 RTYPE ("userreal3", ir->userreal3, 0);
2255 RTYPE ("userreal4", ir->userreal4, 0);
2258 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2259 for (i = 0; (i < ninp); i++)
2262 sfree(inp[i].value);
2266 /* Process options if necessary */
2267 for (m = 0; m < 2; m++)
2269 for (i = 0; i < 2*DIM; i++)
2278 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2280 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2282 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2284 case epctSEMIISOTROPIC:
2285 case epctSURFACETENSION:
2286 if (sscanf(dumstr[m], "%lf%lf",
2287 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2289 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2291 dumdub[m][YY] = dumdub[m][XX];
2293 case epctANISOTROPIC:
2294 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2295 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2296 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2298 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2302 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2303 epcoupltype_names[ir->epct]);
2307 clear_mat(ir->ref_p);
2308 clear_mat(ir->compress);
2309 for (i = 0; i < DIM; i++)
2311 ir->ref_p[i][i] = dumdub[1][i];
2312 ir->compress[i][i] = dumdub[0][i];
2314 if (ir->epct == epctANISOTROPIC)
2316 ir->ref_p[XX][YY] = dumdub[1][3];
2317 ir->ref_p[XX][ZZ] = dumdub[1][4];
2318 ir->ref_p[YY][ZZ] = dumdub[1][5];
2319 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2321 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2323 ir->compress[XX][YY] = dumdub[0][3];
2324 ir->compress[XX][ZZ] = dumdub[0][4];
2325 ir->compress[YY][ZZ] = dumdub[0][5];
2326 for (i = 0; i < DIM; i++)
2328 for (m = 0; m < i; m++)
2330 ir->ref_p[i][m] = ir->ref_p[m][i];
2331 ir->compress[i][m] = ir->compress[m][i];
2336 if (ir->comm_mode == ecmNO)
2341 opts->couple_moltype = NULL;
2342 if (strlen(is->couple_moltype) > 0)
2344 if (ir->efep != efepNO)
2346 opts->couple_moltype = gmx_strdup(is->couple_moltype);
2347 if (opts->couple_lam0 == opts->couple_lam1)
2349 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2351 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2352 opts->couple_lam1 == ecouplamNONE))
2354 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2359 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2362 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2363 if (ir->efep != efepNO)
2365 if (fep->delta_lambda > 0)
2367 ir->efep = efepSLOWGROWTH;
2371 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2373 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2374 warning_note(wi, "Old option for dhdl-print-energy given: "
2375 "changing \"yes\" to \"total\"\n");
2378 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2380 /* always print out the energy to dhdl if we are doing
2381 expanded ensemble, since we need the total energy for
2382 analysis if the temperature is changing. In some
2383 conditions one may only want the potential energy, so
2384 we will allow that if the appropriate mdp setting has
2385 been enabled. Otherwise, total it is:
2387 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2390 if ((ir->efep != efepNO) || ir->bSimTemp)
2392 ir->bExpanded = FALSE;
2393 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2395 ir->bExpanded = TRUE;
2397 do_fep_params(ir, is->fep_lambda, is->lambda_weights);
2398 if (ir->bSimTemp) /* done after fep params */
2400 do_simtemp_params(ir);
2405 ir->fepvals->n_lambda = 0;
2408 /* WALL PARAMETERS */
2410 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2412 /* ORIENTATION RESTRAINT PARAMETERS */
2414 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
2416 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2419 /* DEFORMATION PARAMETERS */
2421 clear_mat(ir->deform);
2422 for (i = 0; i < 6; i++)
2426 m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2427 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2428 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2429 for (i = 0; i < 3; i++)
2431 ir->deform[i][i] = dumdub[0][i];
2433 ir->deform[YY][XX] = dumdub[0][3];
2434 ir->deform[ZZ][XX] = dumdub[0][4];
2435 ir->deform[ZZ][YY] = dumdub[0][5];
2436 if (ir->epc != epcNO)
2438 for (i = 0; i < 3; i++)
2440 for (j = 0; j <= i; j++)
2442 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2444 warning_error(wi, "A box element has deform set and compressibility > 0");
2448 for (i = 0; i < 3; i++)
2450 for (j = 0; j < i; j++)
2452 if (ir->deform[i][j] != 0)
2454 for (m = j; m < DIM; m++)
2456 if (ir->compress[m][j] != 0)
2458 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.");
2459 warning(wi, warn_buf);
2467 /* Ion/water position swapping checks */
2468 if (ir->eSwapCoords != eswapNO)
2470 if (ir->swap->nstswap < 1)
2472 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2474 if (ir->swap->nAverage < 1)
2476 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2478 if (ir->swap->threshold < 1.0)
2480 warning_error(wi, "Ion count threshold must be at least 1.\n");
2488 static int search_QMstring(const char *s, int ng, const char *gn[])
2490 /* same as normal search_string, but this one searches QM strings */
2493 for (i = 0; (i < ng); i++)
2495 if (gmx_strcasecmp(s, gn[i]) == 0)
2501 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2505 } /* search_QMstring */
2507 /* We would like gn to be const as well, but C doesn't allow this */
2508 int search_string(const char *s, int ng, char *gn[])
2512 for (i = 0; (i < ng); i++)
2514 if (gmx_strcasecmp(s, gn[i]) == 0)
2521 "Group %s referenced in the .mdp file was not found in the index file.\n"
2522 "Group names must match either [moleculetype] names or custom index group\n"
2523 "names, in which case you must supply an index file to the '-n' option\n"
2530 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2531 t_blocka *block, char *gnames[],
2532 int gtype, int restnm,
2533 int grptp, gmx_bool bVerbose,
2536 unsigned short *cbuf;
2537 t_grps *grps = &(groups->grps[gtype]);
2538 int i, j, gid, aj, ognr, ntot = 0;
2541 char warn_buf[STRLEN];
2545 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2548 title = gtypes[gtype];
2551 /* Mark all id's as not set */
2552 for (i = 0; (i < natoms); i++)
2557 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2558 for (i = 0; (i < ng); i++)
2560 /* Lookup the group name in the block structure */
2561 gid = search_string(ptrs[i], block->nr, gnames);
2562 if ((grptp != egrptpONE) || (i == 0))
2564 grps->nm_ind[grps->nr++] = gid;
2568 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2571 /* Now go over the atoms in the group */
2572 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2577 /* Range checking */
2578 if ((aj < 0) || (aj >= natoms))
2580 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2582 /* Lookup up the old group number */
2586 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2587 aj+1, title, ognr+1, i+1);
2591 /* Store the group number in buffer */
2592 if (grptp == egrptpONE)
2605 /* Now check whether we have done all atoms */
2609 if (grptp == egrptpALL)
2611 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2612 natoms-ntot, title);
2614 else if (grptp == egrptpPART)
2616 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2617 natoms-ntot, title);
2618 warning_note(wi, warn_buf);
2620 /* Assign all atoms currently unassigned to a rest group */
2621 for (j = 0; (j < natoms); j++)
2623 if (cbuf[j] == NOGID)
2629 if (grptp != egrptpPART)
2634 "Making dummy/rest group for %s containing %d elements\n",
2635 title, natoms-ntot);
2637 /* Add group name "rest" */
2638 grps->nm_ind[grps->nr] = restnm;
2640 /* Assign the rest name to all atoms not currently assigned to a group */
2641 for (j = 0; (j < natoms); j++)
2643 if (cbuf[j] == NOGID)
2652 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2654 /* All atoms are part of one (or no) group, no index required */
2655 groups->ngrpnr[gtype] = 0;
2656 groups->grpnr[gtype] = NULL;
2660 groups->ngrpnr[gtype] = natoms;
2661 snew(groups->grpnr[gtype], natoms);
2662 for (j = 0; (j < natoms); j++)
2664 groups->grpnr[gtype][j] = cbuf[j];
2670 return (bRest && grptp == egrptpPART);
2673 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2676 gmx_groups_t *groups;
2678 int natoms, ai, aj, i, j, d, g, imin, jmin;
2680 int *nrdf2, *na_vcm, na_tot;
2681 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2682 gmx_mtop_atomloop_all_t aloop;
2684 int mb, mol, ftype, as;
2685 gmx_molblock_t *molb;
2686 gmx_moltype_t *molt;
2689 * First calc 3xnr-atoms for each group
2690 * then subtract half a degree of freedom for each constraint
2692 * Only atoms and nuclei contribute to the degrees of freedom...
2697 groups = &mtop->groups;
2698 natoms = mtop->natoms;
2700 /* Allocate one more for a possible rest group */
2701 /* We need to sum degrees of freedom into doubles,
2702 * since floats give too low nrdf's above 3 million atoms.
2704 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2705 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2706 snew(na_vcm, groups->grps[egcVCM].nr+1);
2708 for (i = 0; i < groups->grps[egcTC].nr; i++)
2712 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2717 snew(nrdf2, natoms);
2718 aloop = gmx_mtop_atomloop_all_init(mtop);
2719 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2722 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2724 g = ggrpnr(groups, egcFREEZE, i);
2725 /* Double count nrdf for particle i */
2726 for (d = 0; d < DIM; d++)
2728 if (opts->nFreeze[g][d] == 0)
2733 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2734 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2739 for (mb = 0; mb < mtop->nmolblock; mb++)
2741 molb = &mtop->molblock[mb];
2742 molt = &mtop->moltype[molb->type];
2743 atom = molt->atoms.atom;
2744 for (mol = 0; mol < molb->nmol; mol++)
2746 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2748 ia = molt->ilist[ftype].iatoms;
2749 for (i = 0; i < molt->ilist[ftype].nr; )
2751 /* Subtract degrees of freedom for the constraints,
2752 * if the particles still have degrees of freedom left.
2753 * If one of the particles is a vsite or a shell, then all
2754 * constraint motion will go there, but since they do not
2755 * contribute to the constraints the degrees of freedom do not
2760 if (((atom[ia[1]].ptype == eptNucleus) ||
2761 (atom[ia[1]].ptype == eptAtom)) &&
2762 ((atom[ia[2]].ptype == eptNucleus) ||
2763 (atom[ia[2]].ptype == eptAtom)))
2781 imin = min(imin, nrdf2[ai]);
2782 jmin = min(jmin, nrdf2[aj]);
2785 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2786 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2787 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2788 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2790 ia += interaction_function[ftype].nratoms+1;
2791 i += interaction_function[ftype].nratoms+1;
2794 ia = molt->ilist[F_SETTLE].iatoms;
2795 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2797 /* Subtract 1 dof from every atom in the SETTLE */
2798 for (j = 0; j < 3; j++)
2801 imin = min(2, nrdf2[ai]);
2803 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2804 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2809 as += molt->atoms.nr;
2813 if (ir->ePull == epullCONSTRAINT)
2815 /* Correct nrdf for the COM constraints.
2816 * We correct using the TC and VCM group of the first atom
2817 * in the reference and pull group. If atoms in one pull group
2818 * belong to different TC or VCM groups it is anyhow difficult
2819 * to determine the optimal nrdf assignment.
2823 for (i = 0; i < pull->ncoord; i++)
2827 for (j = 0; j < 2; j++)
2829 const t_pull_group *pgrp;
2831 pgrp = &pull->group[pull->coord[i].group[j]];
2835 /* Subtract 1/2 dof from each group */
2837 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2838 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2839 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2841 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)]]);
2846 /* We need to subtract the whole DOF from group j=1 */
2853 if (ir->nstcomm != 0)
2855 /* Subtract 3 from the number of degrees of freedom in each vcm group
2856 * when com translation is removed and 6 when rotation is removed
2859 switch (ir->comm_mode)
2862 n_sub = ndof_com(ir);
2869 gmx_incons("Checking comm_mode");
2872 for (i = 0; i < groups->grps[egcTC].nr; i++)
2874 /* Count the number of atoms of TC group i for every VCM group */
2875 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2880 for (ai = 0; ai < natoms; ai++)
2882 if (ggrpnr(groups, egcTC, ai) == i)
2884 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2888 /* Correct for VCM removal according to the fraction of each VCM
2889 * group present in this TC group.
2891 nrdf_uc = nrdf_tc[i];
2894 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2898 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2900 if (nrdf_vcm[j] > n_sub)
2902 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2903 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2907 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2908 j, nrdf_vcm[j], nrdf_tc[i]);
2913 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2915 opts->nrdf[i] = nrdf_tc[i];
2916 if (opts->nrdf[i] < 0)
2921 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2922 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2931 static void decode_cos(char *s, t_cosines *cosine)
2934 char format[STRLEN], f1[STRLEN];
2946 sscanf(t, "%d", &(cosine->n));
2953 snew(cosine->a, cosine->n);
2954 snew(cosine->phi, cosine->n);
2956 sprintf(format, "%%*d");
2957 for (i = 0; (i < cosine->n); i++)
2960 strcat(f1, "%lf%lf");
2961 if (sscanf(t, f1, &a, &phi) < 2)
2963 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2966 cosine->phi[i] = phi;
2967 strcat(format, "%*lf%*lf");
2974 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2975 const char *option, const char *val, int flag)
2977 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2978 * But since this is much larger than STRLEN, such a line can not be parsed.
2979 * The real maximum is the number of names that fit in a string: STRLEN/2.
2981 #define EGP_MAX (STRLEN/2)
2982 int nelem, i, j, k, nr;
2983 char *names[EGP_MAX];
2987 gnames = groups->grpname;
2989 nelem = str_nelem(val, EGP_MAX, names);
2992 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2994 nr = groups->grps[egcENER].nr;
2996 for (i = 0; i < nelem/2; i++)
3000 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
3006 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3007 names[2*i], option);
3011 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3017 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3018 names[2*i+1], option);
3020 if ((j < nr) && (k < nr))
3022 ir->opts.egp_flags[nr*j+k] |= flag;
3023 ir->opts.egp_flags[nr*k+j] |= flag;
3032 static void make_swap_groups(
3041 int ig = -1, i = 0, j;
3045 /* Just a quick check here, more thorough checks are in mdrun */
3046 if (strcmp(splitg0name, splitg1name) == 0)
3048 gmx_fatal(FARGS, "The split groups can not both be '%s'.", splitg0name);
3051 /* First get the swap group index atoms */
3052 ig = search_string(swapgname, grps->nr, gnames);
3053 swap->nat = grps->index[ig+1] - grps->index[ig];
3056 fprintf(stderr, "Swap group '%s' contains %d atoms.\n", swapgname, swap->nat);
3057 snew(swap->ind, swap->nat);
3058 for (i = 0; i < swap->nat; i++)
3060 swap->ind[i] = grps->a[grps->index[ig]+i];
3065 gmx_fatal(FARGS, "You defined an empty group of atoms for swapping.");
3068 /* Now do so for the split groups */
3069 for (j = 0; j < 2; j++)
3073 splitg = splitg0name;
3077 splitg = splitg1name;
3080 ig = search_string(splitg, grps->nr, gnames);
3081 swap->nat_split[j] = grps->index[ig+1] - grps->index[ig];
3082 if (swap->nat_split[j] > 0)
3084 fprintf(stderr, "Split group %d '%s' contains %d atom%s.\n",
3085 j, splitg, swap->nat_split[j], (swap->nat_split[j] > 1) ? "s" : "");
3086 snew(swap->ind_split[j], swap->nat_split[j]);
3087 for (i = 0; i < swap->nat_split[j]; i++)
3089 swap->ind_split[j][i] = grps->a[grps->index[ig]+i];
3094 gmx_fatal(FARGS, "Split group %d has to contain at least 1 atom!", j);
3098 /* Now get the solvent group index atoms */
3099 ig = search_string(solgname, grps->nr, gnames);
3100 swap->nat_sol = grps->index[ig+1] - grps->index[ig];
3101 if (swap->nat_sol > 0)
3103 fprintf(stderr, "Solvent group '%s' contains %d atoms.\n", solgname, swap->nat_sol);
3104 snew(swap->ind_sol, swap->nat_sol);
3105 for (i = 0; i < swap->nat_sol; i++)
3107 swap->ind_sol[i] = grps->a[grps->index[ig]+i];
3112 gmx_fatal(FARGS, "You defined an empty group of solvent. Cannot exchange ions.");
3117 void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3122 ig = search_string(IMDgname, grps->nr, gnames);
3123 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3125 if (IMDgroup->nat > 0)
3127 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3128 IMDgname, IMDgroup->nat);
3129 snew(IMDgroup->ind, IMDgroup->nat);
3130 for (i = 0; i < IMDgroup->nat; i++)
3132 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3138 void do_index(const char* mdparin, const char *ndx,
3141 t_inputrec *ir, rvec *v,
3145 gmx_groups_t *groups;
3149 char warnbuf[STRLEN], **gnames;
3150 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3153 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3154 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3155 int i, j, k, restnm;
3157 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3158 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
3159 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
3160 char warn_buf[STRLEN];
3164 fprintf(stderr, "processing index file...\n");
3170 snew(grps->index, 1);
3172 atoms_all = gmx_mtop_global_atoms(mtop);
3173 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3174 free_t_atoms(&atoms_all, FALSE);
3178 grps = init_index(ndx, &gnames);
3181 groups = &mtop->groups;
3182 natoms = mtop->natoms;
3183 symtab = &mtop->symtab;
3185 snew(groups->grpname, grps->nr+1);
3187 for (i = 0; (i < grps->nr); i++)
3189 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3191 groups->grpname[i] = put_symtab(symtab, "rest");
3193 srenew(gnames, grps->nr+1);
3194 gnames[restnm] = *(groups->grpname[i]);
3195 groups->ngrpname = grps->nr+1;
3197 set_warning_line(wi, mdparin, -1);
3199 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3200 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3201 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3202 if ((ntau_t != ntcg) || (nref_t != ntcg))
3204 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3205 "%d tau-t values", ntcg, nref_t, ntau_t);
3208 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3209 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3210 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3211 nr = groups->grps[egcTC].nr;
3213 snew(ir->opts.nrdf, nr);
3214 snew(ir->opts.tau_t, nr);
3215 snew(ir->opts.ref_t, nr);
3216 if (ir->eI == eiBD && ir->bd_fric == 0)
3218 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3225 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3229 for (i = 0; (i < nr); i++)
3231 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
3232 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
3234 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3235 warning_error(wi, warn_buf);
3238 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3240 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.");
3243 if (ir->opts.tau_t[i] >= 0)
3245 tau_min = min(tau_min, ir->opts.tau_t[i]);
3248 if (ir->etc != etcNO && ir->nsttcouple == -1)
3250 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3255 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3257 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");
3259 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
3261 if (ir->nstpcouple != ir->nsttcouple)
3263 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
3264 ir->nstpcouple = ir->nsttcouple = mincouple;
3265 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);
3266 warning_note(wi, warn_buf);
3270 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3271 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3273 if (ETC_ANDERSEN(ir->etc))
3275 if (ir->nsttcouple != 1)
3278 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");
3279 warning_note(wi, warn_buf);
3282 nstcmin = tcouple_min_integration_steps(ir->etc);
3285 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3287 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3288 ETCOUPLTYPE(ir->etc),
3290 ir->nsttcouple*ir->delta_t);
3291 warning(wi, warn_buf);
3294 for (i = 0; (i < nr); i++)
3296 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3297 if (ir->opts.ref_t[i] < 0)
3299 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3302 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3303 if we are in this conditional) if mc_temp is negative */
3304 if (ir->expandedvals->mc_temp < 0)
3306 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3310 /* Simulated annealing for each group. There are nr groups */
3311 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3312 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3316 if (nSA > 0 && nSA != nr)
3318 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3322 snew(ir->opts.annealing, nr);
3323 snew(ir->opts.anneal_npoints, nr);
3324 snew(ir->opts.anneal_time, nr);
3325 snew(ir->opts.anneal_temp, nr);
3326 for (i = 0; i < nr; i++)
3328 ir->opts.annealing[i] = eannNO;
3329 ir->opts.anneal_npoints[i] = 0;
3330 ir->opts.anneal_time[i] = NULL;
3331 ir->opts.anneal_temp[i] = NULL;
3336 for (i = 0; i < nr; i++)
3338 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3340 ir->opts.annealing[i] = eannNO;
3342 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3344 ir->opts.annealing[i] = eannSINGLE;
3347 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3349 ir->opts.annealing[i] = eannPERIODIC;
3355 /* Read the other fields too */
3356 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3357 if (nSA_points != nSA)
3359 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3361 for (k = 0, i = 0; i < nr; i++)
3363 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3364 if (ir->opts.anneal_npoints[i] == 1)
3366 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3368 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3369 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3370 k += ir->opts.anneal_npoints[i];
3373 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3376 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3378 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3381 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3384 for (i = 0, k = 0; i < nr; i++)
3387 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3389 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3390 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3393 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3395 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3401 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3403 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3404 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3407 if (ir->opts.anneal_temp[i][j] < 0)
3409 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3414 /* Print out some summary information, to make sure we got it right */
3415 for (i = 0, k = 0; i < nr; i++)
3417 if (ir->opts.annealing[i] != eannNO)
3419 j = groups->grps[egcTC].nm_ind[i];
3420 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3421 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3422 ir->opts.anneal_npoints[i]);
3423 fprintf(stderr, "Time (ps) Temperature (K)\n");
3424 /* All terms except the last one */
3425 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3427 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3430 /* Finally the last one */
3431 j = ir->opts.anneal_npoints[i]-1;
3432 if (ir->opts.annealing[i] == eannSINGLE)
3434 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3438 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3439 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3441 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3450 if (ir->ePull != epullNO)
3452 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3454 make_pull_coords(ir->pull);
3459 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3462 if (ir->eSwapCoords != eswapNO)
3464 make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
3467 /* Make indices for IMD session */
3470 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3473 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3474 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3475 if (nacg*DIM != nacc)
3477 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3480 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3481 restnm, egrptpALL_GENREST, bVerbose, wi);
3482 nr = groups->grps[egcACC].nr;
3483 snew(ir->opts.acc, nr);
3484 ir->opts.ngacc = nr;
3486 for (i = k = 0; (i < nacg); i++)
3488 for (j = 0; (j < DIM); j++, k++)
3490 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3493 for (; (i < nr); i++)
3495 for (j = 0; (j < DIM); j++)
3497 ir->opts.acc[i][j] = 0;
3501 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3502 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3503 if (nfrdim != DIM*nfreeze)
3505 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3508 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3509 restnm, egrptpALL_GENREST, bVerbose, wi);
3510 nr = groups->grps[egcFREEZE].nr;
3511 ir->opts.ngfrz = nr;
3512 snew(ir->opts.nFreeze, nr);
3513 for (i = k = 0; (i < nfreeze); i++)
3515 for (j = 0; (j < DIM); j++, k++)
3517 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3518 if (!ir->opts.nFreeze[i][j])
3520 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3522 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3523 "(not %s)", ptr1[k]);
3524 warning(wi, warn_buf);
3529 for (; (i < nr); i++)
3531 for (j = 0; (j < DIM); j++)
3533 ir->opts.nFreeze[i][j] = 0;
3537 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3538 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3539 restnm, egrptpALL_GENREST, bVerbose, wi);
3540 add_wall_energrps(groups, ir->nwall, symtab);
3541 ir->opts.ngener = groups->grps[egcENER].nr;
3542 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3544 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3545 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3548 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3549 "This may lead to artifacts.\n"
3550 "In most cases one should use one group for the whole system.");
3553 /* Now we have filled the freeze struct, so we can calculate NRDF */
3554 calc_nrdf(mtop, ir, gnames);
3560 /* Must check per group! */
3561 for (i = 0; (i < ir->opts.ngtc); i++)
3563 ntot += ir->opts.nrdf[i];
3565 if (ntot != (DIM*natoms))
3567 fac = sqrt(ntot/(DIM*natoms));
3570 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3571 "and removal of center of mass motion\n", fac);
3573 for (i = 0; (i < natoms); i++)
3575 svmul(fac, v[i], v[i]);
3580 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3581 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3582 restnm, egrptpALL_GENREST, bVerbose, wi);
3583 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3584 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3585 restnm, egrptpALL_GENREST, bVerbose, wi);
3586 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3587 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3588 restnm, egrptpONE, bVerbose, wi);
3589 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3590 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3591 restnm, egrptpALL_GENREST, bVerbose, wi);
3593 /* QMMM input processing */
3594 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3595 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3596 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3597 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3599 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3600 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3602 /* group rest, if any, is always MM! */
3603 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3604 restnm, egrptpALL_GENREST, bVerbose, wi);
3605 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3606 ir->opts.ngQM = nQMg;
3607 snew(ir->opts.QMmethod, nr);
3608 snew(ir->opts.QMbasis, nr);
3609 for (i = 0; i < nr; i++)
3611 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3612 * converted to the corresponding enum in names.c
3614 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3616 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3620 nQMmult = str_nelem(is->QMmult, MAXPTR, ptr1);
3621 nQMcharge = str_nelem(is->QMcharge, MAXPTR, ptr2);
3622 nbSH = str_nelem(is->bSH, MAXPTR, ptr3);
3623 snew(ir->opts.QMmult, nr);
3624 snew(ir->opts.QMcharge, nr);
3625 snew(ir->opts.bSH, nr);
3627 for (i = 0; i < nr; i++)
3629 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3630 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3631 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3634 nCASelec = str_nelem(is->CASelectrons, MAXPTR, ptr1);
3635 nCASorb = str_nelem(is->CASorbitals, MAXPTR, ptr2);
3636 snew(ir->opts.CASelectrons, nr);
3637 snew(ir->opts.CASorbitals, nr);
3638 for (i = 0; i < nr; i++)
3640 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3641 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3643 /* special optimization options */
3645 nbOPT = str_nelem(is->bOPT, MAXPTR, ptr1);
3646 nbTS = str_nelem(is->bTS, MAXPTR, ptr2);
3647 snew(ir->opts.bOPT, nr);
3648 snew(ir->opts.bTS, nr);
3649 for (i = 0; i < nr; i++)
3651 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3652 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3654 nSAon = str_nelem(is->SAon, MAXPTR, ptr1);
3655 nSAoff = str_nelem(is->SAoff, MAXPTR, ptr2);
3656 nSAsteps = str_nelem(is->SAsteps, MAXPTR, ptr3);
3657 snew(ir->opts.SAon, nr);
3658 snew(ir->opts.SAoff, nr);
3659 snew(ir->opts.SAsteps, nr);
3661 for (i = 0; i < nr; i++)
3663 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3664 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3665 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3667 /* end of QMMM input */
3671 for (i = 0; (i < egcNR); i++)
3673 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3674 for (j = 0; (j < groups->grps[i].nr); j++)
3676 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3678 fprintf(stderr, "\n");
3682 nr = groups->grps[egcENER].nr;
3683 snew(ir->opts.egp_flags, nr*nr);
3685 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3686 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3688 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3690 if (bExcl && EEL_FULL(ir->coulombtype))
3692 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3695 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3696 if (bTable && !(ir->vdwtype == evdwUSER) &&
3697 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3698 !(ir->coulombtype == eelPMEUSERSWITCH))
3700 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3703 decode_cos(is->efield_x, &(ir->ex[XX]));
3704 decode_cos(is->efield_xt, &(ir->et[XX]));
3705 decode_cos(is->efield_y, &(ir->ex[YY]));
3706 decode_cos(is->efield_yt, &(ir->et[YY]));
3707 decode_cos(is->efield_z, &(ir->ex[ZZ]));
3708 decode_cos(is->efield_zt, &(ir->et[ZZ]));
3712 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3715 for (i = 0; (i < grps->nr); i++)
3727 static void check_disre(gmx_mtop_t *mtop)
3729 gmx_ffparams_t *ffparams;
3730 t_functype *functype;
3732 int i, ndouble, ftype;
3733 int label, old_label;
3735 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3737 ffparams = &mtop->ffparams;
3738 functype = ffparams->functype;
3739 ip = ffparams->iparams;
3742 for (i = 0; i < ffparams->ntypes; i++)
3744 ftype = functype[i];
3745 if (ftype == F_DISRES)
3747 label = ip[i].disres.label;
3748 if (label == old_label)
3750 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3758 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3759 "probably the parameters for multiple pairs in one restraint "
3760 "are not identical\n", ndouble);
3765 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3766 gmx_bool posres_only,
3770 gmx_mtop_ilistloop_t iloop;
3780 for (d = 0; d < DIM; d++)
3782 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3784 /* Check for freeze groups */
3785 for (g = 0; g < ir->opts.ngfrz; g++)
3787 for (d = 0; d < DIM; d++)
3789 if (ir->opts.nFreeze[g][d] != 0)
3797 /* Check for position restraints */
3798 iloop = gmx_mtop_ilistloop_init(sys);
3799 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3802 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3804 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3806 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3807 for (d = 0; d < DIM; d++)
3809 if (pr->posres.fcA[d] != 0)
3815 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3817 /* Check for flat-bottom posres */
3818 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3819 if (pr->fbposres.k != 0)
3821 switch (pr->fbposres.geom)
3823 case efbposresSPHERE:
3824 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3826 case efbposresCYLINDER:
3827 AbsRef[XX] = AbsRef[YY] = 1;
3829 case efbposresX: /* d=XX */
3830 case efbposresY: /* d=YY */
3831 case efbposresZ: /* d=ZZ */
3832 d = pr->fbposres.geom - efbposresX;
3836 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3837 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3845 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3849 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3850 gmx_bool *bC6ParametersWorkWithGeometricRules,
3851 gmx_bool *bC6ParametersWorkWithLBRules,
3852 gmx_bool *bLBRulesPossible)
3854 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3857 double geometricdiff, LBdiff;
3858 double c6i, c6j, c12i, c12j;
3859 double c6, c6_geometric, c6_LB;
3860 double sigmai, sigmaj, epsi, epsj;
3861 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3864 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3865 * force-field floating point parameters.
3868 ptr = getenv("GMX_LJCOMB_TOL");
3873 sscanf(ptr, "%lf", &dbl);
3877 *bC6ParametersWorkWithLBRules = TRUE;
3878 *bC6ParametersWorkWithGeometricRules = TRUE;
3879 bCanDoLBRules = TRUE;
3880 bCanDoGeometricRules = TRUE;
3881 ntypes = mtop->ffparams.atnr;
3882 snew(typecount, ntypes);
3883 gmx_mtop_count_atomtypes(mtop, state, typecount);
3884 geometricdiff = LBdiff = 0.0;
3885 *bLBRulesPossible = TRUE;
3886 for (tpi = 0; tpi < ntypes; ++tpi)
3888 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3889 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3890 for (tpj = tpi; tpj < ntypes; ++tpj)
3892 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3893 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3894 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3895 c6_geometric = sqrt(c6i * c6j);
3896 if (!gmx_numzero(c6_geometric))
3898 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3900 sigmai = pow(c12i / c6i, 1.0/6.0);
3901 sigmaj = pow(c12j / c6j, 1.0/6.0);
3902 epsi = c6i * c6i /(4.0 * c12i);
3903 epsj = c6j * c6j /(4.0 * c12j);
3904 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3908 *bLBRulesPossible = FALSE;
3909 c6_LB = c6_geometric;
3911 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3914 if (FALSE == bCanDoLBRules)
3916 *bC6ParametersWorkWithLBRules = FALSE;
3919 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3921 if (FALSE == bCanDoGeometricRules)
3923 *bC6ParametersWorkWithGeometricRules = FALSE;
3931 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3935 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3937 check_combination_rule_differences(mtop, 0,
3938 &bC6ParametersWorkWithGeometricRules,
3939 &bC6ParametersWorkWithLBRules,
3941 if (ir->ljpme_combination_rule == eljpmeLB)
3943 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3945 warning(wi, "You are using arithmetic-geometric combination rules "
3946 "in LJ-PME, but your non-bonded C6 parameters do not "
3947 "follow these rules.");
3952 if (FALSE == bC6ParametersWorkWithGeometricRules)
3954 if (ir->eDispCorr != edispcNO)
3956 warning_note(wi, "You are using geometric combination rules in "
3957 "LJ-PME, but your non-bonded C6 parameters do "
3958 "not follow these rules. "
3959 "This will introduce very small errors in the forces and energies in "
3960 "your simulations. Dispersion correction will correct total energy "
3961 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3965 warning_note(wi, "You are using geometric combination rules in "
3966 "LJ-PME, but your non-bonded C6 parameters do "
3967 "not follow these rules. "
3968 "This will introduce very small errors in the forces and energies in "
3969 "your simulations. If your system is homogeneous, consider using dispersion correction "
3970 "for the total energy and pressure.");
3976 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3979 char err_buf[STRLEN];
3980 int i, m, c, nmol, npct;
3981 gmx_bool bCharge, bAcc;
3982 real gdt_max, *mgrp, mt;
3984 gmx_mtop_atomloop_block_t aloopb;
3985 gmx_mtop_atomloop_all_t aloop;
3988 char warn_buf[STRLEN];
3990 set_warning_line(wi, mdparin, -1);
3992 if (ir->cutoff_scheme == ecutsVERLET &&
3993 ir->verletbuf_tol > 0 &&
3995 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
3996 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
3998 /* Check if a too small Verlet buffer might potentially
3999 * cause more drift than the thermostat can couple off.
4001 /* Temperature error fraction for warning and suggestion */
4002 const real T_error_warn = 0.002;
4003 const real T_error_suggest = 0.001;
4004 /* For safety: 2 DOF per atom (typical with constraints) */
4005 const real nrdf_at = 2;
4006 real T, tau, max_T_error;
4011 for (i = 0; i < ir->opts.ngtc; i++)
4013 T = max(T, ir->opts.ref_t[i]);
4014 tau = max(tau, ir->opts.tau_t[i]);
4018 /* This is a worst case estimate of the temperature error,
4019 * assuming perfect buffer estimation and no cancelation
4020 * of errors. The factor 0.5 is because energy distributes
4021 * equally over Ekin and Epot.
4023 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4024 if (max_T_error > T_error_warn)
4026 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.",
4027 ir->verletbuf_tol, T, tau,
4029 100*T_error_suggest,
4030 ir->verletbuf_tol*T_error_suggest/max_T_error);
4031 warning(wi, warn_buf);
4036 if (ETC_ANDERSEN(ir->etc))
4040 for (i = 0; i < ir->opts.ngtc; i++)
4042 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4043 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4044 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
4045 i, ir->opts.tau_t[i]);
4046 CHECK(ir->opts.tau_t[i] < 0);
4049 for (i = 0; i < ir->opts.ngtc; i++)
4051 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
4052 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);
4053 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
4057 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4058 ir->comm_mode == ecmNO &&
4059 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4060 !ETC_ANDERSEN(ir->etc))
4062 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");
4065 /* Check for pressure coupling with absolute position restraints */
4066 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4068 absolute_reference(ir, sys, TRUE, AbsRef);
4070 for (m = 0; m < DIM; m++)
4072 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4074 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4082 aloopb = gmx_mtop_atomloop_block_init(sys);
4083 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4085 if (atom->q != 0 || atom->qB != 0)
4093 if (EEL_FULL(ir->coulombtype))
4096 "You are using full electrostatics treatment %s for a system without charges.\n"
4097 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4098 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4099 warning(wi, err_buf);
4104 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4107 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4108 "You might want to consider using %s electrostatics.\n",
4110 warning_note(wi, err_buf);
4114 /* Check if combination rules used in LJ-PME are the same as in the force field */
4115 if (EVDW_PME(ir->vdwtype))
4117 check_combination_rules(ir, sys, wi);
4120 /* Generalized reaction field */
4121 if (ir->opts.ngtc == 0)
4123 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4125 CHECK(ir->coulombtype == eelGRF);
4129 sprintf(err_buf, "When using coulombtype = %s"
4130 " ref-t for temperature coupling should be > 0",
4132 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4135 if (ir->eI == eiSD2)
4137 sprintf(warn_buf, "The stochastic dynamics integrator %s is deprecated, since\n"
4138 "it is slower than integrator %s and is slightly less accurate\n"
4139 "with constraints. Use the %s integrator.",
4140 ei_names[ir->eI], ei_names[eiSD1], ei_names[eiSD1]);
4141 warning_note(wi, warn_buf);
4145 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4147 for (m = 0; (m < DIM); m++)
4149 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4158 snew(mgrp, sys->groups.grps[egcACC].nr);
4159 aloop = gmx_mtop_atomloop_all_init(sys);
4160 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4162 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4165 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4167 for (m = 0; (m < DIM); m++)
4169 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4173 for (m = 0; (m < DIM); m++)
4175 if (fabs(acc[m]) > 1e-6)
4177 const char *dim[DIM] = { "X", "Y", "Z" };
4179 "Net Acceleration in %s direction, will %s be corrected\n",
4180 dim[m], ir->nstcomm != 0 ? "" : "not");
4181 if (ir->nstcomm != 0 && m < ndof_com(ir))
4184 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4186 ir->opts.acc[i][m] -= acc[m];
4194 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4195 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4197 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4200 if (ir->ePull != epullNO)
4202 gmx_bool bPullAbsoluteRef;
4204 bPullAbsoluteRef = FALSE;
4205 for (i = 0; i < ir->pull->ncoord; i++)
4207 bPullAbsoluteRef = bPullAbsoluteRef ||
4208 ir->pull->coord[i].group[0] == 0 ||
4209 ir->pull->coord[i].group[1] == 0;
4211 if (bPullAbsoluteRef)
4213 absolute_reference(ir, sys, FALSE, AbsRef);
4214 for (m = 0; m < DIM; m++)
4216 if (ir->pull->dim[m] && !AbsRef[m])
4218 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.");
4224 if (ir->pull->eGeom == epullgDIRPBC)
4226 for (i = 0; i < 3; i++)
4228 for (m = 0; m <= i; m++)
4230 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4231 ir->deform[i][m] != 0)
4233 for (c = 0; c < ir->pull->ncoord; c++)
4235 if (ir->pull->coord[c].vec[m] != 0)
4237 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
4249 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
4253 char warn_buf[STRLEN];
4256 ptr = check_box(ir->ePBC, box);
4259 warning_error(wi, ptr);
4262 if (bConstr && ir->eConstrAlg == econtSHAKE)
4264 if (ir->shake_tol <= 0.0)
4266 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4268 warning_error(wi, warn_buf);
4271 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
4273 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
4274 if (ir->epc == epcNO)
4276 warning(wi, warn_buf);
4280 warning_error(wi, warn_buf);
4285 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
4287 /* If we have Lincs constraints: */
4288 if (ir->eI == eiMD && ir->etc == etcNO &&
4289 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4291 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4292 warning_note(wi, warn_buf);
4295 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4297 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4298 warning_note(wi, warn_buf);
4300 if (ir->epc == epcMTTK)
4302 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4306 if (bConstr && ir->epc == epcMTTK)
4308 warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
4311 if (ir->LincsWarnAngle > 90.0)
4313 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4314 warning(wi, warn_buf);
4315 ir->LincsWarnAngle = 90.0;
4318 if (ir->ePBC != epbcNONE)
4320 if (ir->nstlist == 0)
4322 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4324 bTWIN = (ir->rlistlong > ir->rlist);
4325 if (ir->ns_type == ensGRID)
4327 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
4329 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 %s.\n",
4330 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
4331 warning_error(wi, warn_buf);
4336 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
4337 if (2*ir->rlistlong >= min_size)
4339 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4340 warning_error(wi, warn_buf);
4343 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4350 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4354 real rvdw1, rvdw2, rcoul1, rcoul2;
4355 char warn_buf[STRLEN];
4357 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4361 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4366 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4372 if (rvdw1 + rvdw2 > ir->rlist ||
4373 rcoul1 + rcoul2 > ir->rlist)
4376 "The sum of the two largest charge group radii (%f) "
4377 "is larger than rlist (%f)\n",
4378 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4379 warning(wi, warn_buf);
4383 /* Here we do not use the zero at cut-off macro,
4384 * since user defined interactions might purposely
4385 * not be zero at the cut-off.
4387 if (ir_vdw_is_zero_at_cutoff(ir) &&
4388 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4390 sprintf(warn_buf, "The sum of the two largest charge group "
4391 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4392 "With exact cut-offs, better performance can be "
4393 "obtained with cutoff-scheme = %s, because it "
4394 "does not use charge groups at all.",
4396 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4397 ir->rlistlong, ir->rvdw,
4398 ecutscheme_names[ecutsVERLET]);
4401 warning(wi, warn_buf);
4405 warning_note(wi, warn_buf);
4408 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4409 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4411 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4412 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4414 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4415 ir->rlistlong, ir->rcoulomb,
4416 ecutscheme_names[ecutsVERLET]);
4419 warning(wi, warn_buf);
4423 warning_note(wi, warn_buf);