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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
50 #include "gmx_fatal.h"
63 #include "mtop_util.h"
64 #include "chargegroup.h"
69 #define MAXLAMBDAS 1024
71 /* Resource parameters
72 * Do not change any of these until you read the instruction
73 * in readinp.h. Some cpp's do not take spaces after the backslash
74 * (like the c-shell), which will give you a very weird compiler
78 static char tcgrps[STRLEN], tau_t[STRLEN], ref_t[STRLEN],
79 acc[STRLEN], accgrps[STRLEN], freeze[STRLEN], frdim[STRLEN],
80 energy[STRLEN], user1[STRLEN], user2[STRLEN], vcm[STRLEN], xtc_grps[STRLEN],
81 couple_moltype[STRLEN], orirefitgrp[STRLEN], egptable[STRLEN], egpexcl[STRLEN],
82 wall_atomtype[STRLEN], wall_density[STRLEN], deform[STRLEN], QMMM[STRLEN];
83 static char fep_lambda[efptNR][STRLEN];
84 static char lambda_weights[STRLEN];
85 static char **pull_grp;
86 static char **rot_grp;
87 static char anneal[STRLEN], anneal_npoints[STRLEN],
88 anneal_time[STRLEN], anneal_temp[STRLEN];
89 static char QMmethod[STRLEN], QMbasis[STRLEN], QMcharge[STRLEN], QMmult[STRLEN],
90 bSH[STRLEN], CASorbitals[STRLEN], CASelectrons[STRLEN], SAon[STRLEN],
91 SAoff[STRLEN], SAsteps[STRLEN], bTS[STRLEN], bOPT[STRLEN];
92 static char efield_x[STRLEN], efield_xt[STRLEN], efield_y[STRLEN],
93 efield_yt[STRLEN], efield_z[STRLEN], efield_zt[STRLEN];
96 egrptpALL, /* All particles have to be a member of a group. */
97 egrptpALL_GENREST, /* A rest group with name is generated for particles *
98 * that are not part of any group. */
99 egrptpPART, /* As egrptpALL_GENREST, but no name is generated *
100 * for the rest group. */
101 egrptpONE /* Merge all selected groups into one group, *
102 * make a rest group for the remaining particles. */
106 void init_ir(t_inputrec *ir, t_gromppopts *opts)
108 snew(opts->include, STRLEN);
109 snew(opts->define, STRLEN);
110 snew(ir->fepvals, 1);
111 snew(ir->expandedvals, 1);
112 snew(ir->simtempvals, 1);
115 static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
120 for (i = 0; i < ntemps; i++)
122 /* simple linear scaling -- allows more control */
123 if (simtemp->eSimTempScale == esimtempLINEAR)
125 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
127 else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
129 simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low, (1.0*i)/(ntemps-1));
131 else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
133 simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
138 sprintf(errorstr, "eSimTempScale=%d not defined", simtemp->eSimTempScale);
139 gmx_fatal(FARGS, errorstr);
146 static void _low_check(gmx_bool b, char *s, warninp_t wi)
150 warning_error(wi, s);
154 static void check_nst(const char *desc_nst, int nst,
155 const char *desc_p, int *p,
160 if (*p > 0 && *p % nst != 0)
162 /* Round up to the next multiple of nst */
163 *p = ((*p)/nst + 1)*nst;
164 sprintf(buf, "%s should be a multiple of %s, changing %s to %d\n",
165 desc_p, desc_nst, desc_p, *p);
170 static gmx_bool ir_NVE(const t_inputrec *ir)
172 return ((ir->eI == eiMD || EI_VV(ir->eI)) && ir->etc == etcNO);
175 static int lcd(int n1, int n2)
180 for (i = 2; (i <= n1 && i <= n2); i++)
182 if (n1 % i == 0 && n2 % i == 0)
191 static void process_interaction_modifier(const t_inputrec *ir, int *eintmod)
193 if (*eintmod == eintmodPOTSHIFT_VERLET)
195 if (ir->cutoff_scheme == ecutsVERLET)
197 *eintmod = eintmodPOTSHIFT;
201 *eintmod = eintmodNONE;
206 void check_ir(const char *mdparin, t_inputrec *ir, t_gromppopts *opts,
208 /* Check internal consistency */
210 /* Strange macro: first one fills the err_buf, and then one can check
211 * the condition, which will print the message and increase the error
214 #define CHECK(b) _low_check(b, err_buf, wi)
215 char err_buf[256], warn_buf[STRLEN];
221 t_lambda *fep = ir->fepvals;
222 t_expanded *expand = ir->expandedvals;
224 set_warning_line(wi, mdparin, -1);
226 /* BASIC CUT-OFF STUFF */
227 if (ir->rcoulomb < 0)
229 warning_error(wi, "rcoulomb should be >= 0");
233 warning_error(wi, "rvdw should be >= 0");
236 !(ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_drift > 0))
238 warning_error(wi, "rlist should be >= 0");
241 process_interaction_modifier(ir, &ir->coulomb_modifier);
242 process_interaction_modifier(ir, &ir->vdw_modifier);
244 if (ir->cutoff_scheme == ecutsGROUP)
246 /* BASIC CUT-OFF STUFF */
247 if (ir->rlist == 0 ||
248 !((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
249 (EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist)))
251 /* No switched potential and/or no twin-range:
252 * we can set the long-range cut-off to the maximum of the other cut-offs.
254 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
256 else if (ir->rlistlong < 0)
258 ir->rlistlong = max_cutoff(ir->rlist, max_cutoff(ir->rvdw, ir->rcoulomb));
259 sprintf(warn_buf, "rlistlong was not set, setting it to %g (no buffer)",
261 warning(wi, warn_buf);
263 if (ir->rlistlong == 0 && ir->ePBC != epbcNONE)
265 warning_error(wi, "Can not have an infinite cut-off with PBC");
267 if (ir->rlistlong > 0 && (ir->rlist == 0 || ir->rlistlong < ir->rlist))
269 warning_error(wi, "rlistlong can not be shorter than rlist");
271 if (IR_TWINRANGE(*ir) && ir->nstlist <= 0)
273 warning_error(wi, "Can not have nstlist<=0 with twin-range interactions");
277 if (ir->rlistlong == ir->rlist)
281 else if (ir->rlistlong > ir->rlist && ir->nstcalclr == 0)
283 warning_error(wi, "With different cutoffs for electrostatics and VdW, nstcalclr must be -1 or a positive number");
286 if (ir->cutoff_scheme == ecutsVERLET)
290 /* Normal Verlet type neighbor-list, currently only limited feature support */
291 if (inputrec2nboundeddim(ir) < 3)
293 warning_error(wi, "With Verlet lists only full pbc or pbc=xy with walls is supported");
295 if (ir->rcoulomb != ir->rvdw)
297 warning_error(wi, "With Verlet lists rcoulomb!=rvdw is not supported");
299 if (ir->vdwtype != evdwCUT)
301 warning_error(wi, "With Verlet lists only cut-off LJ interactions are supported");
303 if (!(ir->coulombtype == eelCUT ||
304 (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC) ||
305 EEL_PME(ir->coulombtype) || ir->coulombtype == eelEWALD))
307 warning_error(wi, "With Verlet lists only cut-off, reaction-field, PME and Ewald electrostatics are supported");
310 if (ir->nstlist <= 0)
312 warning_error(wi, "With Verlet lists nstlist should be larger than 0");
315 if (ir->nstlist < 10)
317 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.");
320 rc_max = max(ir->rvdw, ir->rcoulomb);
322 if (ir->verletbuf_drift <= 0)
324 if (ir->verletbuf_drift == 0)
326 warning_error(wi, "Can not have an energy drift of exactly 0");
329 if (ir->rlist < rc_max)
331 warning_error(wi, "With verlet lists rlist can not be smaller than rvdw or rcoulomb");
334 if (ir->rlist == rc_max && ir->nstlist > 1)
336 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.");
341 if (ir->rlist > rc_max)
343 warning_note(wi, "You have set rlist larger than the interaction cut-off, but you also have verlet-buffer-drift > 0. Will set rlist using verlet-buffer-drift.");
346 if (ir->nstlist == 1)
348 /* No buffer required */
353 if (EI_DYNAMICS(ir->eI))
355 if (EI_MD(ir->eI) && ir->etc == etcNO)
357 warning_error(wi, "Temperature coupling is required for calculating rlist using the energy drift with verlet-buffer-drift > 0. Either use temperature coupling or set rlist yourself together with verlet-buffer-drift = -1.");
360 if (inputrec2nboundeddim(ir) < 3)
362 warning_error(wi, "The box volume is required for calculating rlist from the energy drift with verlet-buffer-drift > 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-drift = -1.");
364 /* Set rlist temporarily so we can continue processing */
369 /* Set the buffer to 5% of the cut-off */
370 ir->rlist = 1.05*rc_max;
375 /* No twin-range calculations with Verlet lists */
376 ir->rlistlong = ir->rlist;
379 if (ir->nstcalclr == -1)
381 /* if rlist=rlistlong, this will later be changed to nstcalclr=0 */
382 ir->nstcalclr = ir->nstlist;
384 else if (ir->nstcalclr > 0)
386 if (ir->nstlist > 0 && (ir->nstlist % ir->nstcalclr != 0))
388 warning_error(wi, "nstlist must be evenly divisible by nstcalclr. Use nstcalclr = -1 to automatically follow nstlist");
391 else if (ir->nstcalclr < -1)
393 warning_error(wi, "nstcalclr must be a positive number (divisor of nstcalclr), or -1 to follow nstlist.");
396 if (EEL_PME(ir->coulombtype) && ir->rcoulomb > ir->rvdw && ir->nstcalclr > 1)
398 warning_error(wi, "When used with PME, the long-range component of twin-range interactions must be updated every step (nstcalclr)");
401 /* GENERAL INTEGRATOR STUFF */
402 if (!(ir->eI == eiMD || EI_VV(ir->eI)))
406 if (ir->eI == eiVVAK)
408 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]);
409 warning_note(wi, warn_buf);
411 if (!EI_DYNAMICS(ir->eI))
415 if (EI_DYNAMICS(ir->eI))
417 if (ir->nstcalcenergy < 0)
419 ir->nstcalcenergy = ir_optimal_nstcalcenergy(ir);
420 if (ir->nstenergy != 0 && ir->nstenergy < ir->nstcalcenergy)
422 /* nstcalcenergy larger than nstener does not make sense.
423 * We ideally want nstcalcenergy=nstener.
427 ir->nstcalcenergy = lcd(ir->nstenergy, ir->nstlist);
431 ir->nstcalcenergy = ir->nstenergy;
435 else if ( (ir->nstenergy > 0 && ir->nstcalcenergy > ir->nstenergy) ||
436 (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
437 (ir->nstcalcenergy > ir->fepvals->nstdhdl) ) )
440 const char *nsten = "nstenergy";
441 const char *nstdh = "nstdhdl";
442 const char *min_name = nsten;
443 int min_nst = ir->nstenergy;
445 /* find the smallest of ( nstenergy, nstdhdl ) */
446 if (ir->efep != efepNO && ir->fepvals->nstdhdl > 0 &&
447 (ir->fepvals->nstdhdl < ir->nstenergy) )
449 min_nst = ir->fepvals->nstdhdl;
452 /* If the user sets nstenergy small, we should respect that */
454 "Setting nstcalcenergy (%d) equal to %s (%d)",
455 ir->nstcalcenergy, min_name, min_nst);
456 warning_note(wi, warn_buf);
457 ir->nstcalcenergy = min_nst;
460 if (ir->epc != epcNO)
462 if (ir->nstpcouple < 0)
464 ir->nstpcouple = ir_optimal_nstpcouple(ir);
467 if (IR_TWINRANGE(*ir))
469 check_nst("nstlist", ir->nstlist,
470 "nstcalcenergy", &ir->nstcalcenergy, wi);
471 if (ir->epc != epcNO)
473 check_nst("nstlist", ir->nstlist,
474 "nstpcouple", &ir->nstpcouple, wi);
478 if (ir->nstcalcenergy > 0)
480 if (ir->efep != efepNO)
482 /* nstdhdl should be a multiple of nstcalcenergy */
483 check_nst("nstcalcenergy", ir->nstcalcenergy,
484 "nstdhdl", &ir->fepvals->nstdhdl, wi);
485 /* nstexpanded should be a multiple of nstcalcenergy */
486 check_nst("nstcalcenergy", ir->nstcalcenergy,
487 "nstexpanded", &ir->expandedvals->nstexpanded, wi);
489 /* for storing exact averages nstenergy should be
490 * a multiple of nstcalcenergy
492 check_nst("nstcalcenergy", ir->nstcalcenergy,
493 "nstenergy", &ir->nstenergy, wi);
498 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
499 ir->bContinuation && ir->ld_seed != -1)
501 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)");
507 sprintf(err_buf, "TPI only works with pbc = %s", epbc_names[epbcXYZ]);
508 CHECK(ir->ePBC != epbcXYZ);
509 sprintf(err_buf, "TPI only works with ns = %s", ens_names[ensGRID]);
510 CHECK(ir->ns_type != ensGRID);
511 sprintf(err_buf, "with TPI nstlist should be larger than zero");
512 CHECK(ir->nstlist <= 0);
513 sprintf(err_buf, "TPI does not work with full electrostatics other than PME");
514 CHECK(EEL_FULL(ir->coulombtype) && !EEL_PME(ir->coulombtype));
518 if ( (opts->nshake > 0) && (opts->bMorse) )
521 "Using morse bond-potentials while constraining bonds is useless");
522 warning(wi, warn_buf);
525 if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
526 ir->bContinuation && ir->ld_seed != -1)
528 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)");
530 /* verify simulated tempering options */
534 gmx_bool bAllTempZero = TRUE;
535 for (i = 0; i < fep->n_lambda; i++)
537 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]);
538 CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
539 if (fep->all_lambda[efptTEMPERATURE][i] > 0)
541 bAllTempZero = FALSE;
544 sprintf(err_buf, "if simulated tempering is on, temperature-lambdas may not be all zero");
545 CHECK(bAllTempZero == TRUE);
547 sprintf(err_buf, "Simulated tempering is currently only compatible with md-vv");
548 CHECK(ir->eI != eiVV);
550 /* check compatability of the temperature coupling with simulated tempering */
552 if (ir->etc == etcNOSEHOOVER)
554 sprintf(warn_buf, "Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering", etcoupl_names[ir->etc]);
555 warning_note(wi, warn_buf);
558 /* check that the temperatures make sense */
560 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);
561 CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
563 sprintf(err_buf, "Higher simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_high);
564 CHECK(ir->simtempvals->simtemp_high <= 0);
566 sprintf(err_buf, "Lower simulated tempering temperature (%g) must be >= zero", ir->simtempvals->simtemp_low);
567 CHECK(ir->simtempvals->simtemp_low <= 0);
570 /* verify free energy options */
572 if (ir->efep != efepNO)
575 sprintf(err_buf, "The soft-core power is %d and can only be 1 or 2",
577 CHECK(fep->sc_alpha != 0 && fep->sc_power != 1 && fep->sc_power != 2);
579 sprintf(err_buf, "The soft-core sc-r-power is %d and can only be 6 or 48",
580 (int)fep->sc_r_power);
581 CHECK(fep->sc_alpha != 0 && fep->sc_r_power != 6.0 && fep->sc_r_power != 48.0);
583 sprintf(err_buf, "Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero", fep->delta_lambda);
584 CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state > 0) || (fep->init_lambda > 0)));
586 sprintf(err_buf, "Can't use postive delta-lambda (%g) with expanded ensemble simulations", fep->delta_lambda);
587 CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
589 sprintf(err_buf, "Can only use expanded ensemble with md-vv for now; should be supported for other integrators in 5.0");
590 CHECK(!(EI_VV(ir->eI)) && (ir->efep == efepEXPANDED));
592 sprintf(err_buf, "Free-energy not implemented for Ewald");
593 CHECK(ir->coulombtype == eelEWALD);
595 /* check validty of lambda inputs */
596 if (fep->n_lambda == 0)
598 /* Clear output in case of no states:*/
599 sprintf(err_buf, "init-lambda-state set to %d: no lambda states are defined.", fep->init_fep_state);
600 CHECK((fep->init_fep_state >= 0) && (fep->n_lambda == 0));
604 sprintf(err_buf, "initial thermodynamic state %d does not exist, only goes to %d", fep->init_fep_state, fep->n_lambda-1);
605 CHECK((fep->init_fep_state >= fep->n_lambda));
608 sprintf(err_buf, "Lambda state must be set, either with init-lambda-state or with init-lambda");
609 CHECK((fep->init_fep_state < 0) && (fep->init_lambda < 0));
611 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",
612 fep->init_lambda, fep->init_fep_state);
613 CHECK((fep->init_fep_state >= 0) && (fep->init_lambda >= 0));
617 if ((fep->init_lambda >= 0) && (fep->delta_lambda == 0))
621 for (i = 0; i < efptNR; i++)
623 if (fep->separate_dvdl[i])
628 if (n_lambda_terms > 1)
630 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.");
631 warning(wi, warn_buf);
634 if (n_lambda_terms < 2 && fep->n_lambda > 0)
637 "init-lambda is deprecated for setting lambda state (except for slow growth). Use init-lambda-state instead.");
641 for (j = 0; j < efptNR; j++)
643 for (i = 0; i < fep->n_lambda; i++)
645 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]);
646 CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
650 if ((fep->sc_alpha > 0) && (!fep->bScCoul))
652 for (i = 0; i < fep->n_lambda; i++)
654 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],
655 fep->all_lambda[efptCOUL][i]);
656 CHECK((fep->sc_alpha > 0) &&
657 (((fep->all_lambda[efptCOUL][i] > 0.0) &&
658 (fep->all_lambda[efptCOUL][i] < 1.0)) &&
659 ((fep->all_lambda[efptVDW][i] > 0.0) &&
660 (fep->all_lambda[efptVDW][i] < 1.0))));
664 if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
666 real sigma, lambda, r_sc;
669 /* Maximum estimate for A and B charges equal with lambda power 1 */
671 r_sc = pow(lambda*fep->sc_alpha*pow(sigma/ir->rcoulomb, fep->sc_r_power) + 1.0, 1.0/fep->sc_r_power);
672 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.",
674 sigma, lambda, r_sc - 1.0, ir->ewald_rtol);
675 warning_note(wi, warn_buf);
678 /* Free Energy Checks -- In an ideal world, slow growth and FEP would
679 be treated differently, but that's the next step */
681 for (i = 0; i < efptNR; i++)
683 for (j = 0; j < fep->n_lambda; j++)
685 sprintf(err_buf, "%s[%d] must be between 0 and 1", efpt_names[i], j);
686 CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
691 if ((ir->bSimTemp) || (ir->efep == efepEXPANDED))
694 expand = ir->expandedvals;
696 /* checking equilibration of weights inputs for validity */
698 sprintf(err_buf, "weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
699 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
700 CHECK((expand->equil_n_at_lam > 0) && (expand->elmceq != elmceqNUMATLAM));
702 sprintf(err_buf, "weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
703 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
704 CHECK((expand->equil_samples > 0) && (expand->elmceq != elmceqSAMPLES));
706 sprintf(err_buf, "weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
707 expand->equil_steps, elmceq_names[elmceqSTEPS]);
708 CHECK((expand->equil_steps > 0) && (expand->elmceq != elmceqSTEPS));
710 sprintf(err_buf, "weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
711 expand->equil_samples, elmceq_names[elmceqWLDELTA]);
712 CHECK((expand->equil_wl_delta > 0) && (expand->elmceq != elmceqWLDELTA));
714 sprintf(err_buf, "weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
715 expand->equil_ratio, elmceq_names[elmceqRATIO]);
716 CHECK((expand->equil_ratio > 0) && (expand->elmceq != elmceqRATIO));
718 sprintf(err_buf, "weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
719 expand->equil_n_at_lam, elmceq_names[elmceqNUMATLAM]);
720 CHECK((expand->equil_n_at_lam <= 0) && (expand->elmceq == elmceqNUMATLAM));
722 sprintf(err_buf, "weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
723 expand->equil_samples, elmceq_names[elmceqSAMPLES]);
724 CHECK((expand->equil_samples <= 0) && (expand->elmceq == elmceqSAMPLES));
726 sprintf(err_buf, "weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
727 expand->equil_steps, elmceq_names[elmceqSTEPS]);
728 CHECK((expand->equil_steps <= 0) && (expand->elmceq == elmceqSTEPS));
730 sprintf(err_buf, "weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
731 expand->equil_wl_delta, elmceq_names[elmceqWLDELTA]);
732 CHECK((expand->equil_wl_delta <= 0) && (expand->elmceq == elmceqWLDELTA));
734 sprintf(err_buf, "weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
735 expand->equil_ratio, elmceq_names[elmceqRATIO]);
736 CHECK((expand->equil_ratio <= 0) && (expand->elmceq == elmceqRATIO));
738 sprintf(err_buf, "lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
739 elmceq_names[elmceqWLDELTA], elamstats_names[elamstatsWL], elamstats_names[elamstatsWWL]);
740 CHECK((expand->elmceq == elmceqWLDELTA) && (!EWL(expand->elamstats)));
742 sprintf(err_buf, "lmc-repeats (%d) must be greater than 0", expand->lmc_repeats);
743 CHECK((expand->lmc_repeats <= 0));
744 sprintf(err_buf, "minimum-var-min (%d) must be greater than 0", expand->minvarmin);
745 CHECK((expand->minvarmin <= 0));
746 sprintf(err_buf, "weight-c-range (%d) must be greater or equal to 0", expand->c_range);
747 CHECK((expand->c_range < 0));
748 sprintf(err_buf, "init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
749 fep->init_fep_state, expand->lmc_forced_nstart);
750 CHECK((fep->init_fep_state != 0) && (expand->lmc_forced_nstart > 0) && (expand->elmcmove != elmcmoveNO));
751 sprintf(err_buf, "lmc-forced-nstart (%d) must not be negative", expand->lmc_forced_nstart);
752 CHECK((expand->lmc_forced_nstart < 0));
753 sprintf(err_buf, "init-lambda-state (%d) must be in the interval [0,number of lambdas)", fep->init_fep_state);
754 CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
756 sprintf(err_buf, "init-wl-delta (%f) must be greater than or equal to 0", expand->init_wl_delta);
757 CHECK((expand->init_wl_delta < 0));
758 sprintf(err_buf, "wl-ratio (%f) must be between 0 and 1", expand->wl_ratio);
759 CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
760 sprintf(err_buf, "wl-scale (%f) must be between 0 and 1", expand->wl_scale);
761 CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
763 /* if there is no temperature control, we need to specify an MC temperature */
764 sprintf(err_buf, "If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
765 if (expand->nstTij > 0)
767 sprintf(err_buf, "nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
768 expand->nstTij, ir->nstlog);
769 CHECK((mod(expand->nstTij, ir->nstlog) != 0));
774 sprintf(err_buf, "walls only work with pbc=%s", epbc_names[epbcXY]);
775 CHECK(ir->nwall && ir->ePBC != epbcXY);
778 if (ir->ePBC != epbcXYZ && ir->nwall != 2)
780 if (ir->ePBC == epbcNONE)
782 if (ir->epc != epcNO)
784 warning(wi, "Turning off pressure coupling for vacuum system");
790 sprintf(err_buf, "Can not have pressure coupling with pbc=%s",
791 epbc_names[ir->ePBC]);
792 CHECK(ir->epc != epcNO);
794 sprintf(err_buf, "Can not have Ewald with pbc=%s", epbc_names[ir->ePBC]);
795 CHECK(EEL_FULL(ir->coulombtype));
797 sprintf(err_buf, "Can not have dispersion correction with pbc=%s",
798 epbc_names[ir->ePBC]);
799 CHECK(ir->eDispCorr != edispcNO);
802 if (ir->rlist == 0.0)
804 sprintf(err_buf, "can only have neighborlist cut-off zero (=infinite)\n"
805 "with coulombtype = %s or coulombtype = %s\n"
806 "without periodic boundary conditions (pbc = %s) and\n"
807 "rcoulomb and rvdw set to zero",
808 eel_names[eelCUT], eel_names[eelUSER], epbc_names[epbcNONE]);
809 CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
810 (ir->ePBC != epbcNONE) ||
811 (ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
815 warning_error(wi, "Can not have heuristic neighborlist updates without cut-off");
819 warning_note(wi, "Simulating without cut-offs is usually (slightly) faster with nstlist=0, nstype=simple and particle decomposition");
824 if (ir->nstcomm == 0)
826 ir->comm_mode = ecmNO;
828 if (ir->comm_mode != ecmNO)
832 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");
833 ir->nstcomm = abs(ir->nstcomm);
836 if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy)
838 warning_note(wi, "nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
839 ir->nstcomm = ir->nstcalcenergy;
842 if (ir->comm_mode == ecmANGULAR)
844 sprintf(err_buf, "Can not remove the rotation around the center of mass with periodic molecules");
845 CHECK(ir->bPeriodicMols);
846 if (ir->ePBC != epbcNONE)
848 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).");
853 if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR)
855 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.");
858 sprintf(err_buf, "Twin-range neighbour searching (NS) with simple NS"
859 " algorithm not implemented");
860 CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
861 && (ir->ns_type == ensSIMPLE));
863 /* TEMPERATURE COUPLING */
864 if (ir->etc == etcYES)
866 ir->etc = etcBERENDSEN;
867 warning_note(wi, "Old option for temperature coupling given: "
868 "changing \"yes\" to \"Berendsen\"\n");
871 if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
873 if (ir->opts.nhchainlength < 1)
875 sprintf(warn_buf, "number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n", ir->opts.nhchainlength);
876 ir->opts.nhchainlength = 1;
877 warning(wi, warn_buf);
880 if (ir->etc == etcNOSEHOOVER && !EI_VV(ir->eI) && ir->opts.nhchainlength > 1)
882 warning_note(wi, "leapfrog does not yet support Nose-Hoover chains, nhchainlength reset to 1");
883 ir->opts.nhchainlength = 1;
888 ir->opts.nhchainlength = 0;
891 if (ir->eI == eiVVAK)
893 sprintf(err_buf, "%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
895 CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
898 if (ETC_ANDERSEN(ir->etc))
900 sprintf(err_buf, "%s temperature control not supported for integrator %s.", etcoupl_names[ir->etc], ei_names[ir->eI]);
901 CHECK(!(EI_VV(ir->eI)));
903 for (i = 0; i < ir->opts.ngtc; i++)
905 sprintf(err_buf, "all tau_t must currently be equal using Andersen temperature control, violated for group %d", i);
906 CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
907 sprintf(err_buf, "all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
908 i, ir->opts.tau_t[i]);
909 CHECK(ir->opts.tau_t[i] < 0);
911 if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN))
913 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]);
914 warning_note(wi, warn_buf);
917 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]);
918 CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
920 for (i = 0; i < ir->opts.ngtc; i++)
922 int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
923 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);
924 CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
927 if (ir->etc == etcBERENDSEN)
929 sprintf(warn_buf, "The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
930 ETCOUPLTYPE(ir->etc), ETCOUPLTYPE(etcVRESCALE));
931 warning_note(wi, warn_buf);
934 if ((ir->etc == etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
935 && ir->epc == epcBERENDSEN)
937 sprintf(warn_buf, "Using Berendsen pressure coupling invalidates the "
938 "true ensemble for the thermostat");
939 warning(wi, warn_buf);
942 /* PRESSURE COUPLING */
943 if (ir->epc == epcISOTROPIC)
945 ir->epc = epcBERENDSEN;
946 warning_note(wi, "Old option for pressure coupling given: "
947 "changing \"Isotropic\" to \"Berendsen\"\n");
950 if (ir->epc != epcNO)
952 dt_pcoupl = ir->nstpcouple*ir->delta_t;
954 sprintf(err_buf, "tau-p must be > 0 instead of %g\n", ir->tau_p);
955 CHECK(ir->tau_p <= 0);
957 if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
959 sprintf(warn_buf, "For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
960 EPCOUPLTYPE(ir->epc), ir->tau_p, pcouple_min_integration_steps(ir->epc), dt_pcoupl);
961 warning(wi, warn_buf);
964 sprintf(err_buf, "compressibility must be > 0 when using pressure"
965 " coupling %s\n", EPCOUPLTYPE(ir->epc));
966 CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
967 ir->compress[ZZ][ZZ] < 0 ||
968 (trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
969 ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
971 if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
974 "You are generating velocities so I am assuming you "
975 "are equilibrating a system. You are using "
976 "%s pressure coupling, but this can be "
977 "unstable for equilibration. If your system crashes, try "
978 "equilibrating first with Berendsen pressure coupling. If "
979 "you are not equilibrating the system, you can probably "
980 "ignore this warning.",
981 epcoupl_names[ir->epc]);
982 warning(wi, warn_buf);
990 if ((ir->epc != epcBERENDSEN) && (ir->epc != epcMTTK))
992 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.");
998 if (ir->epc == epcMTTK)
1000 warning_error(wi, "MTTK pressure coupling requires a Velocity-verlet integrator");
1004 /* ELECTROSTATICS */
1005 /* More checks are in triple check (grompp.c) */
1007 if (ir->coulombtype == eelSWITCH)
1009 sprintf(warn_buf, "coulombtype = %s is only for testing purposes and can lead to serious "
1010 "artifacts, advice: use coulombtype = %s",
1011 eel_names[ir->coulombtype],
1012 eel_names[eelRF_ZERO]);
1013 warning(wi, warn_buf);
1016 if (ir->epsilon_r != 1 && ir->implicit_solvent == eisGBSA)
1018 sprintf(warn_buf, "epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric", ir->epsilon_r);
1019 warning_note(wi, warn_buf);
1022 if (EEL_RF(ir->coulombtype) && ir->epsilon_rf == 1 && ir->epsilon_r != 1)
1024 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);
1025 warning(wi, warn_buf);
1026 ir->epsilon_rf = ir->epsilon_r;
1027 ir->epsilon_r = 1.0;
1030 if (getenv("GALACTIC_DYNAMICS") == NULL)
1032 sprintf(err_buf, "epsilon-r must be >= 0 instead of %g\n", ir->epsilon_r);
1033 CHECK(ir->epsilon_r < 0);
1036 if (EEL_RF(ir->coulombtype))
1038 /* reaction field (at the cut-off) */
1040 if (ir->coulombtype == eelRF_ZERO)
1042 sprintf(warn_buf, "With coulombtype = %s, epsilon-rf must be 0, assuming you meant epsilon_rf=0",
1043 eel_names[ir->coulombtype]);
1044 CHECK(ir->epsilon_rf != 0);
1045 ir->epsilon_rf = 0.0;
1048 sprintf(err_buf, "epsilon-rf must be >= epsilon-r");
1049 CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
1050 (ir->epsilon_r == 0));
1051 if (ir->epsilon_rf == ir->epsilon_r)
1053 sprintf(warn_buf, "Using epsilon-rf = epsilon-r with %s does not make sense",
1054 eel_names[ir->coulombtype]);
1055 warning(wi, warn_buf);
1058 /* Allow rlist>rcoulomb for tabulated long range stuff. This just
1059 * means the interaction is zero outside rcoulomb, but it helps to
1060 * provide accurate energy conservation.
1062 if (EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype))
1064 if (EEL_SWITCHED(ir->coulombtype))
1067 "With coulombtype = %s rcoulomb_switch must be < rcoulomb. Or, better: Use the potential modifier options!",
1068 eel_names[ir->coulombtype]);
1069 CHECK(ir->rcoulomb_switch >= ir->rcoulomb);
1072 else if (ir->coulombtype == eelCUT || EEL_RF(ir->coulombtype))
1074 if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1076 sprintf(err_buf, "With coulombtype = %s, rcoulomb should be >= rlist unless you use a potential modifier",
1077 eel_names[ir->coulombtype]);
1078 CHECK(ir->rlist > ir->rcoulomb);
1082 if (ir->coulombtype == eelSWITCH || ir->coulombtype == eelSHIFT ||
1083 ir->vdwtype == evdwSWITCH || ir->vdwtype == evdwSHIFT)
1086 "The switch/shift interaction settings are just for compatibility; you will get better "
1087 "performance from applying potential modifiers to your interactions!\n");
1088 warning_note(wi, warn_buf);
1091 if (ir->coulombtype == eelPMESWITCH || ir->coulomb_modifier == eintmodPOTSWITCH)
1093 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1095 real percentage = 100*(ir->rcoulomb-ir->rcoulomb_switch)/ir->rcoulomb;
1096 sprintf(warn_buf, "The switching range for 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.",
1097 percentage,ir->rcoulomb_switch,ir->rcoulomb,ir->ewald_rtol);
1098 warning(wi, warn_buf);
1102 if (ir->vdwtype == evdwSWITCH || ir->vdw_modifier == eintmodPOTSWITCH)
1104 if (ir->rvdw_switch==0)
1106 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.");
1107 warning(wi, warn_buf);
1111 if (EEL_FULL(ir->coulombtype))
1113 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1114 ir->coulombtype == eelPMEUSERSWITCH)
1116 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1117 eel_names[ir->coulombtype]);
1118 CHECK(ir->rcoulomb > ir->rlist);
1120 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1122 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1125 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1126 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1127 "a potential modifier.", eel_names[ir->coulombtype]);
1128 if (ir->nstcalclr == 1)
1130 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1134 CHECK(ir->rcoulomb != ir->rlist);
1140 if (EEL_PME(ir->coulombtype))
1142 if (ir->pme_order < 3)
1144 warning_error(wi, "pme-order can not be smaller than 3");
1148 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1150 if (ir->ewald_geometry == eewg3D)
1152 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1153 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1154 warning(wi, warn_buf);
1156 /* This check avoids extra pbc coding for exclusion corrections */
1157 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1158 CHECK(ir->wall_ewald_zfac < 2);
1161 if (EVDW_SWITCHED(ir->vdwtype))
1163 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1164 evdw_names[ir->vdwtype]);
1165 CHECK(ir->rvdw_switch >= ir->rvdw);
1167 else if (ir->vdwtype == evdwCUT)
1169 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1171 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1172 CHECK(ir->rlist > ir->rvdw);
1175 if (ir->cutoff_scheme == ecutsGROUP)
1177 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1178 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1181 warning_note(wi, "With exact cut-offs, rlist should be "
1182 "larger than rcoulomb and rvdw, so that there "
1183 "is a buffer region for particle motion "
1184 "between neighborsearch steps");
1187 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1188 && (ir->rlistlong <= ir->rcoulomb))
1190 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1191 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1192 warning_note(wi, warn_buf);
1194 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1196 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1197 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1198 warning_note(wi, warn_buf);
1202 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1204 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.");
1207 if (ir->nstlist == -1)
1209 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1210 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1212 sprintf(err_buf, "nstlist can not be smaller than -1");
1213 CHECK(ir->nstlist < -1);
1215 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1218 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1221 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1223 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1226 /* ENERGY CONSERVATION */
1227 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1229 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1231 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1232 evdw_names[evdwSHIFT]);
1233 warning_note(wi, warn_buf);
1235 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1237 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1238 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1239 warning_note(wi, warn_buf);
1243 if (EI_VV(ir->eI) && IR_TWINRANGE(*ir) && ir->nstlist > 1)
1245 sprintf(warn_buf, "Twin-range multiple time stepping does not work with integrator %s.", ei_names[ir->eI]);
1246 warning_error(wi, warn_buf);
1249 /* IMPLICIT SOLVENT */
1250 if (ir->coulombtype == eelGB_NOTUSED)
1252 ir->coulombtype = eelCUT;
1253 ir->implicit_solvent = eisGBSA;
1254 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1255 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1256 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1259 if (ir->sa_algorithm == esaSTILL)
1261 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1262 CHECK(ir->sa_algorithm == esaSTILL);
1265 if (ir->implicit_solvent == eisGBSA)
1267 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1268 CHECK(ir->rgbradii != ir->rlist);
1270 if (ir->coulombtype != eelCUT)
1272 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1273 CHECK(ir->coulombtype != eelCUT);
1275 if (ir->vdwtype != evdwCUT)
1277 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1278 CHECK(ir->vdwtype != evdwCUT);
1280 if (ir->nstgbradii < 1)
1282 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1283 warning_note(wi, warn_buf);
1286 if (ir->sa_algorithm == esaNO)
1288 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1289 warning_note(wi, warn_buf);
1291 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1293 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");
1294 warning_note(wi, warn_buf);
1296 if (ir->gb_algorithm == egbSTILL)
1298 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1302 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1305 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1307 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1308 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1315 if (ir->cutoff_scheme != ecutsGROUP)
1317 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1321 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1323 if (ir->epc != epcNO)
1325 warning_error(wi, "AdresS simulation does not support pressure coupling");
1327 if (EEL_FULL(ir->coulombtype))
1329 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1334 /* count the number of text elemets separated by whitespace in a string.
1335 str = the input string
1336 maxptr = the maximum number of allowed elements
1337 ptr = the output array of pointers to the first character of each element
1338 returns: the number of elements. */
1339 int str_nelem(const char *str, int maxptr, char *ptr[])
1344 copy0 = strdup(str);
1347 while (*copy != '\0')
1351 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1359 while ((*copy != '\0') && !isspace(*copy))
1378 /* interpret a number of doubles from a string and put them in an array,
1379 after allocating space for them.
1380 str = the input string
1381 n = the (pre-allocated) number of doubles read
1382 r = the output array of doubles. */
1383 static void parse_n_real(char *str, int *n, real **r)
1388 *n = str_nelem(str, MAXPTR, ptr);
1391 for (i = 0; i < *n; i++)
1393 (*r)[i] = strtod(ptr[i], NULL);
1397 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1400 int i, j, max_n_lambda, nweights, nfep[efptNR];
1401 t_lambda *fep = ir->fepvals;
1402 t_expanded *expand = ir->expandedvals;
1403 real **count_fep_lambdas;
1404 gmx_bool bOneLambda = TRUE;
1406 snew(count_fep_lambdas, efptNR);
1408 /* FEP input processing */
1409 /* first, identify the number of lambda values for each type.
1410 All that are nonzero must have the same number */
1412 for (i = 0; i < efptNR; i++)
1414 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1417 /* now, determine the number of components. All must be either zero, or equal. */
1420 for (i = 0; i < efptNR; i++)
1422 if (nfep[i] > max_n_lambda)
1424 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1425 must have the same number if its not zero.*/
1430 for (i = 0; i < efptNR; i++)
1434 ir->fepvals->separate_dvdl[i] = FALSE;
1436 else if (nfep[i] == max_n_lambda)
1438 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1439 respect to the temperature currently */
1441 ir->fepvals->separate_dvdl[i] = TRUE;
1446 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1447 nfep[i], efpt_names[i], max_n_lambda);
1450 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1451 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1453 /* the number of lambdas is the number we've read in, which is either zero
1454 or the same for all */
1455 fep->n_lambda = max_n_lambda;
1457 /* allocate space for the array of lambda values */
1458 snew(fep->all_lambda, efptNR);
1459 /* if init_lambda is defined, we need to set lambda */
1460 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1462 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1464 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1465 for (i = 0; i < efptNR; i++)
1467 snew(fep->all_lambda[i], fep->n_lambda);
1468 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1471 for (j = 0; j < fep->n_lambda; j++)
1473 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1475 sfree(count_fep_lambdas[i]);
1478 sfree(count_fep_lambdas);
1480 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1481 bookkeeping -- for now, init_lambda */
1483 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1485 for (i = 0; i < fep->n_lambda; i++)
1487 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1491 /* check to see if only a single component lambda is defined, and soft core is defined.
1492 In this case, turn on coulomb soft core */
1494 if (max_n_lambda == 0)
1500 for (i = 0; i < efptNR; i++)
1502 if ((nfep[i] != 0) && (i != efptFEP))
1508 if ((bOneLambda) && (fep->sc_alpha > 0))
1510 fep->bScCoul = TRUE;
1513 /* Fill in the others with the efptFEP if they are not explicitly
1514 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1515 they are all zero. */
1517 for (i = 0; i < efptNR; i++)
1519 if ((nfep[i] == 0) && (i != efptFEP))
1521 for (j = 0; j < fep->n_lambda; j++)
1523 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1529 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1530 if (fep->sc_r_power == 48)
1532 if (fep->sc_alpha > 0.1)
1534 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1538 expand = ir->expandedvals;
1539 /* now read in the weights */
1540 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1543 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1545 else if (nweights != fep->n_lambda)
1547 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1548 nweights, fep->n_lambda);
1550 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1552 expand->nstexpanded = fep->nstdhdl;
1553 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1555 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1557 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1558 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1559 2*tau_t just to be careful so it's not to frequent */
1564 static void do_simtemp_params(t_inputrec *ir)
1567 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1568 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1573 static void do_wall_params(t_inputrec *ir,
1574 char *wall_atomtype, char *wall_density,
1578 char *names[MAXPTR];
1581 opts->wall_atomtype[0] = NULL;
1582 opts->wall_atomtype[1] = NULL;
1584 ir->wall_atomtype[0] = -1;
1585 ir->wall_atomtype[1] = -1;
1586 ir->wall_density[0] = 0;
1587 ir->wall_density[1] = 0;
1591 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1592 if (nstr != ir->nwall)
1594 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1597 for (i = 0; i < ir->nwall; i++)
1599 opts->wall_atomtype[i] = strdup(names[i]);
1602 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1604 nstr = str_nelem(wall_density, MAXPTR, names);
1605 if (nstr != ir->nwall)
1607 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1609 for (i = 0; i < ir->nwall; i++)
1611 sscanf(names[i], "%lf", &dbl);
1614 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1616 ir->wall_density[i] = dbl;
1622 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1630 srenew(groups->grpname, groups->ngrpname+nwall);
1631 grps = &(groups->grps[egcENER]);
1632 srenew(grps->nm_ind, grps->nr+nwall);
1633 for (i = 0; i < nwall; i++)
1635 sprintf(str, "wall%d", i);
1636 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1637 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1642 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1643 t_expanded *expand, warninp_t wi)
1645 int ninp, nerror = 0;
1651 /* read expanded ensemble parameters */
1652 CCTYPE ("expanded ensemble variables");
1653 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1654 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1655 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1656 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1657 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1658 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1659 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1660 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1661 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1662 CCTYPE("Seed for Monte Carlo in lambda space");
1663 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1664 RTYPE ("mc-temperature", expand->mc_temp, -1);
1665 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1666 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1667 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1668 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1669 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1670 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1671 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1672 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1673 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1674 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1675 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1683 void get_ir(const char *mdparin, const char *mdparout,
1684 t_inputrec *ir, t_gromppopts *opts,
1688 double dumdub[2][6];
1692 char warn_buf[STRLEN];
1693 t_lambda *fep = ir->fepvals;
1694 t_expanded *expand = ir->expandedvals;
1696 inp = read_inpfile(mdparin, &ninp, NULL, wi);
1698 snew(dumstr[0], STRLEN);
1699 snew(dumstr[1], STRLEN);
1701 /* remove the following deprecated commands */
1704 REM_TYPE("domain-decomposition");
1705 REM_TYPE("andersen-seed");
1707 REM_TYPE("dihre-fc");
1708 REM_TYPE("dihre-tau");
1709 REM_TYPE("nstdihreout");
1710 REM_TYPE("nstcheckpoint");
1712 /* replace the following commands with the clearer new versions*/
1713 REPL_TYPE("unconstrained-start", "continuation");
1714 REPL_TYPE("foreign-lambda", "fep-lambdas");
1716 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1717 CTYPE ("Preprocessor information: use cpp syntax.");
1718 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1719 STYPE ("include", opts->include, NULL);
1720 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1721 STYPE ("define", opts->define, NULL);
1723 CCTYPE ("RUN CONTROL PARAMETERS");
1724 EETYPE("integrator", ir->eI, ei_names);
1725 CTYPE ("Start time and timestep in ps");
1726 RTYPE ("tinit", ir->init_t, 0.0);
1727 RTYPE ("dt", ir->delta_t, 0.001);
1728 STEPTYPE ("nsteps", ir->nsteps, 0);
1729 CTYPE ("For exact run continuation or redoing part of a run");
1730 STEPTYPE ("init-step", ir->init_step, 0);
1731 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1732 ITYPE ("simulation-part", ir->simulation_part, 1);
1733 CTYPE ("mode for center of mass motion removal");
1734 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1735 CTYPE ("number of steps for center of mass motion removal");
1736 ITYPE ("nstcomm", ir->nstcomm, 100);
1737 CTYPE ("group(s) for center of mass motion removal");
1738 STYPE ("comm-grps", vcm, NULL);
1740 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1741 CTYPE ("Friction coefficient (amu/ps) and random seed");
1742 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1743 ITYPE ("ld-seed", ir->ld_seed, 1993);
1746 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1747 CTYPE ("Force tolerance and initial step-size");
1748 RTYPE ("emtol", ir->em_tol, 10.0);
1749 RTYPE ("emstep", ir->em_stepsize, 0.01);
1750 CTYPE ("Max number of iterations in relax-shells");
1751 ITYPE ("niter", ir->niter, 20);
1752 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1753 RTYPE ("fcstep", ir->fc_stepsize, 0);
1754 CTYPE ("Frequency of steepest descents steps when doing CG");
1755 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1756 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1758 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1759 RTYPE ("rtpi", ir->rtpi, 0.05);
1761 /* Output options */
1762 CCTYPE ("OUTPUT CONTROL OPTIONS");
1763 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1764 ITYPE ("nstxout", ir->nstxout, 0);
1765 ITYPE ("nstvout", ir->nstvout, 0);
1766 ITYPE ("nstfout", ir->nstfout, 0);
1767 ir->nstcheckpoint = 1000;
1768 CTYPE ("Output frequency for energies to log file and energy file");
1769 ITYPE ("nstlog", ir->nstlog, 1000);
1770 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1771 ITYPE ("nstenergy", ir->nstenergy, 1000);
1772 CTYPE ("Output frequency and precision for .xtc file");
1773 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1774 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1775 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1776 CTYPE ("select multiple groups. By default all atoms will be written.");
1777 STYPE ("xtc-grps", xtc_grps, NULL);
1778 CTYPE ("Selection of energy groups");
1779 STYPE ("energygrps", energy, NULL);
1781 /* Neighbor searching */
1782 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1783 CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1784 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1785 CTYPE ("nblist update frequency");
1786 ITYPE ("nstlist", ir->nstlist, 10);
1787 CTYPE ("ns algorithm (simple or grid)");
1788 EETYPE("ns-type", ir->ns_type, ens_names);
1789 /* set ndelta to the optimal value of 2 */
1791 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1792 EETYPE("pbc", ir->ePBC, epbc_names);
1793 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1794 CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1795 CTYPE ("a value of -1 means: use rlist");
1796 RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
1797 CTYPE ("nblist cut-off");
1798 RTYPE ("rlist", ir->rlist, 1.0);
1799 CTYPE ("long-range cut-off for switched potentials");
1800 RTYPE ("rlistlong", ir->rlistlong, -1);
1801 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1803 /* Electrostatics */
1804 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1805 CTYPE ("Method for doing electrostatics");
1806 EETYPE("coulombtype", ir->coulombtype, eel_names);
1807 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1808 CTYPE ("cut-off lengths");
1809 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1810 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1811 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1812 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1813 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1814 CTYPE ("Method for doing Van der Waals");
1815 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1816 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1817 CTYPE ("cut-off lengths");
1818 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1819 RTYPE ("rvdw", ir->rvdw, 1.0);
1820 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1821 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1822 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1823 RTYPE ("table-extension", ir->tabext, 1.0);
1824 CTYPE ("Separate tables between energy group pairs");
1825 STYPE ("energygrp-table", egptable, NULL);
1826 CTYPE ("Spacing for the PME/PPPM FFT grid");
1827 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1828 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1829 ITYPE ("fourier-nx", ir->nkx, 0);
1830 ITYPE ("fourier-ny", ir->nky, 0);
1831 ITYPE ("fourier-nz", ir->nkz, 0);
1832 CTYPE ("EWALD/PME/PPPM parameters");
1833 ITYPE ("pme-order", ir->pme_order, 4);
1834 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1835 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1836 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1837 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1839 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1840 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1842 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1843 CTYPE ("Algorithm for calculating Born radii");
1844 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1845 CTYPE ("Frequency of calculating the Born radii inside rlist");
1846 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1847 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1848 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1849 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1850 CTYPE ("Dielectric coefficient of the implicit solvent");
1851 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1852 CTYPE ("Salt concentration in M for Generalized Born models");
1853 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1854 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1855 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1856 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1857 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1858 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1859 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1860 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1861 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1862 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1864 /* Coupling stuff */
1865 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1866 CTYPE ("Temperature coupling");
1867 EETYPE("tcoupl", ir->etc, etcoupl_names);
1868 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1869 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1870 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1871 CTYPE ("Groups to couple separately");
1872 STYPE ("tc-grps", tcgrps, NULL);
1873 CTYPE ("Time constant (ps) and reference temperature (K)");
1874 STYPE ("tau-t", tau_t, NULL);
1875 STYPE ("ref-t", ref_t, NULL);
1876 CTYPE ("pressure coupling");
1877 EETYPE("pcoupl", ir->epc, epcoupl_names);
1878 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1879 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1880 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1881 RTYPE ("tau-p", ir->tau_p, 1.0);
1882 STYPE ("compressibility", dumstr[0], NULL);
1883 STYPE ("ref-p", dumstr[1], NULL);
1884 CTYPE ("Scaling of reference coordinates, No, All or COM");
1885 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1888 CCTYPE ("OPTIONS FOR QMMM calculations");
1889 EETYPE("QMMM", ir->bQMMM, yesno_names);
1890 CTYPE ("Groups treated Quantum Mechanically");
1891 STYPE ("QMMM-grps", QMMM, NULL);
1892 CTYPE ("QM method");
1893 STYPE("QMmethod", QMmethod, NULL);
1894 CTYPE ("QMMM scheme");
1895 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1896 CTYPE ("QM basisset");
1897 STYPE("QMbasis", QMbasis, NULL);
1898 CTYPE ("QM charge");
1899 STYPE ("QMcharge", QMcharge, NULL);
1900 CTYPE ("QM multiplicity");
1901 STYPE ("QMmult", QMmult, NULL);
1902 CTYPE ("Surface Hopping");
1903 STYPE ("SH", bSH, NULL);
1904 CTYPE ("CAS space options");
1905 STYPE ("CASorbitals", CASorbitals, NULL);
1906 STYPE ("CASelectrons", CASelectrons, NULL);
1907 STYPE ("SAon", SAon, NULL);
1908 STYPE ("SAoff", SAoff, NULL);
1909 STYPE ("SAsteps", SAsteps, NULL);
1910 CTYPE ("Scale factor for MM charges");
1911 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1912 CTYPE ("Optimization of QM subsystem");
1913 STYPE ("bOPT", bOPT, NULL);
1914 STYPE ("bTS", bTS, NULL);
1916 /* Simulated annealing */
1917 CCTYPE("SIMULATED ANNEALING");
1918 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1919 STYPE ("annealing", anneal, NULL);
1920 CTYPE ("Number of time points to use for specifying annealing in each group");
1921 STYPE ("annealing-npoints", anneal_npoints, NULL);
1922 CTYPE ("List of times at the annealing points for each group");
1923 STYPE ("annealing-time", anneal_time, NULL);
1924 CTYPE ("Temp. at each annealing point, for each group.");
1925 STYPE ("annealing-temp", anneal_temp, NULL);
1928 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1929 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1930 RTYPE ("gen-temp", opts->tempi, 300.0);
1931 ITYPE ("gen-seed", opts->seed, 173529);
1934 CCTYPE ("OPTIONS FOR BONDS");
1935 EETYPE("constraints", opts->nshake, constraints);
1936 CTYPE ("Type of constraint algorithm");
1937 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1938 CTYPE ("Do not constrain the start configuration");
1939 EETYPE("continuation", ir->bContinuation, yesno_names);
1940 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1941 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1942 CTYPE ("Relative tolerance of shake");
1943 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1944 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1945 ITYPE ("lincs-order", ir->nProjOrder, 4);
1946 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1947 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1948 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1949 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1950 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1951 CTYPE ("rotates over more degrees than");
1952 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1953 CTYPE ("Convert harmonic bonds to morse potentials");
1954 EETYPE("morse", opts->bMorse, yesno_names);
1956 /* Energy group exclusions */
1957 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1958 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1959 STYPE ("energygrp-excl", egpexcl, NULL);
1963 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1964 ITYPE ("nwall", ir->nwall, 0);
1965 EETYPE("wall-type", ir->wall_type, ewt_names);
1966 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1967 STYPE ("wall-atomtype", wall_atomtype, NULL);
1968 STYPE ("wall-density", wall_density, NULL);
1969 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1972 CCTYPE("COM PULLING");
1973 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1974 EETYPE("pull", ir->ePull, epull_names);
1975 if (ir->ePull != epullNO)
1978 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1981 /* Enforced rotation */
1982 CCTYPE("ENFORCED ROTATION");
1983 CTYPE("Enforced rotation: No or Yes");
1984 EETYPE("rotation", ir->bRot, yesno_names);
1988 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1992 CCTYPE("NMR refinement stuff");
1993 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1994 EETYPE("disre", ir->eDisre, edisre_names);
1995 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1996 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1997 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1998 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1999 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
2000 RTYPE ("disre-tau", ir->dr_tau, 0.0);
2001 CTYPE ("Output frequency for pair distances to energy file");
2002 ITYPE ("nstdisreout", ir->nstdisreout, 100);
2003 CTYPE ("Orientation restraints: No or Yes");
2004 EETYPE("orire", opts->bOrire, yesno_names);
2005 CTYPE ("Orientation restraints force constant and tau for time averaging");
2006 RTYPE ("orire-fc", ir->orires_fc, 0.0);
2007 RTYPE ("orire-tau", ir->orires_tau, 0.0);
2008 STYPE ("orire-fitgrp", orirefitgrp, NULL);
2009 CTYPE ("Output frequency for trace(SD) and S to energy file");
2010 ITYPE ("nstorireout", ir->nstorireout, 100);
2012 /* free energy variables */
2013 CCTYPE ("Free energy variables");
2014 EETYPE("free-energy", ir->efep, efep_names);
2015 STYPE ("couple-moltype", couple_moltype, NULL);
2016 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2017 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2018 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2020 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2022 it was not entered */
2023 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2024 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2025 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2026 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2027 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2028 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2029 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2030 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2031 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2032 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2033 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2034 STYPE ("init-lambda-weights", lambda_weights, NULL);
2035 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2036 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2037 ITYPE ("sc-power", fep->sc_power, 1);
2038 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2039 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2040 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2041 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2042 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2043 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2044 separate_dhdl_file_names);
2045 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2046 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2047 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2049 /* Non-equilibrium MD stuff */
2050 CCTYPE("Non-equilibrium MD stuff");
2051 STYPE ("acc-grps", accgrps, NULL);
2052 STYPE ("accelerate", acc, NULL);
2053 STYPE ("freezegrps", freeze, NULL);
2054 STYPE ("freezedim", frdim, NULL);
2055 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2056 STYPE ("deform", deform, NULL);
2058 /* simulated tempering variables */
2059 CCTYPE("simulated tempering variables");
2060 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2061 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2062 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2063 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2065 /* expanded ensemble variables */
2066 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2068 read_expandedparams(&ninp, &inp, expand, wi);
2071 /* Electric fields */
2072 CCTYPE("Electric fields");
2073 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2074 CTYPE ("and a phase angle (real)");
2075 STYPE ("E-x", efield_x, NULL);
2076 STYPE ("E-xt", efield_xt, NULL);
2077 STYPE ("E-y", efield_y, NULL);
2078 STYPE ("E-yt", efield_yt, NULL);
2079 STYPE ("E-z", efield_z, NULL);
2080 STYPE ("E-zt", efield_zt, NULL);
2082 /* AdResS defined thingies */
2083 CCTYPE ("AdResS parameters");
2084 EETYPE("adress", ir->bAdress, yesno_names);
2087 snew(ir->adress, 1);
2088 read_adressparams(&ninp, &inp, ir->adress, wi);
2091 /* User defined thingies */
2092 CCTYPE ("User defined thingies");
2093 STYPE ("user1-grps", user1, NULL);
2094 STYPE ("user2-grps", user2, NULL);
2095 ITYPE ("userint1", ir->userint1, 0);
2096 ITYPE ("userint2", ir->userint2, 0);
2097 ITYPE ("userint3", ir->userint3, 0);
2098 ITYPE ("userint4", ir->userint4, 0);
2099 RTYPE ("userreal1", ir->userreal1, 0);
2100 RTYPE ("userreal2", ir->userreal2, 0);
2101 RTYPE ("userreal3", ir->userreal3, 0);
2102 RTYPE ("userreal4", ir->userreal4, 0);
2105 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2106 for (i = 0; (i < ninp); i++)
2109 sfree(inp[i].value);
2113 /* Process options if necessary */
2114 for (m = 0; m < 2; m++)
2116 for (i = 0; i < 2*DIM; i++)
2125 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2127 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2129 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2131 case epctSEMIISOTROPIC:
2132 case epctSURFACETENSION:
2133 if (sscanf(dumstr[m], "%lf%lf",
2134 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2136 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2138 dumdub[m][YY] = dumdub[m][XX];
2140 case epctANISOTROPIC:
2141 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2142 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2143 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2145 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2149 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2150 epcoupltype_names[ir->epct]);
2154 clear_mat(ir->ref_p);
2155 clear_mat(ir->compress);
2156 for (i = 0; i < DIM; i++)
2158 ir->ref_p[i][i] = dumdub[1][i];
2159 ir->compress[i][i] = dumdub[0][i];
2161 if (ir->epct == epctANISOTROPIC)
2163 ir->ref_p[XX][YY] = dumdub[1][3];
2164 ir->ref_p[XX][ZZ] = dumdub[1][4];
2165 ir->ref_p[YY][ZZ] = dumdub[1][5];
2166 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2168 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2170 ir->compress[XX][YY] = dumdub[0][3];
2171 ir->compress[XX][ZZ] = dumdub[0][4];
2172 ir->compress[YY][ZZ] = dumdub[0][5];
2173 for (i = 0; i < DIM; i++)
2175 for (m = 0; m < i; m++)
2177 ir->ref_p[i][m] = ir->ref_p[m][i];
2178 ir->compress[i][m] = ir->compress[m][i];
2183 if (ir->comm_mode == ecmNO)
2188 opts->couple_moltype = NULL;
2189 if (strlen(couple_moltype) > 0)
2191 if (ir->efep != efepNO)
2193 opts->couple_moltype = strdup(couple_moltype);
2194 if (opts->couple_lam0 == opts->couple_lam1)
2196 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2198 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2199 opts->couple_lam1 == ecouplamNONE))
2201 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2206 warning(wi, "Can not couple a molecule with free_energy = no");
2209 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2210 if (ir->efep != efepNO)
2212 if (fep->delta_lambda > 0)
2214 ir->efep = efepSLOWGROWTH;
2220 fep->bPrintEnergy = TRUE;
2221 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2222 if the temperature is changing. */
2225 if ((ir->efep != efepNO) || ir->bSimTemp)
2227 ir->bExpanded = FALSE;
2228 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2230 ir->bExpanded = TRUE;
2232 do_fep_params(ir, fep_lambda, lambda_weights);
2233 if (ir->bSimTemp) /* done after fep params */
2235 do_simtemp_params(ir);
2240 ir->fepvals->n_lambda = 0;
2243 /* WALL PARAMETERS */
2245 do_wall_params(ir, wall_atomtype, wall_density, opts);
2247 /* ORIENTATION RESTRAINT PARAMETERS */
2249 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2251 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2254 /* DEFORMATION PARAMETERS */
2256 clear_mat(ir->deform);
2257 for (i = 0; i < 6; i++)
2261 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2262 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2263 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2264 for (i = 0; i < 3; i++)
2266 ir->deform[i][i] = dumdub[0][i];
2268 ir->deform[YY][XX] = dumdub[0][3];
2269 ir->deform[ZZ][XX] = dumdub[0][4];
2270 ir->deform[ZZ][YY] = dumdub[0][5];
2271 if (ir->epc != epcNO)
2273 for (i = 0; i < 3; i++)
2275 for (j = 0; j <= i; j++)
2277 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2279 warning_error(wi, "A box element has deform set and compressibility > 0");
2283 for (i = 0; i < 3; i++)
2285 for (j = 0; j < i; j++)
2287 if (ir->deform[i][j] != 0)
2289 for (m = j; m < DIM; m++)
2291 if (ir->compress[m][j] != 0)
2293 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.");
2294 warning(wi, warn_buf);
2306 static int search_QMstring(char *s, int ng, const char *gn[])
2308 /* same as normal search_string, but this one searches QM strings */
2311 for (i = 0; (i < ng); i++)
2313 if (gmx_strcasecmp(s, gn[i]) == 0)
2319 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2323 } /* search_QMstring */
2326 int search_string(char *s, int ng, char *gn[])
2330 for (i = 0; (i < ng); i++)
2332 if (gmx_strcasecmp(s, gn[i]) == 0)
2339 "Group %s referenced in the .mdp file was not found in the index file.\n"
2340 "Group names must match either [moleculetype] names or custom index group\n"
2341 "names, in which case you must supply an index file to the '-n' option\n"
2348 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2349 t_blocka *block, char *gnames[],
2350 int gtype, int restnm,
2351 int grptp, gmx_bool bVerbose,
2354 unsigned short *cbuf;
2355 t_grps *grps = &(groups->grps[gtype]);
2356 int i, j, gid, aj, ognr, ntot = 0;
2359 char warn_buf[STRLEN];
2363 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2366 title = gtypes[gtype];
2369 /* Mark all id's as not set */
2370 for (i = 0; (i < natoms); i++)
2375 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2376 for (i = 0; (i < ng); i++)
2378 /* Lookup the group name in the block structure */
2379 gid = search_string(ptrs[i], block->nr, gnames);
2380 if ((grptp != egrptpONE) || (i == 0))
2382 grps->nm_ind[grps->nr++] = gid;
2386 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2389 /* Now go over the atoms in the group */
2390 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2395 /* Range checking */
2396 if ((aj < 0) || (aj >= natoms))
2398 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2400 /* Lookup up the old group number */
2404 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2405 aj+1, title, ognr+1, i+1);
2409 /* Store the group number in buffer */
2410 if (grptp == egrptpONE)
2423 /* Now check whether we have done all atoms */
2427 if (grptp == egrptpALL)
2429 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2430 natoms-ntot, title);
2432 else if (grptp == egrptpPART)
2434 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2435 natoms-ntot, title);
2436 warning_note(wi, warn_buf);
2438 /* Assign all atoms currently unassigned to a rest group */
2439 for (j = 0; (j < natoms); j++)
2441 if (cbuf[j] == NOGID)
2447 if (grptp != egrptpPART)
2452 "Making dummy/rest group for %s containing %d elements\n",
2453 title, natoms-ntot);
2455 /* Add group name "rest" */
2456 grps->nm_ind[grps->nr] = restnm;
2458 /* Assign the rest name to all atoms not currently assigned to a group */
2459 for (j = 0; (j < natoms); j++)
2461 if (cbuf[j] == NOGID)
2470 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2472 /* All atoms are part of one (or no) group, no index required */
2473 groups->ngrpnr[gtype] = 0;
2474 groups->grpnr[gtype] = NULL;
2478 groups->ngrpnr[gtype] = natoms;
2479 snew(groups->grpnr[gtype], natoms);
2480 for (j = 0; (j < natoms); j++)
2482 groups->grpnr[gtype][j] = cbuf[j];
2488 return (bRest && grptp == egrptpPART);
2491 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2494 gmx_groups_t *groups;
2496 int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
2498 int *nrdf2, *na_vcm, na_tot;
2499 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2500 gmx_mtop_atomloop_all_t aloop;
2502 int mb, mol, ftype, as;
2503 gmx_molblock_t *molb;
2504 gmx_moltype_t *molt;
2507 * First calc 3xnr-atoms for each group
2508 * then subtract half a degree of freedom for each constraint
2510 * Only atoms and nuclei contribute to the degrees of freedom...
2515 groups = &mtop->groups;
2516 natoms = mtop->natoms;
2518 /* Allocate one more for a possible rest group */
2519 /* We need to sum degrees of freedom into doubles,
2520 * since floats give too low nrdf's above 3 million atoms.
2522 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2523 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2524 snew(na_vcm, groups->grps[egcVCM].nr+1);
2526 for (i = 0; i < groups->grps[egcTC].nr; i++)
2530 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2535 snew(nrdf2, natoms);
2536 aloop = gmx_mtop_atomloop_all_init(mtop);
2537 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2540 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2542 g = ggrpnr(groups, egcFREEZE, i);
2543 /* Double count nrdf for particle i */
2544 for (d = 0; d < DIM; d++)
2546 if (opts->nFreeze[g][d] == 0)
2551 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2552 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2557 for (mb = 0; mb < mtop->nmolblock; mb++)
2559 molb = &mtop->molblock[mb];
2560 molt = &mtop->moltype[molb->type];
2561 atom = molt->atoms.atom;
2562 for (mol = 0; mol < molb->nmol; mol++)
2564 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2566 ia = molt->ilist[ftype].iatoms;
2567 for (i = 0; i < molt->ilist[ftype].nr; )
2569 /* Subtract degrees of freedom for the constraints,
2570 * if the particles still have degrees of freedom left.
2571 * If one of the particles is a vsite or a shell, then all
2572 * constraint motion will go there, but since they do not
2573 * contribute to the constraints the degrees of freedom do not
2578 if (((atom[ia[1]].ptype == eptNucleus) ||
2579 (atom[ia[1]].ptype == eptAtom)) &&
2580 ((atom[ia[2]].ptype == eptNucleus) ||
2581 (atom[ia[2]].ptype == eptAtom)))
2599 imin = min(imin, nrdf2[ai]);
2600 jmin = min(jmin, nrdf2[aj]);
2603 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2604 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2605 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2606 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2608 ia += interaction_function[ftype].nratoms+1;
2609 i += interaction_function[ftype].nratoms+1;
2612 ia = molt->ilist[F_SETTLE].iatoms;
2613 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2615 /* Subtract 1 dof from every atom in the SETTLE */
2616 for (j = 0; j < 3; j++)
2619 imin = min(2, nrdf2[ai]);
2621 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2622 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2627 as += molt->atoms.nr;
2631 if (ir->ePull == epullCONSTRAINT)
2633 /* Correct nrdf for the COM constraints.
2634 * We correct using the TC and VCM group of the first atom
2635 * in the reference and pull group. If atoms in one pull group
2636 * belong to different TC or VCM groups it is anyhow difficult
2637 * to determine the optimal nrdf assignment.
2640 if (pull->eGeom == epullgPOS)
2643 for (i = 0; i < DIM; i++)
2655 for (i = 0; i < pull->ngrp; i++)
2658 if (pull->grp[0].nat > 0)
2660 /* Subtract 1/2 dof from the reference group */
2661 ai = pull->grp[0].ind[0];
2662 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2664 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
2665 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2669 /* Subtract 1/2 dof from the pulled group */
2670 ai = pull->grp[1+i].ind[0];
2671 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2672 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2673 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2675 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)]]);
2680 if (ir->nstcomm != 0)
2682 /* Subtract 3 from the number of degrees of freedom in each vcm group
2683 * when com translation is removed and 6 when rotation is removed
2686 switch (ir->comm_mode)
2689 n_sub = ndof_com(ir);
2696 gmx_incons("Checking comm_mode");
2699 for (i = 0; i < groups->grps[egcTC].nr; i++)
2701 /* Count the number of atoms of TC group i for every VCM group */
2702 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2707 for (ai = 0; ai < natoms; ai++)
2709 if (ggrpnr(groups, egcTC, ai) == i)
2711 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2715 /* Correct for VCM removal according to the fraction of each VCM
2716 * group present in this TC group.
2718 nrdf_uc = nrdf_tc[i];
2721 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2725 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2727 if (nrdf_vcm[j] > n_sub)
2729 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2730 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2734 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2735 j, nrdf_vcm[j], nrdf_tc[i]);
2740 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2742 opts->nrdf[i] = nrdf_tc[i];
2743 if (opts->nrdf[i] < 0)
2748 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2749 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2758 static void decode_cos(char *s, t_cosines *cosine, gmx_bool bTime)
2761 char format[STRLEN], f1[STRLEN];
2773 sscanf(t, "%d", &(cosine->n));
2780 snew(cosine->a, cosine->n);
2781 snew(cosine->phi, cosine->n);
2783 sprintf(format, "%%*d");
2784 for (i = 0; (i < cosine->n); i++)
2787 strcat(f1, "%lf%lf");
2788 if (sscanf(t, f1, &a, &phi) < 2)
2790 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2793 cosine->phi[i] = phi;
2794 strcat(format, "%*lf%*lf");
2801 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2802 const char *option, const char *val, int flag)
2804 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2805 * But since this is much larger than STRLEN, such a line can not be parsed.
2806 * The real maximum is the number of names that fit in a string: STRLEN/2.
2808 #define EGP_MAX (STRLEN/2)
2809 int nelem, i, j, k, nr;
2810 char *names[EGP_MAX];
2814 gnames = groups->grpname;
2816 nelem = str_nelem(val, EGP_MAX, names);
2819 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2821 nr = groups->grps[egcENER].nr;
2823 for (i = 0; i < nelem/2; i++)
2827 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2833 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2834 names[2*i], option);
2838 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2844 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2845 names[2*i+1], option);
2847 if ((j < nr) && (k < nr))
2849 ir->opts.egp_flags[nr*j+k] |= flag;
2850 ir->opts.egp_flags[nr*k+j] |= flag;
2858 void do_index(const char* mdparin, const char *ndx,
2861 t_inputrec *ir, rvec *v,
2865 gmx_groups_t *groups;
2869 char warnbuf[STRLEN], **gnames;
2870 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2873 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2874 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2875 int i, j, k, restnm;
2877 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2878 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2879 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2880 char warn_buf[STRLEN];
2884 fprintf(stderr, "processing index file...\n");
2890 snew(grps->index, 1);
2892 atoms_all = gmx_mtop_global_atoms(mtop);
2893 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2894 free_t_atoms(&atoms_all, FALSE);
2898 grps = init_index(ndx, &gnames);
2901 groups = &mtop->groups;
2902 natoms = mtop->natoms;
2903 symtab = &mtop->symtab;
2905 snew(groups->grpname, grps->nr+1);
2907 for (i = 0; (i < grps->nr); i++)
2909 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2911 groups->grpname[i] = put_symtab(symtab, "rest");
2913 srenew(gnames, grps->nr+1);
2914 gnames[restnm] = *(groups->grpname[i]);
2915 groups->ngrpname = grps->nr+1;
2917 set_warning_line(wi, mdparin, -1);
2919 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2920 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2921 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2922 if ((ntau_t != ntcg) || (nref_t != ntcg))
2924 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2925 "%d tau-t values", ntcg, nref_t, ntau_t);
2928 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2929 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2930 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2931 nr = groups->grps[egcTC].nr;
2933 snew(ir->opts.nrdf, nr);
2934 snew(ir->opts.tau_t, nr);
2935 snew(ir->opts.ref_t, nr);
2936 if (ir->eI == eiBD && ir->bd_fric == 0)
2938 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2945 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2949 for (i = 0; (i < nr); i++)
2951 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2952 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2954 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2955 warning_error(wi, warn_buf);
2958 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2960 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.");
2963 if (ir->opts.tau_t[i] >= 0)
2965 tau_min = min(tau_min, ir->opts.tau_t[i]);
2968 if (ir->etc != etcNO && ir->nsttcouple == -1)
2970 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2975 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2977 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");
2979 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2981 if (ir->nstpcouple != ir->nsttcouple)
2983 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2984 ir->nstpcouple = ir->nsttcouple = mincouple;
2985 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);
2986 warning_note(wi, warn_buf);
2990 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2991 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2993 if (ETC_ANDERSEN(ir->etc))
2995 if (ir->nsttcouple != 1)
2998 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");
2999 warning_note(wi, warn_buf);
3002 nstcmin = tcouple_min_integration_steps(ir->etc);
3005 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
3007 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
3008 ETCOUPLTYPE(ir->etc),
3010 ir->nsttcouple*ir->delta_t);
3011 warning(wi, warn_buf);
3014 for (i = 0; (i < nr); i++)
3016 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3017 if (ir->opts.ref_t[i] < 0)
3019 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3022 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3023 if we are in this conditional) if mc_temp is negative */
3024 if (ir->expandedvals->mc_temp < 0)
3026 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3030 /* Simulated annealing for each group. There are nr groups */
3031 nSA = str_nelem(anneal, MAXPTR, ptr1);
3032 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3036 if (nSA > 0 && nSA != nr)
3038 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3042 snew(ir->opts.annealing, nr);
3043 snew(ir->opts.anneal_npoints, nr);
3044 snew(ir->opts.anneal_time, nr);
3045 snew(ir->opts.anneal_temp, nr);
3046 for (i = 0; i < nr; i++)
3048 ir->opts.annealing[i] = eannNO;
3049 ir->opts.anneal_npoints[i] = 0;
3050 ir->opts.anneal_time[i] = NULL;
3051 ir->opts.anneal_temp[i] = NULL;
3056 for (i = 0; i < nr; i++)
3058 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3060 ir->opts.annealing[i] = eannNO;
3062 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3064 ir->opts.annealing[i] = eannSINGLE;
3067 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3069 ir->opts.annealing[i] = eannPERIODIC;
3075 /* Read the other fields too */
3076 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3077 if (nSA_points != nSA)
3079 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3081 for (k = 0, i = 0; i < nr; i++)
3083 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3084 if (ir->opts.anneal_npoints[i] == 1)
3086 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3088 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3089 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3090 k += ir->opts.anneal_npoints[i];
3093 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3096 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3098 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3101 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3104 for (i = 0, k = 0; i < nr; i++)
3107 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3109 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3110 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3113 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3115 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3121 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3123 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3124 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3127 if (ir->opts.anneal_temp[i][j] < 0)
3129 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3134 /* Print out some summary information, to make sure we got it right */
3135 for (i = 0, k = 0; i < nr; i++)
3137 if (ir->opts.annealing[i] != eannNO)
3139 j = groups->grps[egcTC].nm_ind[i];
3140 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3141 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3142 ir->opts.anneal_npoints[i]);
3143 fprintf(stderr, "Time (ps) Temperature (K)\n");
3144 /* All terms except the last one */
3145 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3147 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3150 /* Finally the last one */
3151 j = ir->opts.anneal_npoints[i]-1;
3152 if (ir->opts.annealing[i] == eannSINGLE)
3154 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3158 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3159 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3161 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3170 if (ir->ePull != epullNO)
3172 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3177 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3180 nacc = str_nelem(acc, MAXPTR, ptr1);
3181 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3182 if (nacg*DIM != nacc)
3184 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3187 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3188 restnm, egrptpALL_GENREST, bVerbose, wi);
3189 nr = groups->grps[egcACC].nr;
3190 snew(ir->opts.acc, nr);
3191 ir->opts.ngacc = nr;
3193 for (i = k = 0; (i < nacg); i++)
3195 for (j = 0; (j < DIM); j++, k++)
3197 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3200 for (; (i < nr); i++)
3202 for (j = 0; (j < DIM); j++)
3204 ir->opts.acc[i][j] = 0;
3208 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3209 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3210 if (nfrdim != DIM*nfreeze)
3212 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3215 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3216 restnm, egrptpALL_GENREST, bVerbose, wi);
3217 nr = groups->grps[egcFREEZE].nr;
3218 ir->opts.ngfrz = nr;
3219 snew(ir->opts.nFreeze, nr);
3220 for (i = k = 0; (i < nfreeze); i++)
3222 for (j = 0; (j < DIM); j++, k++)
3224 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3225 if (!ir->opts.nFreeze[i][j])
3227 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3229 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3230 "(not %s)", ptr1[k]);
3231 warning(wi, warn_buf);
3236 for (; (i < nr); i++)
3238 for (j = 0; (j < DIM); j++)
3240 ir->opts.nFreeze[i][j] = 0;
3244 nenergy = str_nelem(energy, MAXPTR, ptr1);
3245 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3246 restnm, egrptpALL_GENREST, bVerbose, wi);
3247 add_wall_energrps(groups, ir->nwall, symtab);
3248 ir->opts.ngener = groups->grps[egcENER].nr;
3249 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3251 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3252 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3255 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3256 "This may lead to artifacts.\n"
3257 "In most cases one should use one group for the whole system.");
3260 /* Now we have filled the freeze struct, so we can calculate NRDF */
3261 calc_nrdf(mtop, ir, gnames);
3267 /* Must check per group! */
3268 for (i = 0; (i < ir->opts.ngtc); i++)
3270 ntot += ir->opts.nrdf[i];
3272 if (ntot != (DIM*natoms))
3274 fac = sqrt(ntot/(DIM*natoms));
3277 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3278 "and removal of center of mass motion\n", fac);
3280 for (i = 0; (i < natoms); i++)
3282 svmul(fac, v[i], v[i]);
3287 nuser = str_nelem(user1, MAXPTR, ptr1);
3288 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3289 restnm, egrptpALL_GENREST, bVerbose, wi);
3290 nuser = str_nelem(user2, MAXPTR, ptr1);
3291 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3292 restnm, egrptpALL_GENREST, bVerbose, wi);
3293 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3294 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3295 restnm, egrptpONE, bVerbose, wi);
3296 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3297 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3298 restnm, egrptpALL_GENREST, bVerbose, wi);
3300 /* QMMM input processing */
3301 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3302 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3303 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3304 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3306 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3307 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3309 /* group rest, if any, is always MM! */
3310 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3311 restnm, egrptpALL_GENREST, bVerbose, wi);
3312 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3313 ir->opts.ngQM = nQMg;
3314 snew(ir->opts.QMmethod, nr);
3315 snew(ir->opts.QMbasis, nr);
3316 for (i = 0; i < nr; i++)
3318 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3319 * converted to the corresponding enum in names.c
3321 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3323 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3327 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3328 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3329 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3330 snew(ir->opts.QMmult, nr);
3331 snew(ir->opts.QMcharge, nr);
3332 snew(ir->opts.bSH, nr);
3334 for (i = 0; i < nr; i++)
3336 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3337 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3338 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3341 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3342 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3343 snew(ir->opts.CASelectrons, nr);
3344 snew(ir->opts.CASorbitals, nr);
3345 for (i = 0; i < nr; i++)
3347 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3348 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3350 /* special optimization options */
3352 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3353 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3354 snew(ir->opts.bOPT, nr);
3355 snew(ir->opts.bTS, nr);
3356 for (i = 0; i < nr; i++)
3358 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3359 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3361 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3362 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3363 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3364 snew(ir->opts.SAon, nr);
3365 snew(ir->opts.SAoff, nr);
3366 snew(ir->opts.SAsteps, nr);
3368 for (i = 0; i < nr; i++)
3370 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3371 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3372 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3374 /* end of QMMM input */
3378 for (i = 0; (i < egcNR); i++)
3380 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3381 for (j = 0; (j < groups->grps[i].nr); j++)
3383 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3385 fprintf(stderr, "\n");
3389 nr = groups->grps[egcENER].nr;
3390 snew(ir->opts.egp_flags, nr*nr);
3392 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3393 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3395 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3397 if (bExcl && EEL_FULL(ir->coulombtype))
3399 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3402 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3403 if (bTable && !(ir->vdwtype == evdwUSER) &&
3404 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3405 !(ir->coulombtype == eelPMEUSERSWITCH))
3407 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3410 decode_cos(efield_x, &(ir->ex[XX]), FALSE);
3411 decode_cos(efield_xt, &(ir->et[XX]), TRUE);
3412 decode_cos(efield_y, &(ir->ex[YY]), FALSE);
3413 decode_cos(efield_yt, &(ir->et[YY]), TRUE);
3414 decode_cos(efield_z, &(ir->ex[ZZ]), FALSE);
3415 decode_cos(efield_zt, &(ir->et[ZZ]), TRUE);
3419 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3422 for (i = 0; (i < grps->nr); i++)
3434 static void check_disre(gmx_mtop_t *mtop)
3436 gmx_ffparams_t *ffparams;
3437 t_functype *functype;
3439 int i, ndouble, ftype;
3440 int label, old_label;
3442 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3444 ffparams = &mtop->ffparams;
3445 functype = ffparams->functype;
3446 ip = ffparams->iparams;
3449 for (i = 0; i < ffparams->ntypes; i++)
3451 ftype = functype[i];
3452 if (ftype == F_DISRES)
3454 label = ip[i].disres.label;
3455 if (label == old_label)
3457 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3465 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3466 "probably the parameters for multiple pairs in one restraint "
3467 "are not identical\n", ndouble);
3472 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3473 gmx_bool posres_only,
3477 gmx_mtop_ilistloop_t iloop;
3487 for (d = 0; d < DIM; d++)
3489 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3491 /* Check for freeze groups */
3492 for (g = 0; g < ir->opts.ngfrz; g++)
3494 for (d = 0; d < DIM; d++)
3496 if (ir->opts.nFreeze[g][d] != 0)
3504 /* Check for position restraints */
3505 iloop = gmx_mtop_ilistloop_init(sys);
3506 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3509 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3511 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3513 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3514 for (d = 0; d < DIM; d++)
3516 if (pr->posres.fcA[d] != 0)
3525 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3528 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3532 int i, m, g, nmol, npct;
3533 gmx_bool bCharge, bAcc;
3534 real gdt_max, *mgrp, mt;
3536 gmx_mtop_atomloop_block_t aloopb;
3537 gmx_mtop_atomloop_all_t aloop;
3540 char warn_buf[STRLEN];
3542 set_warning_line(wi, mdparin, -1);
3544 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3545 ir->comm_mode == ecmNO &&
3546 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
3547 !ETC_ANDERSEN(ir->etc))
3549 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");
3552 /* Check for pressure coupling with absolute position restraints */
3553 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3555 absolute_reference(ir, sys, TRUE, AbsRef);
3557 for (m = 0; m < DIM; m++)
3559 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3561 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3569 aloopb = gmx_mtop_atomloop_block_init(sys);
3570 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3572 if (atom->q != 0 || atom->qB != 0)
3580 if (EEL_FULL(ir->coulombtype))
3583 "You are using full electrostatics treatment %s for a system without charges.\n"
3584 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3585 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3586 warning(wi, err_buf);
3591 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3594 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3595 "You might want to consider using %s electrostatics.\n",
3597 warning_note(wi, err_buf);
3601 /* Generalized reaction field */
3602 if (ir->opts.ngtc == 0)
3604 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3606 CHECK(ir->coulombtype == eelGRF);
3610 sprintf(err_buf, "When using coulombtype = %s"
3611 " ref-t for temperature coupling should be > 0",
3613 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3616 if (ir->eI == eiSD1 &&
3617 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3618 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3620 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3621 warning_note(wi, warn_buf);
3625 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3627 for (m = 0; (m < DIM); m++)
3629 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3638 snew(mgrp, sys->groups.grps[egcACC].nr);
3639 aloop = gmx_mtop_atomloop_all_init(sys);
3640 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3642 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3645 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3647 for (m = 0; (m < DIM); m++)
3649 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3653 for (m = 0; (m < DIM); m++)
3655 if (fabs(acc[m]) > 1e-6)
3657 const char *dim[DIM] = { "X", "Y", "Z" };
3659 "Net Acceleration in %s direction, will %s be corrected\n",
3660 dim[m], ir->nstcomm != 0 ? "" : "not");
3661 if (ir->nstcomm != 0 && m < ndof_com(ir))
3664 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3666 ir->opts.acc[i][m] -= acc[m];
3674 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3675 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3677 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3680 if (ir->ePull != epullNO)
3682 if (ir->pull->grp[0].nat == 0)
3684 absolute_reference(ir, sys, FALSE, AbsRef);
3685 for (m = 0; m < DIM; m++)
3687 if (ir->pull->dim[m] && !AbsRef[m])
3689 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.");
3695 if (ir->pull->eGeom == epullgDIRPBC)
3697 for (i = 0; i < 3; i++)
3699 for (m = 0; m <= i; m++)
3701 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3702 ir->deform[i][m] != 0)
3704 for (g = 1; g < ir->pull->ngrp; g++)
3706 if (ir->pull->grp[g].vec[m] != 0)
3708 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3720 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3724 char warn_buf[STRLEN];
3727 ptr = check_box(ir->ePBC, box);
3730 warning_error(wi, ptr);
3733 if (bConstr && ir->eConstrAlg == econtSHAKE)
3735 if (ir->shake_tol <= 0.0)
3737 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3739 warning_error(wi, warn_buf);
3742 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3744 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3745 if (ir->epc == epcNO)
3747 warning(wi, warn_buf);
3751 warning_error(wi, warn_buf);
3756 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3758 /* If we have Lincs constraints: */
3759 if (ir->eI == eiMD && ir->etc == etcNO &&
3760 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3762 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3763 warning_note(wi, warn_buf);
3766 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3768 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3769 warning_note(wi, warn_buf);
3771 if (ir->epc == epcMTTK)
3773 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3777 if (bConstr && ir->epc == epcMTTK)
3779 warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
3782 if (ir->LincsWarnAngle > 90.0)
3784 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3785 warning(wi, warn_buf);
3786 ir->LincsWarnAngle = 90.0;
3789 if (ir->ePBC != epbcNONE)
3791 if (ir->nstlist == 0)
3793 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3795 bTWIN = (ir->rlistlong > ir->rlist);
3796 if (ir->ns_type == ensGRID)
3798 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3800 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",
3801 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3802 warning_error(wi, warn_buf);
3807 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3808 if (2*ir->rlistlong >= min_size)
3810 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3811 warning_error(wi, warn_buf);
3814 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3821 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3825 real rvdw1, rvdw2, rcoul1, rcoul2;
3826 char warn_buf[STRLEN];
3828 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
3832 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3837 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
3843 if (rvdw1 + rvdw2 > ir->rlist ||
3844 rcoul1 + rcoul2 > ir->rlist)
3847 "The sum of the two largest charge group radii (%f) "
3848 "is larger than rlist (%f)\n",
3849 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
3850 warning(wi, warn_buf);
3854 /* Here we do not use the zero at cut-off macro,
3855 * since user defined interactions might purposely
3856 * not be zero at the cut-off.
3858 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
3859 ir->vdw_modifier != eintmodNONE) &&
3860 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
3862 sprintf(warn_buf, "The sum of the two largest charge group "
3863 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
3864 "With exact cut-offs, better performance can be "
3865 "obtained with cutoff-scheme = %s, because it "
3866 "does not use charge groups at all.",
3868 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3869 ir->rlistlong, ir->rvdw,
3870 ecutscheme_names[ecutsVERLET]);
3873 warning(wi, warn_buf);
3877 warning_note(wi, warn_buf);
3880 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
3881 ir->coulomb_modifier != eintmodNONE) &&
3882 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3884 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
3885 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
3887 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3888 ir->rlistlong, ir->rcoulomb,
3889 ecutscheme_names[ecutsVERLET]);
3892 warning(wi, warn_buf);
3896 warning_note(wi, warn_buf);