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.
45 #include "gromacs/utility/smalloc.h"
49 #include "gmx_fatal.h"
53 #include "gromacs/utility/cstringutil.h"
62 #include "mtop_util.h"
63 #include "chargegroup.h"
65 #include "calc_verletbuf.h"
70 /* Resource parameters
71 * Do not change any of these until you read the instruction
72 * in readinp.h. Some cpp's do not take spaces after the backslash
73 * (like the c-shell), which will give you a very weird compiler
77 typedef struct t_inputrec_strings
79 char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
80 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
81 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], x_compressed_groups[STRLEN],
82 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
83 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN],
85 char fep_lambda[efptNR][STRLEN];
86 char lambda_weights[STRLEN];
89 char anneal[STRLEN], anneal_npoints[STRLEN],
90 anneal_time[STRLEN], anneal_temp[STRLEN];
91 char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
92 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
93 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
94 char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
95 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
97 } gmx_inputrec_strings;
99 static gmx_inputrec_strings *is = NULL;
101 void init_inputrec_strings()
105 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.");
110 void done_inputrec_strings()
116 static char swapgrp[STRLEN], splitgrp0[STRLEN], splitgrp1[STRLEN], solgrp[STRLEN];
119 egrptpALL, /* All particles have to be a member of a group. */
120 egrptpALL_GENREST, /* A rest group with name is generated for particles *
121 * that are not part of any group. */
122 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
123 * for the rest group. */
124 egrptpONE /* Merge all selected groups into one group, *
125 * make a rest group for the remaining particles. */
128 static const char *constraints[eshNR+1] = {
129 "none", "h-bonds", "all-bonds", "h-angles", "all-angles", NULL
132 static const char *couple_lam[ecouplamNR+1] = {
133 "vdw-q", "vdw", "q", "none", NULL
136 void init_ir(t_inputrec *ir, t_gromppopts *opts)
138 snew(opts->include, STRLEN);
139 snew(opts->define, STRLEN);
140 snew(ir->fepvals, 1);
141 snew(ir->expandedvals, 1);
142 snew(ir->simtempvals, 1);
145 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
150 for (i = 0; i < ntemps; i++)
152 /* simple linear scaling -- allows more control */
153 if (simtemp->eSimTempScale == esimtempLINEAR)
155 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
157 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
159 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
161 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
163 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
168 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
169 gmx_fatal(FARGS, errorstr);
176 static void _low_check(gmx_bool b, char *s, warninp_t wi)
180 warning_error(wi, s);
184 static void check_nst(const char *desc_nst, int nst,
185 const char *desc_p, int *p,
190 if (*p > 0 && *p % nst != 0)
192 /* Round up to the next multiple of nst */
193 *p = ((*p)/nst + 1)*nst;
194 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
195 desc_p, desc_nst, desc_p, *p);
200 static gmx_bool ir_NVE(const t_inputrec *ir)
202 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
205 static int lcd(int n1, int n2)
210 for (i = 2; (i <= n1 && i <= n2); i++)
212 if (n1 % i == 0 && n2 % i == 0)
221 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
223 if (*eintmod == eintmodPOTSHIFT_VERLET)
225 if (ir->cutoff_scheme == ecutsVERLET)
227 *eintmod = eintmodPOTSHIFT;
231 *eintmod = eintmodNONE;
236 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
238 /* Check internal consistency.
239 * NOTE: index groups are not set here yet, don't check things
240 * like temperature coupling group options here, but in triple_check
243 /* Strange macro: first one fills the err_buf, and then one can check
244 * the condition, which will print the message and increase the error
247 #define CHECK(b) _low_check(b, err_buf, wi)
248 char err_buf[256], warn_buf[STRLEN];
254 t_lambda *fep = ir->fepvals;
255 t_expanded *expand = ir->expandedvals;
257 set_warning_line(wi, mdparin, -1);
259 /* BASIC CUT-OFF STUFF */
260 if (ir->rcoulomb < 0)
262 warning_error(wi, "rcoulomb should be >= 0");
266 warning_error(wi, "rvdw should be >= 0");
269 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0))
271 warning_error(wi, "rlist should be >= 0");
274 process_interaction_modifier(ir, &ir->coulomb_modifier);
275 process_interaction_modifier(ir, &ir->vdw_modifier);
277 if (ir->cutoff_scheme == ecutsGROUP)
280 "The group cutoff scheme is deprecated in Gromacs 5.0 and will be removed in a future "
281 "release when all interaction forms are supported for the verlet scheme. The verlet "
282 "scheme already scales better, and it is compatible with GPUs and other accelerators.");
284 /* BASIC CUT-OFF STUFF */
285 if (ir->rlist == 0 ||
286 !((ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > ir->rlist) ||
287 (ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > ir->rlist)))
289 /* No switched potential and/or no twin-range:
290 * we can set the long-range cut-off to the maximum of the other cut-offs.
292 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
294 else if (ir->rlistlong < 0)
296 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
297 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
299 warning(wi, warn_buf);
301 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
303 warning_error(wi, "Can not have an infinite cut-off with PBC");
305 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
307 warning_error(wi, "rlistlong can not be shorter than rlist");
309 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
311 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
315 if (ir->rlistlong == ir->rlist)
319 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
321 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
324 if (ir->cutoff_scheme == ecutsVERLET)
328 /* Normal Verlet type neighbor-list, currently only limited feature support */
329 if (inputrec2nboundeddim(ir) < 3)
331 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
333 if (ir->rcoulomb != ir->rvdw)
335 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
337 if (ir->vdwtype == evdwSHIFT || ir->vdwtype == evdwSWITCH)
339 if (ir->vdw_modifier == eintmodNONE ||
340 ir->vdw_modifier == eintmodPOTSHIFT)
342 ir->vdw_modifier = (ir->vdwtype == evdwSHIFT ? eintmodFORCESWITCH : eintmodPOTSWITCH);
344 sprintf(warn_buf, "Replacing vdwtype=%s by the equivalent combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], evdw_names[evdwCUT], eintmod_names[ir->vdw_modifier]);
345 warning_note(wi, warn_buf);
347 ir->vdwtype = evdwCUT;
351 sprintf(warn_buf, "Unsupported combination of vdwtype=%s and vdw_modifier=%s", evdw_names[ir->vdwtype], eintmod_names[ir->vdw_modifier]);
352 warning_error(wi, warn_buf);
356 if (!(ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME))
358 warning_error(wi, "With Verlet lists only cut-off and PME LJ interactions are supported");
360 if (!(ir->coulombtype == eelCUT ||
361 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
362 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
364 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
366 if (!(ir->coulomb_modifier == eintmodNONE ||
367 ir->coulomb_modifier == eintmodPOTSHIFT))
369 sprintf(warn_buf, "coulomb_modifier=%s is not supported with the Verlet cut-off scheme", eintmod_names[ir->coulomb_modifier]);
370 warning_error(wi, warn_buf);
373 if (ir->nstlist <= 0)
375 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
378 if (ir->nstlist < 10)
380 warning_note(wi, "With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note that with the Verlet scheme, nstlist has no effect on the accuracy of your simulation.");
383 rc_max = max(ir->rvdw, ir->rcoulomb);
385 if (ir->verletbuf_tol <= 0)
387 if (ir->verletbuf_tol == 0)
389 warning_error(wi, "Can not have Verlet buffer tolerance of exactly 0");
392 if (ir->rlist < rc_max)
394 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
397 if (ir->rlist == rc_max && ir->nstlist > 1)
399 warning_note(wi, "rlist is equal to rvdw and/or rcoulomb: there is no explicit Verlet buffer. The cluster pair list does have a buffering effect, but choosing a larger rlist might be necessary for good energy conservation.");
404 if (ir->rlist > rc_max)
406 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-tolerance > 0. Will set rlist using verlet-buffer-tolerance.");
409 if (ir->nstlist == 1)
411 /* No buffer required */
416 if (EI_DYNAMICS(ir->eI))
418 if (inputrec2nboundeddim(ir) < 3)
420 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-tolerance > 0. You are using at least one unbounded dimension, so no volume can be computed. Either use a finite box, or set rlist yourself together with verlet-buffer-tolerance = -1.");
422 /* Set rlist temporarily so we can continue processing */
427 /* Set the buffer to 5% of the cut-off */
428 ir->rlist = (1.0 + verlet_buffer_ratio_nodynamics)*rc_max;
433 /* No twin-range calculations with Verlet lists */
434 ir->rlistlong = ir->rlist;
437 if (ir->nstcalclr == -1)
439 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
440 ir->nstcalclr = ir->nstlist;
442 else if (ir->nstcalclr > 0)
444 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
446 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
449 else if (ir->nstcalclr < -1)
451 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
454 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
456 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
459 /* GENERAL INTEGRATOR STUFF */
460 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
464 if (ir->eI == eiVVAK)
466 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]);
467 warning_note(wi, warn_buf);
469 if (!EI_DYNAMICS(ir->eI))
473 if (EI_DYNAMICS(ir->eI))
475 if (ir->nstcalcenergy < 0)
477 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
478 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
480 /* nstcalcenergy larger than nstener does not make sense.
481 * We ideally want nstcalcenergy=nstener.
485 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
489 ir->nstcalcenergy = ir->nstenergy;
493 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
494 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
495 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
498 const char *nsten = "nstenergy";
499 const char *nstdh = "nstdhdl";
500 const char *min_name = nsten;
501 int min_nst = ir->nstenergy;
503 /* find the smallest of ( nstenergy, nstdhdl ) */
504 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
505 (ir->nstenergy == 0 || ir->fepvals->nstdhdl < ir->nstenergy))
507 min_nst = ir->fepvals->nstdhdl;
510 /* If the user sets nstenergy small, we should respect that */
512 "Setting nstcalcenergy (%d) equal to %s (%d)",
513 ir->nstcalcenergy, min_name, min_nst);
514 warning_note(wi, warn_buf);
515 ir->nstcalcenergy = min_nst;
518 if (ir->epc != epcNO)
520 if (ir->nstpcouple < 0)
522 ir->nstpcouple = ir_optimal_nstpcouple(ir);
525 if (IR_TWINRANGE(*ir))
527 check_nst("nstlist", ir->nstlist,
528 "nstcalcenergy", &ir->nstcalcenergy, wi);
529 if (ir->epc != epcNO)
531 check_nst("nstlist", ir->nstlist,
532 "nstpcouple", &ir->nstpcouple, wi);
536 if (ir->nstcalcenergy > 0)
538 if (ir->efep != efepNO)
540 /* nstdhdl should be a multiple of nstcalcenergy */
541 check_nst("nstcalcenergy", ir->nstcalcenergy,
542 "nstdhdl", &ir->fepvals->nstdhdl, wi);
543 /* nstexpanded should be a multiple of nstcalcenergy */
544 check_nst("nstcalcenergy", ir->nstcalcenergy,
545 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
547 /* for storing exact averages nstenergy should be
548 * a multiple of nstcalcenergy
550 check_nst("nstcalcenergy", ir->nstcalcenergy,
551 "nstenergy", &ir->nstenergy, wi);
555 if (ir->nsteps == 0 && !ir->bContinuation)
557 warning_note(wi, "For a correct single-point energy evaluation with nsteps = 0, use continuation = yes to avoid constraining the input coordinates.");
561 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
562 ir->bContinuation && ir->ld_seed != -1)
564 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)");
570 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
571 CHECK(ir->ePBC != epbcXYZ);
572 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
573 CHECK(ir->ns_type != ensGRID);
574 sprintf(err_buf, "with TPI nstlist should be larger than zero");
575 CHECK(ir->nstlist <= 0);
576 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
577 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
578 sprintf(err_buf, "TPI does not work (yet) with the Verlet cut-off scheme");
579 CHECK(ir->cutoff_scheme == ecutsVERLET);
583 if ( (opts->nshake > 0) && (opts->bMorse) )
586 "Using morse bond-potentials while constraining bonds is useless");
587 warning(wi, warn_buf);
590 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
591 ir->bContinuation && ir->ld_seed != -1)
593 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)");
595 /* verify simulated tempering options */
599 gmx_bool bAllTempZero = TRUE;
600 for (i = 0; i < fep->n_lambda; i++)
602 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]);
603 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
604 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
606 bAllTempZero = FALSE;
609 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
610 CHECK(bAllTempZero == TRUE);
612 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
613 CHECK(ir->eI != eiVV);
615 /* check compatability of the temperature coupling with simulated tempering */
617 if (ir->etc == etcNOSEHOOVER)
619 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
620 warning_note(wi, warn_buf);
623 /* check that the temperatures make sense */
625 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);
626 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
628 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
629 CHECK(ir->simtempvals->simtemp_high <= 0);
631 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
632 CHECK(ir->simtempvals->simtemp_low <= 0);
635 /* verify free energy options */
637 if (ir->efep != efepNO)
640 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
642 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
644 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
645 (int)fep->sc_r_power);
646 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
648 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
649 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
651 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
652 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
654 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
655 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
657 sprintf(err_buf, "Free-energy not implemented for Ewald");
658 CHECK(ir->coulombtype == eelEWALD);
660 /* check validty of lambda inputs */
661 if (fep->n_lambda == 0)
663 /* Clear output in case of no states:*/
664 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
665 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
669 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
670 CHECK((fep->init_fep_state >= fep->n_lambda));
673 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
674 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
676 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",
677 fep->init_lambda, fep->init_fep_state);
678 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
682 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
686 for (i = 0; i < efptNR; i++)
688 if (fep->separate_dvdl[i])
693 if (n_lambda_terms > 1)
695 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.");
696 warning(wi, warn_buf);
699 if (n_lambda_terms < 2 && fep->n_lambda > 0)
702 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
706 for (j = 0; j < efptNR; j++)
708 for (i = 0; i < fep->n_lambda; i++)
710 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]);
711 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
715 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
717 for (i = 0; i < fep->n_lambda; i++)
719 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],
720 fep->all_lambda[efptCOUL][i]);
721 CHECK((fep->sc_alpha > 0) &&
722 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
723 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
724 ((fep->all_lambda[efptVDW][i] > 0.0) &&
725 (fep->all_lambda[efptVDW][i] < 1.0))));
729 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
731 real sigma, lambda, r_sc;
734 /* Maximum estimate for A and B charges equal with lambda power 1 */
736 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
737 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.",
739 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
740 warning_note(wi, warn_buf);
743 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
744 be treated differently, but that's the next step */
746 for (i = 0; i < efptNR; i++)
748 for (j = 0; j < fep->n_lambda; j++)
750 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
751 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
756 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
759 expand = ir->expandedvals;
761 /* checking equilibration of weights inputs for validity */
763 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
764 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
765 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
767 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
768 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
769 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
771 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
772 expand->equil_steps, elmceq_names[elmceqSTEPS]);
773 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
775 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
776 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
777 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
779 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
780 expand->equil_ratio, elmceq_names[elmceqRATIO]);
781 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
783 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
784 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
785 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
787 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
788 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
789 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
791 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
792 expand->equil_steps, elmceq_names[elmceqSTEPS]);
793 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
795 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
796 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
797 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
799 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
800 expand->equil_ratio, elmceq_names[elmceqRATIO]);
801 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
803 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
804 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
805 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
807 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
808 CHECK((expand->lmc_repeats <= 0));
809 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
810 CHECK((expand->minvarmin <= 0));
811 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
812 CHECK((expand->c_range < 0));
813 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
814 fep->init_fep_state, expand->lmc_forced_nstart);
815 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
816 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
817 CHECK((expand->lmc_forced_nstart < 0));
818 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
819 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
821 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
822 CHECK((expand->init_wl_delta < 0));
823 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
824 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
825 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
826 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
828 /* if there is no temperature control, we need to specify an MC temperature */
829 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
830 if (expand->nstTij > 0)
832 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
833 expand->nstTij, ir->nstlog);
834 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
839 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
840 CHECK(ir->nwall && ir->ePBC != epbcXY);
843 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
845 if (ir->ePBC == epbcNONE)
847 if (ir->epc != epcNO)
849 warning(wi, "Turning off pressure coupling for vacuum system");
855 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
856 epbc_names[ir->ePBC]);
857 CHECK(ir->epc != epcNO);
859 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
860 CHECK(EEL_FULL(ir->coulombtype));
862 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
863 epbc_names[ir->ePBC]);
864 CHECK(ir->eDispCorr != edispcNO);
867 if (ir->rlist == 0.0)
869 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
870 "with coulombtype = %s or coulombtype = %s\n"
871 "without periodic boundary conditions (pbc = %s) and\n"
872 "rcoulomb and rvdw set to zero",
873 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
874 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
875 (ir->ePBC != epbcNONE) ||
876 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
880 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
884 warning_note(wi, "Simulating without cut-offs can be (slightly) faster with nstlist=0, nstype=simple and only one MPI rank");
889 if (ir->nstcomm == 0)
891 ir->comm_mode = ecmNO;
893 if (ir->comm_mode != ecmNO)
897 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");
898 ir->nstcomm = abs(ir->nstcomm);
901 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
903 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
904 ir->nstcomm = ir->nstcalcenergy;
907 if (ir->comm_mode == ecmANGULAR)
909 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
910 CHECK(ir->bPeriodicMols);
911 if (ir->ePBC != epbcNONE)
913 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).");
918 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
920 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.");
923 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
924 " algorithm not implemented");
925 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
926 && (ir->ns_type == ensSIMPLE));
928 /* TEMPERATURE COUPLING */
929 if (ir->etc == etcYES)
931 ir->etc = etcBERENDSEN;
932 warning_note(wi, "Old option for temperature coupling given: "
933 "changing \"yes\" to \"Berendsen\"\n");
936 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
938 if (ir->opts.nhchainlength < 1)
940 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
941 ir->opts.nhchainlength = 1;
942 warning(wi, warn_buf);
945 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
947 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
948 ir->opts.nhchainlength = 1;
953 ir->opts.nhchainlength = 0;
956 if (ir->eI == eiVVAK)
958 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
960 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
963 if (ETC_ANDERSEN(ir->etc))
965 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
966 CHECK(!(EI_VV(ir->eI)));
968 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
970 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]);
971 warning_note(wi, warn_buf);
974 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]);
975 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
978 if (ir->etc == etcBERENDSEN)
980 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
981 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
982 warning_note(wi, warn_buf);
985 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
986 && ir->epc == epcBERENDSEN)
988 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
989 "true ensemble for the thermostat");
990 warning(wi, warn_buf);
993 /* PRESSURE COUPLING */
994 if (ir->epc == epcISOTROPIC)
996 ir->epc = epcBERENDSEN;
997 warning_note(wi, "Old option for pressure coupling given: "
998 "changing \"Isotropic\" to \"Berendsen\"\n");
1001 if (ir->epc != epcNO)
1003 dt_pcoupl = ir->nstpcouple*ir->delta_t;
1005 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
1006 CHECK(ir->tau_p <= 0);
1008 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc) - 10*GMX_REAL_EPS)
1010 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
1011 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
1012 warning(wi, warn_buf);
1015 sprintf(err_buf, "compressibility must be > 0 when using pressure"
1016 " coupling %s\n", EPCOUPLTYPE(ir->epc));
1017 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
1018 ir->compress[ZZ][ZZ] < 0 ||
1019 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
1020 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
1022 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
1025 "You are generating velocities so I am assuming you "
1026 "are equilibrating a system. You are using "
1027 "%s pressure coupling, but this can be "
1028 "unstable for equilibration. If your system crashes, try "
1029 "equilibrating first with Berendsen pressure coupling. If "
1030 "you are not equilibrating the system, you can probably "
1031 "ignore this warning.",
1032 epcoupl_names[ir->epc]);
1033 warning(wi, warn_buf);
1039 if (ir->epc > epcNO)
1041 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
1043 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.");
1049 if (ir->epc == epcMTTK)
1051 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1055 /* ELECTROSTATICS */
1056 /* More checks are in triple check (grompp.c) */
1058 if (ir->coulombtype == eelSWITCH)
1060 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1061 "artifacts, advice: use coulombtype = %s",
1062 eel_names[ir->coulombtype],
1063 eel_names[eelRF_ZERO]);
1064 warning(wi, warn_buf);
1067 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1069 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1070 warning_note(wi, warn_buf);
1073 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1075 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);
1076 warning(wi, warn_buf);
1077 ir->epsilon_rf = ir->epsilon_r;
1078 ir->epsilon_r = 1.0;
1081 if (getenv("GMX_DO_GALACTIC_DYNAMICS") == NULL)
1083 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1084 CHECK(ir->epsilon_r < 0);
1087 if (EEL_RF(ir->coulombtype))
1089 /* reaction field (at the cut-off) */
1091 if (ir->coulombtype == eelRF_ZERO)
1093 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1094 eel_names[ir->coulombtype]);
1095 CHECK(ir->epsilon_rf != 0);
1096 ir->epsilon_rf = 0.0;
1099 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1100 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1101 (ir->epsilon_r == 0));
1102 if (ir->epsilon_rf == ir->epsilon_r)
1104 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1105 eel_names[ir->coulombtype]);
1106 warning(wi, warn_buf);
1109 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1110 * means the interaction is zero outside rcoulomb, but it helps to
1111 * provide accurate energy conservation.
1113 if (ir_coulomb_might_be_zero_at_cutoff(ir))
1115 if (ir_coulomb_switched(ir))
1118 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1119 eel_names[ir->coulombtype]);
1120 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1123 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1125 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1127 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1128 eel_names[ir->coulombtype]);
1129 CHECK(ir->rlist > ir->rcoulomb);
1133 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT)
1136 "Explicit switch/shift coulomb interactions cannot be used in combination with a secondary coulomb-modifier.");
1137 CHECK( ir->coulomb_modifier != eintmodNONE);
1139 if (ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1142 "Explicit switch/shift vdw interactions cannot be used in combination with a secondary vdw-modifier.");
1143 CHECK( ir->vdw_modifier != eintmodNONE);
1146 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1147 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1150 "The switch/shift interaction settings are just for compatibility; you will get better "
1151 "performance from applying potential modifiers to your interactions!\n");
1152 warning_note(wi, warn_buf);
1155 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1157 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1159 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1160 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.",
1161 percentage, ir->rcoulomb_switch, ir->rcoulomb, ir->ewald_rtol);
1162 warning(wi, warn_buf);
1166 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1168 if (ir->rvdw_switch == 0)
1170 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.");
1171 warning(wi, warn_buf);
1175 if (EEL_FULL(ir->coulombtype))
1177 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1178 ir->coulombtype == eelPMEUSERSWITCH)
1180 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1181 eel_names[ir->coulombtype]);
1182 CHECK(ir->rcoulomb > ir->rlist);
1184 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1186 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1189 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1190 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1191 "a potential modifier.", eel_names[ir->coulombtype]);
1192 if (ir->nstcalclr == 1)
1194 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1198 CHECK(ir->rcoulomb != ir->rlist);
1204 if (EEL_PME(ir->coulombtype) || EVDW_PME(ir->vdwtype))
1206 if (ir->pme_order < 3)
1208 warning_error(wi, "pme-order can not be smaller than 3");
1212 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1214 if (ir->ewald_geometry == eewg3D)
1216 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1217 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1218 warning(wi, warn_buf);
1220 /* This check avoids extra pbc coding for exclusion corrections */
1221 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1222 CHECK(ir->wall_ewald_zfac < 2);
1225 if (ir_vdw_switched(ir))
1227 sprintf(err_buf, "With switched vdw forces or potentials, rvdw-switch must be < rvdw");
1228 CHECK(ir->rvdw_switch >= ir->rvdw);
1230 if (ir->rvdw_switch < 0.5*ir->rvdw)
1232 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.",
1233 ir->rvdw_switch, ir->rvdw);
1234 warning_note(wi, warn_buf);
1237 else if (ir->vdwtype == evdwCUT || ir->vdwtype == evdwPME)
1239 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1241 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1242 CHECK(ir->rlist > ir->rvdw);
1246 if (ir->vdwtype == evdwPME)
1248 if (!(ir->vdw_modifier == eintmodNONE || ir->vdw_modifier == eintmodPOTSHIFT))
1250 sprintf(err_buf, "With vdwtype = %s, the only supported modifiers are %s a\
1252 evdw_names[ir->vdwtype],
1253 eintmod_names[eintmodPOTSHIFT],
1254 eintmod_names[eintmodNONE]);
1258 if (ir->cutoff_scheme == ecutsGROUP)
1260 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1261 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1264 warning_note(wi, "With exact cut-offs, rlist should be "
1265 "larger than rcoulomb and rvdw, so that there "
1266 "is a buffer region for particle motion "
1267 "between neighborsearch steps");
1270 if (ir_coulomb_is_zero_at_cutoff(ir) && ir->rlistlong <= ir->rcoulomb)
1272 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1273 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1274 warning_note(wi, warn_buf);
1276 if (ir_vdw_switched(ir) && (ir->rlistlong <= ir->rvdw))
1278 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1279 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1280 warning_note(wi, warn_buf);
1284 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1286 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.");
1289 if (ir->nstlist == -1)
1291 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1292 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1294 sprintf(err_buf, "nstlist can not be smaller than -1");
1295 CHECK(ir->nstlist < -1);
1297 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1300 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1303 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1305 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1308 /* ENERGY CONSERVATION */
1309 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1311 if (!ir_vdw_might_be_zero_at_cutoff(ir) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1313 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1314 evdw_names[evdwSHIFT]);
1315 warning_note(wi, warn_buf);
1317 if (!ir_coulomb_might_be_zero_at_cutoff(ir) && ir->rcoulomb > 0)
1319 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1320 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1321 warning_note(wi, warn_buf);
1325 if (EI_VV(ir->eI) && IR_TWINRANGE(*ir) && ir->nstlist > 1)
1327 sprintf(warn_buf, "Twin-range multiple time stepping does not work with integrator %s.", ei_names[ir->eI]);
1328 warning_error(wi, warn_buf);
1331 /* IMPLICIT SOLVENT */
1332 if (ir->coulombtype == eelGB_NOTUSED)
1334 ir->coulombtype = eelCUT;
1335 ir->implicit_solvent = eisGBSA;
1336 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1337 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1338 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1341 if (ir->sa_algorithm == esaSTILL)
1343 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1344 CHECK(ir->sa_algorithm == esaSTILL);
1347 if (ir->implicit_solvent == eisGBSA)
1349 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1350 CHECK(ir->rgbradii != ir->rlist);
1352 if (ir->coulombtype != eelCUT)
1354 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1355 CHECK(ir->coulombtype != eelCUT);
1357 if (ir->vdwtype != evdwCUT)
1359 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1360 CHECK(ir->vdwtype != evdwCUT);
1362 if (ir->nstgbradii < 1)
1364 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1365 warning_note(wi, warn_buf);
1368 if (ir->sa_algorithm == esaNO)
1370 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1371 warning_note(wi, warn_buf);
1373 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1375 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");
1376 warning_note(wi, warn_buf);
1378 if (ir->gb_algorithm == egbSTILL)
1380 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1384 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1387 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1389 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1390 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1397 if (ir->cutoff_scheme != ecutsGROUP)
1399 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1403 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1405 if (ir->epc != epcNO)
1407 warning_error(wi, "AdresS simulation does not support pressure coupling");
1409 if (EEL_FULL(ir->coulombtype))
1411 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1416 /* count the number of text elemets separated by whitespace in a string.
1417 str = the input string
1418 maxptr = the maximum number of allowed elements
1419 ptr = the output array of pointers to the first character of each element
1420 returns: the number of elements. */
1421 int str_nelem(const char *str, int maxptr, char *ptr[])
1426 copy0 = strdup(str);
1429 while (*copy != '\0')
1433 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1441 while ((*copy != '\0') && !isspace(*copy))
1460 /* interpret a number of doubles from a string and put them in an array,
1461 after allocating space for them.
1462 str = the input string
1463 n = the (pre-allocated) number of doubles read
1464 r = the output array of doubles. */
1465 static void parse_n_real(char *str, int *n, real **r)
1470 *n = str_nelem(str, MAXPTR, ptr);
1473 for (i = 0; i < *n; i++)
1475 (*r)[i] = strtod(ptr[i], NULL);
1479 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1482 int i, j, max_n_lambda, nweights, nfep[efptNR];
1483 t_lambda *fep = ir->fepvals;
1484 t_expanded *expand = ir->expandedvals;
1485 real **count_fep_lambdas;
1486 gmx_bool bOneLambda = TRUE;
1488 snew(count_fep_lambdas, efptNR);
1490 /* FEP input processing */
1491 /* first, identify the number of lambda values for each type.
1492 All that are nonzero must have the same number */
1494 for (i = 0; i < efptNR; i++)
1496 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1499 /* now, determine the number of components. All must be either zero, or equal. */
1502 for (i = 0; i < efptNR; i++)
1504 if (nfep[i] > max_n_lambda)
1506 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1507 must have the same number if its not zero.*/
1512 for (i = 0; i < efptNR; i++)
1516 ir->fepvals->separate_dvdl[i] = FALSE;
1518 else if (nfep[i] == max_n_lambda)
1520 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1521 respect to the temperature currently */
1523 ir->fepvals->separate_dvdl[i] = TRUE;
1528 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1529 nfep[i], efpt_names[i], max_n_lambda);
1532 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1533 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1535 /* the number of lambdas is the number we've read in, which is either zero
1536 or the same for all */
1537 fep->n_lambda = max_n_lambda;
1539 /* allocate space for the array of lambda values */
1540 snew(fep->all_lambda, efptNR);
1541 /* if init_lambda is defined, we need to set lambda */
1542 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1544 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1546 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1547 for (i = 0; i < efptNR; i++)
1549 snew(fep->all_lambda[i], fep->n_lambda);
1550 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1553 for (j = 0; j < fep->n_lambda; j++)
1555 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1557 sfree(count_fep_lambdas[i]);
1560 sfree(count_fep_lambdas);
1562 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1563 bookkeeping -- for now, init_lambda */
1565 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1567 for (i = 0; i < fep->n_lambda; i++)
1569 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1573 /* check to see if only a single component lambda is defined, and soft core is defined.
1574 In this case, turn on coulomb soft core */
1576 if (max_n_lambda == 0)
1582 for (i = 0; i < efptNR; i++)
1584 if ((nfep[i] != 0) && (i != efptFEP))
1590 if ((bOneLambda) && (fep->sc_alpha > 0))
1592 fep->bScCoul = TRUE;
1595 /* Fill in the others with the efptFEP if they are not explicitly
1596 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1597 they are all zero. */
1599 for (i = 0; i < efptNR; i++)
1601 if ((nfep[i] == 0) && (i != efptFEP))
1603 for (j = 0; j < fep->n_lambda; j++)
1605 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1611 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1612 if (fep->sc_r_power == 48)
1614 if (fep->sc_alpha > 0.1)
1616 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1620 expand = ir->expandedvals;
1621 /* now read in the weights */
1622 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1625 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1627 else if (nweights != fep->n_lambda)
1629 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1630 nweights, fep->n_lambda);
1632 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1634 expand->nstexpanded = fep->nstdhdl;
1635 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1637 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1639 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1640 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1641 2*tau_t just to be careful so it's not to frequent */
1646 static void do_simtemp_params(t_inputrec *ir)
1649 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1650 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1655 static void do_wall_params(t_inputrec *ir,
1656 char *wall_atomtype, char *wall_density,
1660 char *names[MAXPTR];
1663 opts->wall_atomtype[0] = NULL;
1664 opts->wall_atomtype[1] = NULL;
1666 ir->wall_atomtype[0] = -1;
1667 ir->wall_atomtype[1] = -1;
1668 ir->wall_density[0] = 0;
1669 ir->wall_density[1] = 0;
1673 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1674 if (nstr != ir->nwall)
1676 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1679 for (i = 0; i < ir->nwall; i++)
1681 opts->wall_atomtype[i] = strdup(names[i]);
1684 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1686 nstr = str_nelem(wall_density, MAXPTR, names);
1687 if (nstr != ir->nwall)
1689 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1691 for (i = 0; i < ir->nwall; i++)
1693 sscanf(names[i], "%lf", &dbl);
1696 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1698 ir->wall_density[i] = dbl;
1704 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1712 srenew(groups->grpname, groups->ngrpname+nwall);
1713 grps = &(groups->grps[egcENER]);
1714 srenew(grps->nm_ind, grps->nr+nwall);
1715 for (i = 0; i < nwall; i++)
1717 sprintf(str, "wall%d", i);
1718 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1719 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1724 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1725 t_expanded *expand, warninp_t wi)
1727 int ninp, nerror = 0;
1733 /* read expanded ensemble parameters */
1734 CCTYPE ("expanded ensemble variables");
1735 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1736 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1737 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1738 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1739 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1740 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1741 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1742 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1743 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1744 CCTYPE("Seed for Monte Carlo in lambda space");
1745 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1746 RTYPE ("mc-temperature", expand->mc_temp, -1);
1747 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1748 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1749 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1750 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1751 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1752 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1753 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1754 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1755 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1756 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1757 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1765 void get_ir(const char *mdparin, const char *mdparout,
1766 t_inputrec *ir, t_gromppopts *opts,
1770 double dumdub[2][6];
1774 char warn_buf[STRLEN];
1775 t_lambda *fep = ir->fepvals;
1776 t_expanded *expand = ir->expandedvals;
1778 init_inputrec_strings();
1779 inp = read_inpfile(mdparin, &ninp, wi);
1781 snew(dumstr[0], STRLEN);
1782 snew(dumstr[1], STRLEN);
1784 if (-1 == search_einp(ninp, inp, "cutoff-scheme"))
1787 "%s did not specify a value for the .mdp option "
1788 "\"cutoff-scheme\". Probably it was first intended for use "
1789 "with GROMACS before 4.6. In 4.6, the Verlet scheme was "
1790 "introduced, but the group scheme was still the default. "
1791 "The default is now the Verlet scheme, so you will observe "
1792 "different behaviour.", mdparin);
1793 warning_note(wi, warn_buf);
1796 /* ignore the following deprecated commands */
1799 REM_TYPE("domain-decomposition");
1800 REM_TYPE("andersen-seed");
1802 REM_TYPE("dihre-fc");
1803 REM_TYPE("dihre-tau");
1804 REM_TYPE("nstdihreout");
1805 REM_TYPE("nstcheckpoint");
1806 REM_TYPE("optimize-fft");
1808 /* replace the following commands with the clearer new versions*/
1809 REPL_TYPE("unconstrained-start", "continuation");
1810 REPL_TYPE("foreign-lambda", "fep-lambdas");
1811 REPL_TYPE("verlet-buffer-drift", "verlet-buffer-tolerance");
1812 REPL_TYPE("nstxtcout", "nstxout-compressed");
1813 REPL_TYPE("xtc-grps", "compressed-x-grps");
1814 REPL_TYPE("xtc-precision", "compressed-x-precision");
1816 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1817 CTYPE ("Preprocessor information: use cpp syntax.");
1818 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1819 STYPE ("include", opts->include, NULL);
1820 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1821 STYPE ("define", opts->define, NULL);
1823 CCTYPE ("RUN CONTROL PARAMETERS");
1824 EETYPE("integrator", ir->eI, ei_names);
1825 CTYPE ("Start time and timestep in ps");
1826 RTYPE ("tinit", ir->init_t, 0.0);
1827 RTYPE ("dt", ir->delta_t, 0.001);
1828 STEPTYPE ("nsteps", ir->nsteps, 0);
1829 CTYPE ("For exact run continuation or redoing part of a run");
1830 STEPTYPE ("init-step", ir->init_step, 0);
1831 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1832 ITYPE ("simulation-part", ir->simulation_part, 1);
1833 CTYPE ("mode for center of mass motion removal");
1834 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1835 CTYPE ("number of steps for center of mass motion removal");
1836 ITYPE ("nstcomm", ir->nstcomm, 100);
1837 CTYPE ("group(s) for center of mass motion removal");
1838 STYPE ("comm-grps", is->vcm, NULL);
1840 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1841 CTYPE ("Friction coefficient (amu/ps) and random seed");
1842 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1843 STEPTYPE ("ld-seed", ir->ld_seed, -1);
1846 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1847 CTYPE ("Force tolerance and initial step-size");
1848 RTYPE ("emtol", ir->em_tol, 10.0);
1849 RTYPE ("emstep", ir->em_stepsize, 0.01);
1850 CTYPE ("Max number of iterations in relax-shells");
1851 ITYPE ("niter", ir->niter, 20);
1852 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1853 RTYPE ("fcstep", ir->fc_stepsize, 0);
1854 CTYPE ("Frequency of steepest descents steps when doing CG");
1855 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1856 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1858 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1859 RTYPE ("rtpi", ir->rtpi, 0.05);
1861 /* Output options */
1862 CCTYPE ("OUTPUT CONTROL OPTIONS");
1863 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1864 ITYPE ("nstxout", ir->nstxout, 0);
1865 ITYPE ("nstvout", ir->nstvout, 0);
1866 ITYPE ("nstfout", ir->nstfout, 0);
1867 CTYPE ("Output frequency for energies to log file and energy file");
1868 ITYPE ("nstlog", ir->nstlog, 1000);
1869 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1870 ITYPE ("nstenergy", ir->nstenergy, 1000);
1871 CTYPE ("Output frequency and precision for .xtc file");
1872 ITYPE ("nstxout-compressed", ir->nstxout_compressed, 0);
1873 RTYPE ("compressed-x-precision", ir->x_compression_precision, 1000.0);
1874 CTYPE ("This selects the subset of atoms for the compressed");
1875 CTYPE ("trajectory file. You can select multiple groups. By");
1876 CTYPE ("default, all atoms will be written.");
1877 STYPE ("compressed-x-grps", is->x_compressed_groups, NULL);
1878 CTYPE ("Selection of energy groups");
1879 STYPE ("energygrps", is->energy, NULL);
1881 /* Neighbor searching */
1882 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1883 CTYPE ("cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)");
1884 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1885 CTYPE ("nblist update frequency");
1886 ITYPE ("nstlist", ir->nstlist, 10);
1887 CTYPE ("ns algorithm (simple or grid)");
1888 EETYPE("ns-type", ir->ns_type, ens_names);
1889 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1890 EETYPE("pbc", ir->ePBC, epbc_names);
1891 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1892 CTYPE ("Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,");
1893 CTYPE ("a value of -1 means: use rlist");
1894 RTYPE("verlet-buffer-tolerance", ir->verletbuf_tol, 0.005);
1895 CTYPE ("nblist cut-off");
1896 RTYPE ("rlist", ir->rlist, 1.0);
1897 CTYPE ("long-range cut-off for switched potentials");
1898 RTYPE ("rlistlong", ir->rlistlong, -1);
1899 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1901 /* Electrostatics */
1902 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1903 CTYPE ("Method for doing electrostatics");
1904 EETYPE("coulombtype", ir->coulombtype, eel_names);
1905 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1906 CTYPE ("cut-off lengths");
1907 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1908 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1909 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1910 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1911 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1912 CTYPE ("Method for doing Van der Waals");
1913 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1914 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1915 CTYPE ("cut-off lengths");
1916 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1917 RTYPE ("rvdw", ir->rvdw, 1.0);
1918 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1919 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1920 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1921 RTYPE ("table-extension", ir->tabext, 1.0);
1922 CTYPE ("Separate tables between energy group pairs");
1923 STYPE ("energygrp-table", is->egptable, NULL);
1924 CTYPE ("Spacing for the PME/PPPM FFT grid");
1925 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1926 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1927 ITYPE ("fourier-nx", ir->nkx, 0);
1928 ITYPE ("fourier-ny", ir->nky, 0);
1929 ITYPE ("fourier-nz", ir->nkz, 0);
1930 CTYPE ("EWALD/PME/PPPM parameters");
1931 ITYPE ("pme-order", ir->pme_order, 4);
1932 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1933 RTYPE ("ewald-rtol-lj", ir->ewald_rtol_lj, 0.001);
1934 EETYPE("lj-pme-comb-rule", ir->ljpme_combination_rule, eljpme_names);
1935 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1936 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1938 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1939 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1941 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1942 CTYPE ("Algorithm for calculating Born radii");
1943 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1944 CTYPE ("Frequency of calculating the Born radii inside rlist");
1945 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1946 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1947 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1948 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1949 CTYPE ("Dielectric coefficient of the implicit solvent");
1950 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1951 CTYPE ("Salt concentration in M for Generalized Born models");
1952 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1953 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1954 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1955 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1956 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1957 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1958 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1959 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1960 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1961 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1963 /* Coupling stuff */
1964 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1965 CTYPE ("Temperature coupling");
1966 EETYPE("tcoupl", ir->etc, etcoupl_names);
1967 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1968 ITYPE("nh-chain-length", ir->opts.nhchainlength, 10);
1969 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1970 CTYPE ("Groups to couple separately");
1971 STYPE ("tc-grps", is->tcgrps, NULL);
1972 CTYPE ("Time constant (ps) and reference temperature (K)");
1973 STYPE ("tau-t", is->tau_t, NULL);
1974 STYPE ("ref-t", is->ref_t, NULL);
1975 CTYPE ("pressure coupling");
1976 EETYPE("pcoupl", ir->epc, epcoupl_names);
1977 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1978 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1979 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1980 RTYPE ("tau-p", ir->tau_p, 1.0);
1981 STYPE ("compressibility", dumstr[0], NULL);
1982 STYPE ("ref-p", dumstr[1], NULL);
1983 CTYPE ("Scaling of reference coordinates, No, All or COM");
1984 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1987 CCTYPE ("OPTIONS FOR QMMM calculations");
1988 EETYPE("QMMM", ir->bQMMM, yesno_names);
1989 CTYPE ("Groups treated Quantum Mechanically");
1990 STYPE ("QMMM-grps", is->QMMM, NULL);
1991 CTYPE ("QM method");
1992 STYPE("QMmethod", is->QMmethod, NULL);
1993 CTYPE ("QMMM scheme");
1994 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1995 CTYPE ("QM basisset");
1996 STYPE("QMbasis", is->QMbasis, NULL);
1997 CTYPE ("QM charge");
1998 STYPE ("QMcharge", is->QMcharge, NULL);
1999 CTYPE ("QM multiplicity");
2000 STYPE ("QMmult", is->QMmult, NULL);
2001 CTYPE ("Surface Hopping");
2002 STYPE ("SH", is->bSH, NULL);
2003 CTYPE ("CAS space options");
2004 STYPE ("CASorbitals", is->CASorbitals, NULL);
2005 STYPE ("CASelectrons", is->CASelectrons, NULL);
2006 STYPE ("SAon", is->SAon, NULL);
2007 STYPE ("SAoff", is->SAoff, NULL);
2008 STYPE ("SAsteps", is->SAsteps, NULL);
2009 CTYPE ("Scale factor for MM charges");
2010 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
2011 CTYPE ("Optimization of QM subsystem");
2012 STYPE ("bOPT", is->bOPT, NULL);
2013 STYPE ("bTS", is->bTS, NULL);
2015 /* Simulated annealing */
2016 CCTYPE("SIMULATED ANNEALING");
2017 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
2018 STYPE ("annealing", is->anneal, NULL);
2019 CTYPE ("Number of time points to use for specifying annealing in each group");
2020 STYPE ("annealing-npoints", is->anneal_npoints, NULL);
2021 CTYPE ("List of times at the annealing points for each group");
2022 STYPE ("annealing-time", is->anneal_time, NULL);
2023 CTYPE ("Temp. at each annealing point, for each group.");
2024 STYPE ("annealing-temp", is->anneal_temp, NULL);
2027 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
2028 EETYPE("gen-vel", opts->bGenVel, yesno_names);
2029 RTYPE ("gen-temp", opts->tempi, 300.0);
2030 ITYPE ("gen-seed", opts->seed, -1);
2033 CCTYPE ("OPTIONS FOR BONDS");
2034 EETYPE("constraints", opts->nshake, constraints);
2035 CTYPE ("Type of constraint algorithm");
2036 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
2037 CTYPE ("Do not constrain the start configuration");
2038 EETYPE("continuation", ir->bContinuation, yesno_names);
2039 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
2040 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
2041 CTYPE ("Relative tolerance of shake");
2042 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
2043 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
2044 ITYPE ("lincs-order", ir->nProjOrder, 4);
2045 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
2046 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
2047 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
2048 ITYPE ("lincs-iter", ir->nLincsIter, 1);
2049 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
2050 CTYPE ("rotates over more degrees than");
2051 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
2052 CTYPE ("Convert harmonic bonds to morse potentials");
2053 EETYPE("morse", opts->bMorse, yesno_names);
2055 /* Energy group exclusions */
2056 CCTYPE ("ENERGY GROUP EXCLUSIONS");
2057 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
2058 STYPE ("energygrp-excl", is->egpexcl, NULL);
2062 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
2063 ITYPE ("nwall", ir->nwall, 0);
2064 EETYPE("wall-type", ir->wall_type, ewt_names);
2065 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
2066 STYPE ("wall-atomtype", is->wall_atomtype, NULL);
2067 STYPE ("wall-density", is->wall_density, NULL);
2068 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
2071 CCTYPE("COM PULLING");
2072 CTYPE("Pull type: no, umbrella, constraint or constant-force");
2073 EETYPE("pull", ir->ePull, epull_names);
2074 if (ir->ePull != epullNO)
2077 is->pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
2080 /* Enforced rotation */
2081 CCTYPE("ENFORCED ROTATION");
2082 CTYPE("Enforced rotation: No or Yes");
2083 EETYPE("rotation", ir->bRot, yesno_names);
2087 is->rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
2090 /* Interactive MD */
2092 CCTYPE("Group to display and/or manipulate in interactive MD session");
2093 STYPE ("IMD-group", is->imd_grp, NULL);
2094 if (is->imd_grp[0] != '\0')
2101 CCTYPE("NMR refinement stuff");
2102 CTYPE ("Distance restraints type: No, Simple or Ensemble");
2103 EETYPE("disre", ir->eDisre, edisre_names);
2104 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
2105 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
2106 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
2107 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
2108 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2109 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2110 CTYPE ("Output frequency for pair distances to energy file");
2111 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2112 CTYPE ("Orientation restraints: No or Yes");
2113 EETYPE("orire", opts->bOrire, yesno_names);
2114 CTYPE ("Orientation restraints force constant and tau for time averaging");
2115 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2116 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2117 STYPE ("orire-fitgrp", is->orirefitgrp, NULL);
2118 CTYPE ("Output frequency for trace(SD) and S to energy file");
2119 ITYPE ("nstorireout", ir->nstorireout, 100);
2121 /* free energy variables */
2122 CCTYPE ("Free energy variables");
2123 EETYPE("free-energy", ir->efep, efep_names);
2124 STYPE ("couple-moltype", is->couple_moltype, NULL);
2125 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2126 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2127 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2129 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2131 it was not entered */
2132 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2133 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2134 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2135 STYPE ("fep-lambdas", is->fep_lambda[efptFEP], NULL);
2136 STYPE ("mass-lambdas", is->fep_lambda[efptMASS], NULL);
2137 STYPE ("coul-lambdas", is->fep_lambda[efptCOUL], NULL);
2138 STYPE ("vdw-lambdas", is->fep_lambda[efptVDW], NULL);
2139 STYPE ("bonded-lambdas", is->fep_lambda[efptBONDED], NULL);
2140 STYPE ("restraint-lambdas", is->fep_lambda[efptRESTRAINT], NULL);
2141 STYPE ("temperature-lambdas", is->fep_lambda[efptTEMPERATURE], NULL);
2142 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2143 STYPE ("init-lambda-weights", is->lambda_weights, NULL);
2144 EETYPE("dhdl-print-energy", fep->edHdLPrintEnergy, edHdLPrintEnergy_names);
2145 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2146 ITYPE ("sc-power", fep->sc_power, 1);
2147 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2148 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2149 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2150 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2151 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2152 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2153 separate_dhdl_file_names);
2154 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2155 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2156 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2158 /* Non-equilibrium MD stuff */
2159 CCTYPE("Non-equilibrium MD stuff");
2160 STYPE ("acc-grps", is->accgrps, NULL);
2161 STYPE ("accelerate", is->acc, NULL);
2162 STYPE ("freezegrps", is->freeze, NULL);
2163 STYPE ("freezedim", is->frdim, NULL);
2164 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2165 STYPE ("deform", is->deform, NULL);
2167 /* simulated tempering variables */
2168 CCTYPE("simulated tempering variables");
2169 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2170 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2171 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2172 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2174 /* expanded ensemble variables */
2175 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2177 read_expandedparams(&ninp, &inp, expand, wi);
2180 /* Electric fields */
2181 CCTYPE("Electric fields");
2182 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2183 CTYPE ("and a phase angle (real)");
2184 STYPE ("E-x", is->efield_x, NULL);
2185 STYPE ("E-xt", is->efield_xt, NULL);
2186 STYPE ("E-y", is->efield_y, NULL);
2187 STYPE ("E-yt", is->efield_yt, NULL);
2188 STYPE ("E-z", is->efield_z, NULL);
2189 STYPE ("E-zt", is->efield_zt, NULL);
2191 CCTYPE("Ion/water position swapping for computational electrophysiology setups");
2192 CTYPE("Swap positions along direction: no, X, Y, Z");
2193 EETYPE("swapcoords", ir->eSwapCoords, eSwapTypes_names);
2194 if (ir->eSwapCoords != eswapNO)
2197 CTYPE("Swap attempt frequency");
2198 ITYPE("swap-frequency", ir->swap->nstswap, 1);
2199 CTYPE("Two index groups that contain the compartment-partitioning atoms");
2200 STYPE("split-group0", splitgrp0, NULL);
2201 STYPE("split-group1", splitgrp1, NULL);
2202 CTYPE("Use center of mass of split groups (yes/no), otherwise center of geometry is used");
2203 EETYPE("massw-split0", ir->swap->massw_split[0], yesno_names);
2204 EETYPE("massw-split1", ir->swap->massw_split[1], yesno_names);
2206 CTYPE("Group name of ions that can be exchanged with solvent molecules");
2207 STYPE("swap-group", swapgrp, NULL);
2208 CTYPE("Group name of solvent molecules");
2209 STYPE("solvent-group", solgrp, NULL);
2211 CTYPE("Split cylinder: radius, upper and lower extension (nm) (this will define the channels)");
2212 CTYPE("Note that the split cylinder settings do not have an influence on the swapping protocol,");
2213 CTYPE("however, if correctly defined, the ion permeation events are counted per channel");
2214 RTYPE("cyl0-r", ir->swap->cyl0r, 2.0);
2215 RTYPE("cyl0-up", ir->swap->cyl0u, 1.0);
2216 RTYPE("cyl0-down", ir->swap->cyl0l, 1.0);
2217 RTYPE("cyl1-r", ir->swap->cyl1r, 2.0);
2218 RTYPE("cyl1-up", ir->swap->cyl1u, 1.0);
2219 RTYPE("cyl1-down", ir->swap->cyl1l, 1.0);
2221 CTYPE("Average the number of ions per compartment over these many swap attempt steps");
2222 ITYPE("coupl-steps", ir->swap->nAverage, 10);
2223 CTYPE("Requested number of anions and cations for each of the two compartments");
2224 CTYPE("-1 means fix the numbers as found in time step 0");
2225 ITYPE("anionsA", ir->swap->nanions[0], -1);
2226 ITYPE("cationsA", ir->swap->ncations[0], -1);
2227 ITYPE("anionsB", ir->swap->nanions[1], -1);
2228 ITYPE("cationsB", ir->swap->ncations[1], -1);
2229 CTYPE("Start to swap ions if threshold difference to requested count is reached");
2230 RTYPE("threshold", ir->swap->threshold, 1.0);
2233 /* AdResS defined thingies */
2234 CCTYPE ("AdResS parameters");
2235 EETYPE("adress", ir->bAdress, yesno_names);
2238 snew(ir->adress, 1);
2239 read_adressparams(&ninp, &inp, ir->adress, wi);
2242 /* User defined thingies */
2243 CCTYPE ("User defined thingies");
2244 STYPE ("user1-grps", is->user1, NULL);
2245 STYPE ("user2-grps", is->user2, NULL);
2246 ITYPE ("userint1", ir->userint1, 0);
2247 ITYPE ("userint2", ir->userint2, 0);
2248 ITYPE ("userint3", ir->userint3, 0);
2249 ITYPE ("userint4", ir->userint4, 0);
2250 RTYPE ("userreal1", ir->userreal1, 0);
2251 RTYPE ("userreal2", ir->userreal2, 0);
2252 RTYPE ("userreal3", ir->userreal3, 0);
2253 RTYPE ("userreal4", ir->userreal4, 0);
2256 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2257 for (i = 0; (i < ninp); i++)
2260 sfree(inp[i].value);
2264 /* Process options if necessary */
2265 for (m = 0; m < 2; m++)
2267 for (i = 0; i < 2*DIM; i++)
2276 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2278 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2280 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2282 case epctSEMIISOTROPIC:
2283 case epctSURFACETENSION:
2284 if (sscanf(dumstr[m], "%lf%lf",
2285 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2287 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2289 dumdub[m][YY] = dumdub[m][XX];
2291 case epctANISOTROPIC:
2292 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2293 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2294 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2296 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2300 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2301 epcoupltype_names[ir->epct]);
2305 clear_mat(ir->ref_p);
2306 clear_mat(ir->compress);
2307 for (i = 0; i < DIM; i++)
2309 ir->ref_p[i][i] = dumdub[1][i];
2310 ir->compress[i][i] = dumdub[0][i];
2312 if (ir->epct == epctANISOTROPIC)
2314 ir->ref_p[XX][YY] = dumdub[1][3];
2315 ir->ref_p[XX][ZZ] = dumdub[1][4];
2316 ir->ref_p[YY][ZZ] = dumdub[1][5];
2317 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2319 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2321 ir->compress[XX][YY] = dumdub[0][3];
2322 ir->compress[XX][ZZ] = dumdub[0][4];
2323 ir->compress[YY][ZZ] = dumdub[0][5];
2324 for (i = 0; i < DIM; i++)
2326 for (m = 0; m < i; m++)
2328 ir->ref_p[i][m] = ir->ref_p[m][i];
2329 ir->compress[i][m] = ir->compress[m][i];
2334 if (ir->comm_mode == ecmNO)
2339 opts->couple_moltype = NULL;
2340 if (strlen(is->couple_moltype) > 0)
2342 if (ir->efep != efepNO)
2344 opts->couple_moltype = strdup(is->couple_moltype);
2345 if (opts->couple_lam0 == opts->couple_lam1)
2347 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2349 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2350 opts->couple_lam1 == ecouplamNONE))
2352 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2357 warning_note(wi, "Free energy is turned off, so we will not decouple the molecule listed in your input.");
2360 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2361 if (ir->efep != efepNO)
2363 if (fep->delta_lambda > 0)
2365 ir->efep = efepSLOWGROWTH;
2369 if (fep->edHdLPrintEnergy == edHdLPrintEnergyYES)
2371 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2372 warning_note(wi, "Old option for dhdl-print-energy given: "
2373 "changing \"yes\" to \"total\"\n");
2376 if (ir->bSimTemp && (fep->edHdLPrintEnergy == edHdLPrintEnergyNO))
2378 /* always print out the energy to dhdl if we are doing
2379 expanded ensemble, since we need the total energy for
2380 analysis if the temperature is changing. In some
2381 conditions one may only want the potential energy, so
2382 we will allow that if the appropriate mdp setting has
2383 been enabled. Otherwise, total it is:
2385 fep->edHdLPrintEnergy = edHdLPrintEnergyTOTAL;
2388 if ((ir->efep != efepNO) || ir->bSimTemp)
2390 ir->bExpanded = FALSE;
2391 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2393 ir->bExpanded = TRUE;
2395 do_fep_params(ir, is->fep_lambda, is->lambda_weights);
2396 if (ir->bSimTemp) /* done after fep params */
2398 do_simtemp_params(ir);
2403 ir->fepvals->n_lambda = 0;
2406 /* WALL PARAMETERS */
2408 do_wall_params(ir, is->wall_atomtype, is->wall_density, opts);
2410 /* ORIENTATION RESTRAINT PARAMETERS */
2412 if (opts->bOrire && str_nelem(is->orirefitgrp, MAXPTR, NULL) != 1)
2414 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2417 /* DEFORMATION PARAMETERS */
2419 clear_mat(ir->deform);
2420 for (i = 0; i < 6; i++)
2424 m = sscanf(is->deform, "%lf %lf %lf %lf %lf %lf",
2425 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2426 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2427 for (i = 0; i < 3; i++)
2429 ir->deform[i][i] = dumdub[0][i];
2431 ir->deform[YY][XX] = dumdub[0][3];
2432 ir->deform[ZZ][XX] = dumdub[0][4];
2433 ir->deform[ZZ][YY] = dumdub[0][5];
2434 if (ir->epc != epcNO)
2436 for (i = 0; i < 3; i++)
2438 for (j = 0; j <= i; j++)
2440 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2442 warning_error(wi, "A box element has deform set and compressibility > 0");
2446 for (i = 0; i < 3; i++)
2448 for (j = 0; j < i; j++)
2450 if (ir->deform[i][j] != 0)
2452 for (m = j; m < DIM; m++)
2454 if (ir->compress[m][j] != 0)
2456 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.");
2457 warning(wi, warn_buf);
2465 /* Ion/water position swapping checks */
2466 if (ir->eSwapCoords != eswapNO)
2468 if (ir->swap->nstswap < 1)
2470 warning_error(wi, "swap_frequency must be 1 or larger when ion swapping is requested");
2472 if (ir->swap->nAverage < 1)
2474 warning_error(wi, "coupl_steps must be 1 or larger.\n");
2476 if (ir->swap->threshold < 1.0)
2478 warning_error(wi, "Ion count threshold must be at least 1.\n");
2486 static int search_QMstring(const char *s, int ng, const char *gn[])
2488 /* same as normal search_string, but this one searches QM strings */
2491 for (i = 0; (i < ng); i++)
2493 if (gmx_strcasecmp(s, gn[i]) == 0)
2499 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2503 } /* search_QMstring */
2505 /* We would like gn to be const as well, but C doesn't allow this */
2506 int search_string(const char *s, int ng, char *gn[])
2510 for (i = 0; (i < ng); i++)
2512 if (gmx_strcasecmp(s, gn[i]) == 0)
2519 "Group %s referenced in the .mdp file was not found in the index file.\n"
2520 "Group names must match either [moleculetype] names or custom index group\n"
2521 "names, in which case you must supply an index file to the '-n' option\n"
2528 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2529 t_blocka *block, char *gnames[],
2530 int gtype, int restnm,
2531 int grptp, gmx_bool bVerbose,
2534 unsigned short *cbuf;
2535 t_grps *grps = &(groups->grps[gtype]);
2536 int i, j, gid, aj, ognr, ntot = 0;
2539 char warn_buf[STRLEN];
2543 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2546 title = gtypes[gtype];
2549 /* Mark all id's as not set */
2550 for (i = 0; (i < natoms); i++)
2555 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2556 for (i = 0; (i < ng); i++)
2558 /* Lookup the group name in the block structure */
2559 gid = search_string(ptrs[i], block->nr, gnames);
2560 if ((grptp != egrptpONE) || (i == 0))
2562 grps->nm_ind[grps->nr++] = gid;
2566 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2569 /* Now go over the atoms in the group */
2570 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2575 /* Range checking */
2576 if ((aj < 0) || (aj >= natoms))
2578 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2580 /* Lookup up the old group number */
2584 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2585 aj+1, title, ognr+1, i+1);
2589 /* Store the group number in buffer */
2590 if (grptp == egrptpONE)
2603 /* Now check whether we have done all atoms */
2607 if (grptp == egrptpALL)
2609 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2610 natoms-ntot, title);
2612 else if (grptp == egrptpPART)
2614 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2615 natoms-ntot, title);
2616 warning_note(wi, warn_buf);
2618 /* Assign all atoms currently unassigned to a rest group */
2619 for (j = 0; (j < natoms); j++)
2621 if (cbuf[j] == NOGID)
2627 if (grptp != egrptpPART)
2632 "Making dummy/rest group for %s containing %d elements\n",
2633 title, natoms-ntot);
2635 /* Add group name "rest" */
2636 grps->nm_ind[grps->nr] = restnm;
2638 /* Assign the rest name to all atoms not currently assigned to a group */
2639 for (j = 0; (j < natoms); j++)
2641 if (cbuf[j] == NOGID)
2650 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2652 /* All atoms are part of one (or no) group, no index required */
2653 groups->ngrpnr[gtype] = 0;
2654 groups->grpnr[gtype] = NULL;
2658 groups->ngrpnr[gtype] = natoms;
2659 snew(groups->grpnr[gtype], natoms);
2660 for (j = 0; (j < natoms); j++)
2662 groups->grpnr[gtype][j] = cbuf[j];
2668 return (bRest && grptp == egrptpPART);
2671 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2674 gmx_groups_t *groups;
2676 int natoms, ai, aj, i, j, d, g, imin, jmin;
2678 int *nrdf2, *na_vcm, na_tot;
2679 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2680 gmx_mtop_atomloop_all_t aloop;
2682 int mb, mol, ftype, as;
2683 gmx_molblock_t *molb;
2684 gmx_moltype_t *molt;
2687 * First calc 3xnr-atoms for each group
2688 * then subtract half a degree of freedom for each constraint
2690 * Only atoms and nuclei contribute to the degrees of freedom...
2695 groups = &mtop->groups;
2696 natoms = mtop->natoms;
2698 /* Allocate one more for a possible rest group */
2699 /* We need to sum degrees of freedom into doubles,
2700 * since floats give too low nrdf's above 3 million atoms.
2702 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2703 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2704 snew(na_vcm, groups->grps[egcVCM].nr+1);
2706 for (i = 0; i < groups->grps[egcTC].nr; i++)
2710 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2715 snew(nrdf2, natoms);
2716 aloop = gmx_mtop_atomloop_all_init(mtop);
2717 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2720 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2722 g = ggrpnr(groups, egcFREEZE, i);
2723 /* Double count nrdf for particle i */
2724 for (d = 0; d < DIM; d++)
2726 if (opts->nFreeze[g][d] == 0)
2731 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2732 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2737 for (mb = 0; mb < mtop->nmolblock; mb++)
2739 molb = &mtop->molblock[mb];
2740 molt = &mtop->moltype[molb->type];
2741 atom = molt->atoms.atom;
2742 for (mol = 0; mol < molb->nmol; mol++)
2744 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2746 ia = molt->ilist[ftype].iatoms;
2747 for (i = 0; i < molt->ilist[ftype].nr; )
2749 /* Subtract degrees of freedom for the constraints,
2750 * if the particles still have degrees of freedom left.
2751 * If one of the particles is a vsite or a shell, then all
2752 * constraint motion will go there, but since they do not
2753 * contribute to the constraints the degrees of freedom do not
2758 if (((atom[ia[1]].ptype == eptNucleus) ||
2759 (atom[ia[1]].ptype == eptAtom)) &&
2760 ((atom[ia[2]].ptype == eptNucleus) ||
2761 (atom[ia[2]].ptype == eptAtom)))
2779 imin = min(imin, nrdf2[ai]);
2780 jmin = min(jmin, nrdf2[aj]);
2783 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2784 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2785 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2786 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2788 ia += interaction_function[ftype].nratoms+1;
2789 i += interaction_function[ftype].nratoms+1;
2792 ia = molt->ilist[F_SETTLE].iatoms;
2793 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2795 /* Subtract 1 dof from every atom in the SETTLE */
2796 for (j = 0; j < 3; j++)
2799 imin = min(2, nrdf2[ai]);
2801 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2802 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2807 as += molt->atoms.nr;
2811 if (ir->ePull == epullCONSTRAINT)
2813 /* Correct nrdf for the COM constraints.
2814 * We correct using the TC and VCM group of the first atom
2815 * in the reference and pull group. If atoms in one pull group
2816 * belong to different TC or VCM groups it is anyhow difficult
2817 * to determine the optimal nrdf assignment.
2821 for (i = 0; i < pull->ncoord; i++)
2825 for (j = 0; j < 2; j++)
2827 const t_pull_group *pgrp;
2829 pgrp = &pull->group[pull->coord[i].group[j]];
2833 /* Subtract 1/2 dof from each group */
2835 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2836 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2837 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2839 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)]]);
2844 /* We need to subtract the whole DOF from group j=1 */
2851 if (ir->nstcomm != 0)
2853 /* Subtract 3 from the number of degrees of freedom in each vcm group
2854 * when com translation is removed and 6 when rotation is removed
2857 switch (ir->comm_mode)
2860 n_sub = ndof_com(ir);
2867 gmx_incons("Checking comm_mode");
2870 for (i = 0; i < groups->grps[egcTC].nr; i++)
2872 /* Count the number of atoms of TC group i for every VCM group */
2873 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2878 for (ai = 0; ai < natoms; ai++)
2880 if (ggrpnr(groups, egcTC, ai) == i)
2882 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2886 /* Correct for VCM removal according to the fraction of each VCM
2887 * group present in this TC group.
2889 nrdf_uc = nrdf_tc[i];
2892 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2896 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2898 if (nrdf_vcm[j] > n_sub)
2900 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2901 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2905 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2906 j, nrdf_vcm[j], nrdf_tc[i]);
2911 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2913 opts->nrdf[i] = nrdf_tc[i];
2914 if (opts->nrdf[i] < 0)
2919 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2920 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2929 static void decode_cos(char *s, t_cosines *cosine)
2932 char format[STRLEN], f1[STRLEN];
2944 sscanf(t, "%d", &(cosine->n));
2951 snew(cosine->a, cosine->n);
2952 snew(cosine->phi, cosine->n);
2954 sprintf(format, "%%*d");
2955 for (i = 0; (i < cosine->n); i++)
2958 strcat(f1, "%lf%lf");
2959 if (sscanf(t, f1, &a, &phi) < 2)
2961 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2964 cosine->phi[i] = phi;
2965 strcat(format, "%*lf%*lf");
2972 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2973 const char *option, const char *val, int flag)
2975 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2976 * But since this is much larger than STRLEN, such a line can not be parsed.
2977 * The real maximum is the number of names that fit in a string: STRLEN/2.
2979 #define EGP_MAX (STRLEN/2)
2980 int nelem, i, j, k, nr;
2981 char *names[EGP_MAX];
2985 gnames = groups->grpname;
2987 nelem = str_nelem(val, EGP_MAX, names);
2990 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2992 nr = groups->grps[egcENER].nr;
2994 for (i = 0; i < nelem/2; i++)
2998 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
3004 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3005 names[2*i], option);
3009 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
3015 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
3016 names[2*i+1], option);
3018 if ((j < nr) && (k < nr))
3020 ir->opts.egp_flags[nr*j+k] |= flag;
3021 ir->opts.egp_flags[nr*k+j] |= flag;
3030 static void make_swap_groups(
3039 int ig = -1, i = 0, j;
3043 /* Just a quick check here, more thorough checks are in mdrun */
3044 if (strcmp(splitg0name, splitg1name) == 0)
3046 gmx_fatal(FARGS, "The split groups can not both be '%s'.", splitg0name);
3049 /* First get the swap group index atoms */
3050 ig = search_string(swapgname, grps->nr, gnames);
3051 swap->nat = grps->index[ig+1] - grps->index[ig];
3054 fprintf(stderr, "Swap group '%s' contains %d atoms.\n", swapgname, swap->nat);
3055 snew(swap->ind, swap->nat);
3056 for (i = 0; i < swap->nat; i++)
3058 swap->ind[i] = grps->a[grps->index[ig]+i];
3063 gmx_fatal(FARGS, "You defined an empty group of atoms for swapping.");
3066 /* Now do so for the split groups */
3067 for (j = 0; j < 2; j++)
3071 splitg = splitg0name;
3075 splitg = splitg1name;
3078 ig = search_string(splitg, grps->nr, gnames);
3079 swap->nat_split[j] = grps->index[ig+1] - grps->index[ig];
3080 if (swap->nat_split[j] > 0)
3082 fprintf(stderr, "Split group %d '%s' contains %d atom%s.\n",
3083 j, splitg, swap->nat_split[j], (swap->nat_split[j] > 1) ? "s" : "");
3084 snew(swap->ind_split[j], swap->nat_split[j]);
3085 for (i = 0; i < swap->nat_split[j]; i++)
3087 swap->ind_split[j][i] = grps->a[grps->index[ig]+i];
3092 gmx_fatal(FARGS, "Split group %d has to contain at least 1 atom!", j);
3096 /* Now get the solvent group index atoms */
3097 ig = search_string(solgname, grps->nr, gnames);
3098 swap->nat_sol = grps->index[ig+1] - grps->index[ig];
3099 if (swap->nat_sol > 0)
3101 fprintf(stderr, "Solvent group '%s' contains %d atoms.\n", solgname, swap->nat_sol);
3102 snew(swap->ind_sol, swap->nat_sol);
3103 for (i = 0; i < swap->nat_sol; i++)
3105 swap->ind_sol[i] = grps->a[grps->index[ig]+i];
3110 gmx_fatal(FARGS, "You defined an empty group of solvent. Cannot exchange ions.");
3115 void make_IMD_group(t_IMD *IMDgroup, char *IMDgname, t_blocka *grps, char **gnames)
3120 ig = search_string(IMDgname, grps->nr, gnames);
3121 IMDgroup->nat = grps->index[ig+1] - grps->index[ig];
3123 if (IMDgroup->nat > 0)
3125 fprintf(stderr, "Group '%s' with %d atoms can be activated for interactive molecular dynamics (IMD).\n",
3126 IMDgname, IMDgroup->nat);
3127 snew(IMDgroup->ind, IMDgroup->nat);
3128 for (i = 0; i < IMDgroup->nat; i++)
3130 IMDgroup->ind[i] = grps->a[grps->index[ig]+i];
3136 void do_index(const char* mdparin, const char *ndx,
3139 t_inputrec *ir, rvec *v,
3143 gmx_groups_t *groups;
3147 char warnbuf[STRLEN], **gnames;
3148 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
3151 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
3152 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
3153 int i, j, k, restnm;
3155 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
3156 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
3157 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
3158 char warn_buf[STRLEN];
3162 fprintf(stderr, "processing index file...\n");
3168 snew(grps->index, 1);
3170 atoms_all = gmx_mtop_global_atoms(mtop);
3171 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
3172 free_t_atoms(&atoms_all, FALSE);
3176 grps = init_index(ndx, &gnames);
3179 groups = &mtop->groups;
3180 natoms = mtop->natoms;
3181 symtab = &mtop->symtab;
3183 snew(groups->grpname, grps->nr+1);
3185 for (i = 0; (i < grps->nr); i++)
3187 groups->grpname[i] = put_symtab(symtab, gnames[i]);
3189 groups->grpname[i] = put_symtab(symtab, "rest");
3191 srenew(gnames, grps->nr+1);
3192 gnames[restnm] = *(groups->grpname[i]);
3193 groups->ngrpname = grps->nr+1;
3195 set_warning_line(wi, mdparin, -1);
3197 ntau_t = str_nelem(is->tau_t, MAXPTR, ptr1);
3198 nref_t = str_nelem(is->ref_t, MAXPTR, ptr2);
3199 ntcg = str_nelem(is->tcgrps, MAXPTR, ptr3);
3200 if ((ntau_t != ntcg) || (nref_t != ntcg))
3202 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
3203 "%d tau-t values", ntcg, nref_t, ntau_t);
3206 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
3207 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
3208 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
3209 nr = groups->grps[egcTC].nr;
3211 snew(ir->opts.nrdf, nr);
3212 snew(ir->opts.tau_t, nr);
3213 snew(ir->opts.ref_t, nr);
3214 if (ir->eI == eiBD && ir->bd_fric == 0)
3216 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
3223 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
3227 for (i = 0; (i < nr); i++)
3229 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
3230 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
3232 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
3233 warning_error(wi, warn_buf);
3236 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
3238 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.");
3241 if (ir->opts.tau_t[i] >= 0)
3243 tau_min = min(tau_min, ir->opts.tau_t[i]);
3246 if (ir->etc != etcNO && ir->nsttcouple == -1)
3248 ir->nsttcouple = ir_optimal_nsttcouple(ir);
3253 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
3255 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");
3257 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
3259 if (ir->nstpcouple != ir->nsttcouple)
3261 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
3262 ir->nstpcouple = ir->nsttcouple = mincouple;
3263 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);
3264 warning_note(wi, warn_buf);
3268 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
3269 primarily for testing purposes, and does not work with temperature coupling other than 1 */
3271 if (ETC_ANDERSEN(ir->etc))
3273 if (ir->nsttcouple != 1)
3276 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");
3277 warning_note(wi, warn_buf);
3280 nstcmin = tcouple_min_integration_steps(ir->etc);
3283 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin - 10*GMX_REAL_EPS)
3285 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3286 ETCOUPLTYPE(ir->etc),
3288 ir->nsttcouple*ir->delta_t);
3289 warning(wi, warn_buf);
3292 for (i = 0; (i < nr); i++)
3294 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3295 if (ir->opts.ref_t[i] < 0)
3297 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3300 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3301 if we are in this conditional) if mc_temp is negative */
3302 if (ir->expandedvals->mc_temp < 0)
3304 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3308 /* Simulated annealing for each group. There are nr groups */
3309 nSA = str_nelem(is->anneal, MAXPTR, ptr1);
3310 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3314 if (nSA > 0 && nSA != nr)
3316 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3320 snew(ir->opts.annealing, nr);
3321 snew(ir->opts.anneal_npoints, nr);
3322 snew(ir->opts.anneal_time, nr);
3323 snew(ir->opts.anneal_temp, nr);
3324 for (i = 0; i < nr; i++)
3326 ir->opts.annealing[i] = eannNO;
3327 ir->opts.anneal_npoints[i] = 0;
3328 ir->opts.anneal_time[i] = NULL;
3329 ir->opts.anneal_temp[i] = NULL;
3334 for (i = 0; i < nr; i++)
3336 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3338 ir->opts.annealing[i] = eannNO;
3340 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3342 ir->opts.annealing[i] = eannSINGLE;
3345 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3347 ir->opts.annealing[i] = eannPERIODIC;
3353 /* Read the other fields too */
3354 nSA_points = str_nelem(is->anneal_npoints, MAXPTR, ptr1);
3355 if (nSA_points != nSA)
3357 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3359 for (k = 0, i = 0; i < nr; i++)
3361 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3362 if (ir->opts.anneal_npoints[i] == 1)
3364 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3366 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3367 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3368 k += ir->opts.anneal_npoints[i];
3371 nSA_time = str_nelem(is->anneal_time, MAXPTR, ptr1);
3374 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3376 nSA_temp = str_nelem(is->anneal_temp, MAXPTR, ptr2);
3379 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3382 for (i = 0, k = 0; i < nr; i++)
3385 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3387 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3388 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3391 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3393 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3399 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3401 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3402 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3405 if (ir->opts.anneal_temp[i][j] < 0)
3407 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3412 /* Print out some summary information, to make sure we got it right */
3413 for (i = 0, k = 0; i < nr; i++)
3415 if (ir->opts.annealing[i] != eannNO)
3417 j = groups->grps[egcTC].nm_ind[i];
3418 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3419 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3420 ir->opts.anneal_npoints[i]);
3421 fprintf(stderr, "Time (ps) Temperature (K)\n");
3422 /* All terms except the last one */
3423 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3425 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3428 /* Finally the last one */
3429 j = ir->opts.anneal_npoints[i]-1;
3430 if (ir->opts.annealing[i] == eannSINGLE)
3432 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3436 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3437 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3439 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3448 if (ir->ePull != epullNO)
3450 make_pull_groups(ir->pull, is->pull_grp, grps, gnames);
3452 make_pull_coords(ir->pull);
3457 make_rotation_groups(ir->rot, is->rot_grp, grps, gnames);
3460 if (ir->eSwapCoords != eswapNO)
3462 make_swap_groups(ir->swap, swapgrp, splitgrp0, splitgrp1, solgrp, grps, gnames);
3465 /* Make indices for IMD session */
3468 make_IMD_group(ir->imd, is->imd_grp, grps, gnames);
3471 nacc = str_nelem(is->acc, MAXPTR, ptr1);
3472 nacg = str_nelem(is->accgrps, MAXPTR, ptr2);
3473 if (nacg*DIM != nacc)
3475 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3478 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3479 restnm, egrptpALL_GENREST, bVerbose, wi);
3480 nr = groups->grps[egcACC].nr;
3481 snew(ir->opts.acc, nr);
3482 ir->opts.ngacc = nr;
3484 for (i = k = 0; (i < nacg); i++)
3486 for (j = 0; (j < DIM); j++, k++)
3488 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3491 for (; (i < nr); i++)
3493 for (j = 0; (j < DIM); j++)
3495 ir->opts.acc[i][j] = 0;
3499 nfrdim = str_nelem(is->frdim, MAXPTR, ptr1);
3500 nfreeze = str_nelem(is->freeze, MAXPTR, ptr2);
3501 if (nfrdim != DIM*nfreeze)
3503 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3506 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3507 restnm, egrptpALL_GENREST, bVerbose, wi);
3508 nr = groups->grps[egcFREEZE].nr;
3509 ir->opts.ngfrz = nr;
3510 snew(ir->opts.nFreeze, nr);
3511 for (i = k = 0; (i < nfreeze); i++)
3513 for (j = 0; (j < DIM); j++, k++)
3515 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3516 if (!ir->opts.nFreeze[i][j])
3518 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3520 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3521 "(not %s)", ptr1[k]);
3522 warning(wi, warn_buf);
3527 for (; (i < nr); i++)
3529 for (j = 0; (j < DIM); j++)
3531 ir->opts.nFreeze[i][j] = 0;
3535 nenergy = str_nelem(is->energy, MAXPTR, ptr1);
3536 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3537 restnm, egrptpALL_GENREST, bVerbose, wi);
3538 add_wall_energrps(groups, ir->nwall, symtab);
3539 ir->opts.ngener = groups->grps[egcENER].nr;
3540 nvcm = str_nelem(is->vcm, MAXPTR, ptr1);
3542 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3543 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3546 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3547 "This may lead to artifacts.\n"
3548 "In most cases one should use one group for the whole system.");
3551 /* Now we have filled the freeze struct, so we can calculate NRDF */
3552 calc_nrdf(mtop, ir, gnames);
3558 /* Must check per group! */
3559 for (i = 0; (i < ir->opts.ngtc); i++)
3561 ntot += ir->opts.nrdf[i];
3563 if (ntot != (DIM*natoms))
3565 fac = sqrt(ntot/(DIM*natoms));
3568 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3569 "and removal of center of mass motion\n", fac);
3571 for (i = 0; (i < natoms); i++)
3573 svmul(fac, v[i], v[i]);
3578 nuser = str_nelem(is->user1, MAXPTR, ptr1);
3579 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3580 restnm, egrptpALL_GENREST, bVerbose, wi);
3581 nuser = str_nelem(is->user2, MAXPTR, ptr1);
3582 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3583 restnm, egrptpALL_GENREST, bVerbose, wi);
3584 nuser = str_nelem(is->x_compressed_groups, MAXPTR, ptr1);
3585 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcCompressedX,
3586 restnm, egrptpONE, bVerbose, wi);
3587 nofg = str_nelem(is->orirefitgrp, MAXPTR, ptr1);
3588 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3589 restnm, egrptpALL_GENREST, bVerbose, wi);
3591 /* QMMM input processing */
3592 nQMg = str_nelem(is->QMMM, MAXPTR, ptr1);
3593 nQMmethod = str_nelem(is->QMmethod, MAXPTR, ptr2);
3594 nQMbasis = str_nelem(is->QMbasis, MAXPTR, ptr3);
3595 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3597 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3598 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3600 /* group rest, if any, is always MM! */
3601 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3602 restnm, egrptpALL_GENREST, bVerbose, wi);
3603 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3604 ir->opts.ngQM = nQMg;
3605 snew(ir->opts.QMmethod, nr);
3606 snew(ir->opts.QMbasis, nr);
3607 for (i = 0; i < nr; i++)
3609 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3610 * converted to the corresponding enum in names.c
3612 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3614 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3618 nQMmult = str_nelem(is->QMmult, MAXPTR, ptr1);
3619 nQMcharge = str_nelem(is->QMcharge, MAXPTR, ptr2);
3620 nbSH = str_nelem(is->bSH, MAXPTR, ptr3);
3621 snew(ir->opts.QMmult, nr);
3622 snew(ir->opts.QMcharge, nr);
3623 snew(ir->opts.bSH, nr);
3625 for (i = 0; i < nr; i++)
3627 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3628 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3629 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3632 nCASelec = str_nelem(is->CASelectrons, MAXPTR, ptr1);
3633 nCASorb = str_nelem(is->CASorbitals, MAXPTR, ptr2);
3634 snew(ir->opts.CASelectrons, nr);
3635 snew(ir->opts.CASorbitals, nr);
3636 for (i = 0; i < nr; i++)
3638 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3639 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3641 /* special optimization options */
3643 nbOPT = str_nelem(is->bOPT, MAXPTR, ptr1);
3644 nbTS = str_nelem(is->bTS, MAXPTR, ptr2);
3645 snew(ir->opts.bOPT, nr);
3646 snew(ir->opts.bTS, nr);
3647 for (i = 0; i < nr; i++)
3649 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3650 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3652 nSAon = str_nelem(is->SAon, MAXPTR, ptr1);
3653 nSAoff = str_nelem(is->SAoff, MAXPTR, ptr2);
3654 nSAsteps = str_nelem(is->SAsteps, MAXPTR, ptr3);
3655 snew(ir->opts.SAon, nr);
3656 snew(ir->opts.SAoff, nr);
3657 snew(ir->opts.SAsteps, nr);
3659 for (i = 0; i < nr; i++)
3661 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3662 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3663 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3665 /* end of QMMM input */
3669 for (i = 0; (i < egcNR); i++)
3671 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3672 for (j = 0; (j < groups->grps[i].nr); j++)
3674 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3676 fprintf(stderr, "\n");
3680 nr = groups->grps[egcENER].nr;
3681 snew(ir->opts.egp_flags, nr*nr);
3683 bExcl = do_egp_flag(ir, groups, "energygrp-excl", is->egpexcl, EGP_EXCL);
3684 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3686 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3688 if (bExcl && EEL_FULL(ir->coulombtype))
3690 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3693 bTable = do_egp_flag(ir, groups, "energygrp-table", is->egptable, EGP_TABLE);
3694 if (bTable && !(ir->vdwtype == evdwUSER) &&
3695 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3696 !(ir->coulombtype == eelPMEUSERSWITCH))
3698 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3701 decode_cos(is->efield_x, &(ir->ex[XX]));
3702 decode_cos(is->efield_xt, &(ir->et[XX]));
3703 decode_cos(is->efield_y, &(ir->ex[YY]));
3704 decode_cos(is->efield_yt, &(ir->et[YY]));
3705 decode_cos(is->efield_z, &(ir->ex[ZZ]));
3706 decode_cos(is->efield_zt, &(ir->et[ZZ]));
3710 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3713 for (i = 0; (i < grps->nr); i++)
3725 static void check_disre(gmx_mtop_t *mtop)
3727 gmx_ffparams_t *ffparams;
3728 t_functype *functype;
3730 int i, ndouble, ftype;
3731 int label, old_label;
3733 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3735 ffparams = &mtop->ffparams;
3736 functype = ffparams->functype;
3737 ip = ffparams->iparams;
3740 for (i = 0; i < ffparams->ntypes; i++)
3742 ftype = functype[i];
3743 if (ftype == F_DISRES)
3745 label = ip[i].disres.label;
3746 if (label == old_label)
3748 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3756 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3757 "probably the parameters for multiple pairs in one restraint "
3758 "are not identical\n", ndouble);
3763 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3764 gmx_bool posres_only,
3768 gmx_mtop_ilistloop_t iloop;
3778 for (d = 0; d < DIM; d++)
3780 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3782 /* Check for freeze groups */
3783 for (g = 0; g < ir->opts.ngfrz; g++)
3785 for (d = 0; d < DIM; d++)
3787 if (ir->opts.nFreeze[g][d] != 0)
3795 /* Check for position restraints */
3796 iloop = gmx_mtop_ilistloop_init(sys);
3797 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3800 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3802 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3804 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3805 for (d = 0; d < DIM; d++)
3807 if (pr->posres.fcA[d] != 0)
3813 for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
3815 /* Check for flat-bottom posres */
3816 pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
3817 if (pr->fbposres.k != 0)
3819 switch (pr->fbposres.geom)
3821 case efbposresSPHERE:
3822 AbsRef[XX] = AbsRef[YY] = AbsRef[ZZ] = 1;
3824 case efbposresCYLINDER:
3825 AbsRef[XX] = AbsRef[YY] = 1;
3827 case efbposresX: /* d=XX */
3828 case efbposresY: /* d=YY */
3829 case efbposresZ: /* d=ZZ */
3830 d = pr->fbposres.geom - efbposresX;
3834 gmx_fatal(FARGS, " Invalid geometry for flat-bottom position restraint.\n"
3835 "Expected nr between 1 and %d. Found %d\n", efbposresNR-1,
3843 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3847 check_combination_rule_differences(const gmx_mtop_t *mtop, int state,
3848 gmx_bool *bC6ParametersWorkWithGeometricRules,
3849 gmx_bool *bC6ParametersWorkWithLBRules,
3850 gmx_bool *bLBRulesPossible)
3852 int ntypes, tpi, tpj, thisLBdiff, thisgeomdiff;
3855 double geometricdiff, LBdiff;
3856 double c6i, c6j, c12i, c12j;
3857 double c6, c6_geometric, c6_LB;
3858 double sigmai, sigmaj, epsi, epsj;
3859 gmx_bool bCanDoLBRules, bCanDoGeometricRules;
3862 /* A tolerance of 1e-5 seems reasonable for (possibly hand-typed)
3863 * force-field floating point parameters.
3866 ptr = getenv("GMX_LJCOMB_TOL");
3871 sscanf(ptr, "%lf", &dbl);
3875 *bC6ParametersWorkWithLBRules = TRUE;
3876 *bC6ParametersWorkWithGeometricRules = TRUE;
3877 bCanDoLBRules = TRUE;
3878 bCanDoGeometricRules = TRUE;
3879 ntypes = mtop->ffparams.atnr;
3880 snew(typecount, ntypes);
3881 gmx_mtop_count_atomtypes(mtop, state, typecount);
3882 geometricdiff = LBdiff = 0.0;
3883 *bLBRulesPossible = TRUE;
3884 for (tpi = 0; tpi < ntypes; ++tpi)
3886 c6i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c6;
3887 c12i = mtop->ffparams.iparams[(ntypes + 1) * tpi].lj.c12;
3888 for (tpj = tpi; tpj < ntypes; ++tpj)
3890 c6j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c6;
3891 c12j = mtop->ffparams.iparams[(ntypes + 1) * tpj].lj.c12;
3892 c6 = mtop->ffparams.iparams[ntypes * tpi + tpj].lj.c6;
3893 c6_geometric = sqrt(c6i * c6j);
3894 if (!gmx_numzero(c6_geometric))
3896 if (!gmx_numzero(c12i) && !gmx_numzero(c12j))
3898 sigmai = pow(c12i / c6i, 1.0/6.0);
3899 sigmaj = pow(c12j / c6j, 1.0/6.0);
3900 epsi = c6i * c6i /(4.0 * c12i);
3901 epsj = c6j * c6j /(4.0 * c12j);
3902 c6_LB = 4.0 * pow(epsi * epsj, 1.0/2.0) * pow(0.5 * (sigmai + sigmaj), 6);
3906 *bLBRulesPossible = FALSE;
3907 c6_LB = c6_geometric;
3909 bCanDoLBRules = gmx_within_tol(c6_LB, c6, tol);
3912 if (FALSE == bCanDoLBRules)
3914 *bC6ParametersWorkWithLBRules = FALSE;
3917 bCanDoGeometricRules = gmx_within_tol(c6_geometric, c6, tol);
3919 if (FALSE == bCanDoGeometricRules)
3921 *bC6ParametersWorkWithGeometricRules = FALSE;
3929 check_combination_rules(const t_inputrec *ir, const gmx_mtop_t *mtop,
3933 gmx_bool bLBRulesPossible, bC6ParametersWorkWithGeometricRules, bC6ParametersWorkWithLBRules;
3935 check_combination_rule_differences(mtop, 0,
3936 &bC6ParametersWorkWithGeometricRules,
3937 &bC6ParametersWorkWithLBRules,
3939 if (ir->ljpme_combination_rule == eljpmeLB)
3941 if (FALSE == bC6ParametersWorkWithLBRules || FALSE == bLBRulesPossible)
3943 warning(wi, "You are using arithmetic-geometric combination rules "
3944 "in LJ-PME, but your non-bonded C6 parameters do not "
3945 "follow these rules.");
3950 if (FALSE == bC6ParametersWorkWithGeometricRules)
3952 if (ir->eDispCorr != edispcNO)
3954 warning_note(wi, "You are using geometric combination rules in "
3955 "LJ-PME, but your non-bonded C6 parameters do "
3956 "not follow these rules. "
3957 "This will introduce very small errors in the forces and energies in "
3958 "your simulations. Dispersion correction will correct total energy "
3959 "and/or pressure for isotropic systems, but not forces or surface tensions.");
3963 warning_note(wi, "You are using geometric combination rules in "
3964 "LJ-PME, but your non-bonded C6 parameters do "
3965 "not follow these rules. "
3966 "This will introduce very small errors in the forces and energies in "
3967 "your simulations. If your system is homogeneous, consider using dispersion correction "
3968 "for the total energy and pressure.");
3974 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3977 char err_buf[STRLEN];
3978 int i, m, c, nmol, npct;
3979 gmx_bool bCharge, bAcc;
3980 real gdt_max, *mgrp, mt;
3982 gmx_mtop_atomloop_block_t aloopb;
3983 gmx_mtop_atomloop_all_t aloop;
3986 char warn_buf[STRLEN];
3988 set_warning_line(wi, mdparin, -1);
3990 if (ir->cutoff_scheme == ecutsVERLET &&
3991 ir->verletbuf_tol > 0 &&
3993 ((EI_MD(ir->eI) || EI_SD(ir->eI)) &&
3994 (ir->etc == etcVRESCALE || ir->etc == etcBERENDSEN)))
3996 /* Check if a too small Verlet buffer might potentially
3997 * cause more drift than the thermostat can couple off.
3999 /* Temperature error fraction for warning and suggestion */
4000 const real T_error_warn = 0.002;
4001 const real T_error_suggest = 0.001;
4002 /* For safety: 2 DOF per atom (typical with constraints) */
4003 const real nrdf_at = 2;
4004 real T, tau, max_T_error;
4009 for (i = 0; i < ir->opts.ngtc; i++)
4011 T = max(T, ir->opts.ref_t[i]);
4012 tau = max(tau, ir->opts.tau_t[i]);
4016 /* This is a worst case estimate of the temperature error,
4017 * assuming perfect buffer estimation and no cancelation
4018 * of errors. The factor 0.5 is because energy distributes
4019 * equally over Ekin and Epot.
4021 max_T_error = 0.5*tau*ir->verletbuf_tol/(nrdf_at*BOLTZ*T);
4022 if (max_T_error > T_error_warn)
4024 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.",
4025 ir->verletbuf_tol, T, tau,
4027 100*T_error_suggest,
4028 ir->verletbuf_tol*T_error_suggest/max_T_error);
4029 warning(wi, warn_buf);
4034 if (ETC_ANDERSEN(ir->etc))
4038 for (i = 0; i < ir->opts.ngtc; i++)
4040 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
4041 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
4042 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
4043 i, ir->opts.tau_t[i]);
4044 CHECK(ir->opts.tau_t[i] < 0);
4047 for (i = 0; i < ir->opts.ngtc; i++)
4049 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
4050 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);
4051 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
4055 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
4056 ir->comm_mode == ecmNO &&
4057 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
4058 !ETC_ANDERSEN(ir->etc))
4060 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");
4063 /* Check for pressure coupling with absolute position restraints */
4064 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
4066 absolute_reference(ir, sys, TRUE, AbsRef);
4068 for (m = 0; m < DIM; m++)
4070 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
4072 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
4080 aloopb = gmx_mtop_atomloop_block_init(sys);
4081 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
4083 if (atom->q != 0 || atom->qB != 0)
4091 if (EEL_FULL(ir->coulombtype))
4094 "You are using full electrostatics treatment %s for a system without charges.\n"
4095 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
4096 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
4097 warning(wi, err_buf);
4102 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
4105 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
4106 "You might want to consider using %s electrostatics.\n",
4108 warning_note(wi, err_buf);
4112 /* Check if combination rules used in LJ-PME are the same as in the force field */
4113 if (EVDW_PME(ir->vdwtype))
4115 check_combination_rules(ir, sys, wi);
4118 /* Generalized reaction field */
4119 if (ir->opts.ngtc == 0)
4121 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
4123 CHECK(ir->coulombtype == eelGRF);
4127 sprintf(err_buf, "When using coulombtype = %s"
4128 " ref-t for temperature coupling should be > 0",
4130 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
4133 if (ir->eI == eiSD2)
4135 sprintf(warn_buf, "The stochastic dynamics integrator %s is deprecated, since\n"
4136 "it is slower than integrator %s and is slightly less accurate\n"
4137 "with constraints. Use the %s integrator.",
4138 ei_names[ir->eI], ei_names[eiSD1], ei_names[eiSD1]);
4139 warning_note(wi, warn_buf);
4143 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4145 for (m = 0; (m < DIM); m++)
4147 if (fabs(ir->opts.acc[i][m]) > 1e-6)
4156 snew(mgrp, sys->groups.grps[egcACC].nr);
4157 aloop = gmx_mtop_atomloop_all_init(sys);
4158 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
4160 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
4163 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4165 for (m = 0; (m < DIM); m++)
4167 acc[m] += ir->opts.acc[i][m]*mgrp[i];
4171 for (m = 0; (m < DIM); m++)
4173 if (fabs(acc[m]) > 1e-6)
4175 const char *dim[DIM] = { "X", "Y", "Z" };
4177 "Net Acceleration in %s direction, will %s be corrected\n",
4178 dim[m], ir->nstcomm != 0 ? "" : "not");
4179 if (ir->nstcomm != 0 && m < ndof_com(ir))
4182 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
4184 ir->opts.acc[i][m] -= acc[m];
4192 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
4193 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
4195 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
4198 if (ir->ePull != epullNO)
4200 gmx_bool bPullAbsoluteRef;
4202 bPullAbsoluteRef = FALSE;
4203 for (i = 0; i < ir->pull->ncoord; i++)
4205 bPullAbsoluteRef = bPullAbsoluteRef ||
4206 ir->pull->coord[i].group[0] == 0 ||
4207 ir->pull->coord[i].group[1] == 0;
4209 if (bPullAbsoluteRef)
4211 absolute_reference(ir, sys, FALSE, AbsRef);
4212 for (m = 0; m < DIM; m++)
4214 if (ir->pull->dim[m] && !AbsRef[m])
4216 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.");
4222 if (ir->pull->eGeom == epullgDIRPBC)
4224 for (i = 0; i < 3; i++)
4226 for (m = 0; m <= i; m++)
4228 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
4229 ir->deform[i][m] != 0)
4231 for (c = 0; c < ir->pull->ncoord; c++)
4233 if (ir->pull->coord[c].vec[m] != 0)
4235 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
4247 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
4251 char warn_buf[STRLEN];
4254 ptr = check_box(ir->ePBC, box);
4257 warning_error(wi, ptr);
4260 if (bConstr && ir->eConstrAlg == econtSHAKE)
4262 if (ir->shake_tol <= 0.0)
4264 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
4266 warning_error(wi, warn_buf);
4269 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
4271 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
4272 if (ir->epc == epcNO)
4274 warning(wi, warn_buf);
4278 warning_error(wi, warn_buf);
4283 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
4285 /* If we have Lincs constraints: */
4286 if (ir->eI == eiMD && ir->etc == etcNO &&
4287 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
4289 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
4290 warning_note(wi, warn_buf);
4293 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
4295 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
4296 warning_note(wi, warn_buf);
4298 if (ir->epc == epcMTTK)
4300 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
4304 if (bConstr && ir->epc == epcMTTK)
4306 warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
4309 if (ir->LincsWarnAngle > 90.0)
4311 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
4312 warning(wi, warn_buf);
4313 ir->LincsWarnAngle = 90.0;
4316 if (ir->ePBC != epbcNONE)
4318 if (ir->nstlist == 0)
4320 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
4322 bTWIN = (ir->rlistlong > ir->rlist);
4323 if (ir->ns_type == ensGRID)
4325 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
4327 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",
4328 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
4329 warning_error(wi, warn_buf);
4334 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
4335 if (2*ir->rlistlong >= min_size)
4337 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
4338 warning_error(wi, warn_buf);
4341 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
4348 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
4352 real rvdw1, rvdw2, rcoul1, rcoul2;
4353 char warn_buf[STRLEN];
4355 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
4359 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
4364 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
4370 if (rvdw1 + rvdw2 > ir->rlist ||
4371 rcoul1 + rcoul2 > ir->rlist)
4374 "The sum of the two largest charge group radii (%f) "
4375 "is larger than rlist (%f)\n",
4376 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
4377 warning(wi, warn_buf);
4381 /* Here we do not use the zero at cut-off macro,
4382 * since user defined interactions might purposely
4383 * not be zero at the cut-off.
4385 if (ir_vdw_is_zero_at_cutoff(ir) &&
4386 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
4388 sprintf(warn_buf, "The sum of the two largest charge group "
4389 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
4390 "With exact cut-offs, better performance can be "
4391 "obtained with cutoff-scheme = %s, because it "
4392 "does not use charge groups at all.",
4394 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4395 ir->rlistlong, ir->rvdw,
4396 ecutscheme_names[ecutsVERLET]);
4399 warning(wi, warn_buf);
4403 warning_note(wi, warn_buf);
4406 if (ir_coulomb_is_zero_at_cutoff(ir) &&
4407 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
4409 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
4410 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
4412 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
4413 ir->rlistlong, ir->rcoulomb,
4414 ecutscheme_names[ecutsVERLET]);
4417 warning(wi, warn_buf);
4421 warning_note(wi, warn_buf);