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)
1093 if (ir->rcoulomb_switch/ir->rcoulomb < 0.9499)
1095 sprintf(warn_buf, "The switching range for %s should be 5%% or less, energy conservation will be good anyhow, since ewald_rtol = %g",
1096 eel_names[ir->coulombtype],
1098 warning(wi, warn_buf);
1102 if (EEL_FULL(ir->coulombtype))
1104 if (ir->coulombtype == eelPMESWITCH || ir->coulombtype == eelPMEUSER ||
1105 ir->coulombtype == eelPMEUSERSWITCH)
1107 sprintf(err_buf, "With coulombtype = %s, rcoulomb must be <= rlist",
1108 eel_names[ir->coulombtype]);
1109 CHECK(ir->rcoulomb > ir->rlist);
1111 else if (ir->cutoff_scheme == ecutsGROUP && ir->coulomb_modifier == eintmodNONE)
1113 if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD)
1116 "With coulombtype = %s (without modifier), rcoulomb must be equal to rlist,\n"
1117 "or rlistlong if nstcalclr=1. For optimal energy conservation,consider using\n"
1118 "a potential modifier.", eel_names[ir->coulombtype]);
1119 if (ir->nstcalclr == 1)
1121 CHECK(ir->rcoulomb != ir->rlist && ir->rcoulomb != ir->rlistlong);
1125 CHECK(ir->rcoulomb != ir->rlist);
1131 if (EEL_PME(ir->coulombtype))
1133 if (ir->pme_order < 3)
1135 warning_error(wi, "pme-order can not be smaller than 3");
1139 if (ir->nwall == 2 && EEL_FULL(ir->coulombtype))
1141 if (ir->ewald_geometry == eewg3D)
1143 sprintf(warn_buf, "With pbc=%s you should use ewald-geometry=%s",
1144 epbc_names[ir->ePBC], eewg_names[eewg3DC]);
1145 warning(wi, warn_buf);
1147 /* This check avoids extra pbc coding for exclusion corrections */
1148 sprintf(err_buf, "wall-ewald-zfac should be >= 2");
1149 CHECK(ir->wall_ewald_zfac < 2);
1152 if (EVDW_SWITCHED(ir->vdwtype))
1154 sprintf(err_buf, "With vdwtype = %s rvdw-switch must be < rvdw. Or, better - use a potential modifier.",
1155 evdw_names[ir->vdwtype]);
1156 CHECK(ir->rvdw_switch >= ir->rvdw);
1158 else if (ir->vdwtype == evdwCUT)
1160 if (ir->cutoff_scheme == ecutsGROUP && ir->vdw_modifier == eintmodNONE)
1162 sprintf(err_buf, "With vdwtype = %s, rvdw must be >= rlist unless you use a potential modifier", evdw_names[ir->vdwtype]);
1163 CHECK(ir->rlist > ir->rvdw);
1166 if (ir->cutoff_scheme == ecutsGROUP)
1168 if (((ir->coulomb_modifier != eintmodNONE && ir->rcoulomb == ir->rlist) ||
1169 (ir->vdw_modifier != eintmodNONE && ir->rvdw == ir->rlist)) &&
1172 warning_note(wi, "With exact cut-offs, rlist should be "
1173 "larger than rcoulomb and rvdw, so that there "
1174 "is a buffer region for particle motion "
1175 "between neighborsearch steps");
1178 if (EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype)
1179 && (ir->rlistlong <= ir->rcoulomb))
1181 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rcoulomb.",
1182 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1183 warning_note(wi, warn_buf);
1185 if (EVDW_SWITCHED(ir->vdwtype) && (ir->rlistlong <= ir->rvdw))
1187 sprintf(warn_buf, "For energy conservation with switch/shift potentials, %s should be 0.1 to 0.3 nm larger than rvdw.",
1188 IR_TWINRANGE(*ir) ? "rlistlong" : "rlist");
1189 warning_note(wi, warn_buf);
1193 if (ir->vdwtype == evdwUSER && ir->eDispCorr != edispcNO)
1195 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.");
1198 if (ir->nstlist == -1)
1200 sprintf(err_buf, "With nstlist=-1 rvdw and rcoulomb should be smaller than rlist to account for diffusion and possibly charge-group radii");
1201 CHECK(ir->rvdw >= ir->rlist || ir->rcoulomb >= ir->rlist);
1203 sprintf(err_buf, "nstlist can not be smaller than -1");
1204 CHECK(ir->nstlist < -1);
1206 if (ir->eI == eiLBFGS && (ir->coulombtype == eelCUT || ir->vdwtype == evdwCUT)
1209 warning(wi, "For efficient BFGS minimization, use switch/shift/pme instead of cut-off.");
1212 if (ir->eI == eiLBFGS && ir->nbfgscorr <= 0)
1214 warning(wi, "Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
1217 /* ENERGY CONSERVATION */
1218 if (ir_NVE(ir) && ir->cutoff_scheme == ecutsGROUP)
1220 if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0 && ir->vdw_modifier == eintmodNONE)
1222 sprintf(warn_buf, "You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
1223 evdw_names[evdwSHIFT]);
1224 warning_note(wi, warn_buf);
1226 if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0 && ir->coulomb_modifier == eintmodNONE)
1228 sprintf(warn_buf, "You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
1229 eel_names[eelPMESWITCH], eel_names[eelRF_ZERO]);
1230 warning_note(wi, warn_buf);
1234 /* IMPLICIT SOLVENT */
1235 if (ir->coulombtype == eelGB_NOTUSED)
1237 ir->coulombtype = eelCUT;
1238 ir->implicit_solvent = eisGBSA;
1239 fprintf(stderr, "Note: Old option for generalized born electrostatics given:\n"
1240 "Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
1241 "setting implicit-solvent value to \"GBSA\" in input section.\n");
1244 if (ir->sa_algorithm == esaSTILL)
1246 sprintf(err_buf, "Still SA algorithm not available yet, use %s or %s instead\n", esa_names[esaAPPROX], esa_names[esaNO]);
1247 CHECK(ir->sa_algorithm == esaSTILL);
1250 if (ir->implicit_solvent == eisGBSA)
1252 sprintf(err_buf, "With GBSA implicit solvent, rgbradii must be equal to rlist.");
1253 CHECK(ir->rgbradii != ir->rlist);
1255 if (ir->coulombtype != eelCUT)
1257 sprintf(err_buf, "With GBSA, coulombtype must be equal to %s\n", eel_names[eelCUT]);
1258 CHECK(ir->coulombtype != eelCUT);
1260 if (ir->vdwtype != evdwCUT)
1262 sprintf(err_buf, "With GBSA, vdw-type must be equal to %s\n", evdw_names[evdwCUT]);
1263 CHECK(ir->vdwtype != evdwCUT);
1265 if (ir->nstgbradii < 1)
1267 sprintf(warn_buf, "Using GBSA with nstgbradii<1, setting nstgbradii=1");
1268 warning_note(wi, warn_buf);
1271 if (ir->sa_algorithm == esaNO)
1273 sprintf(warn_buf, "No SA (non-polar) calculation requested together with GB. Are you sure this is what you want?\n");
1274 warning_note(wi, warn_buf);
1276 if (ir->sa_surface_tension < 0 && ir->sa_algorithm != esaNO)
1278 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");
1279 warning_note(wi, warn_buf);
1281 if (ir->gb_algorithm == egbSTILL)
1283 ir->sa_surface_tension = 0.0049 * CAL2JOULE * 100;
1287 ir->sa_surface_tension = 0.0054 * CAL2JOULE * 100;
1290 if (ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO)
1292 sprintf(err_buf, "Surface tension set to 0 while SA-calculation requested\n");
1293 CHECK(ir->sa_surface_tension == 0 && ir->sa_algorithm != esaNO);
1300 if (ir->cutoff_scheme != ecutsGROUP)
1302 warning_error(wi, "AdresS simulation supports only cutoff-scheme=group");
1306 warning_error(wi, "AdresS simulation supports only stochastic dynamics");
1308 if (ir->epc != epcNO)
1310 warning_error(wi, "AdresS simulation does not support pressure coupling");
1312 if (EEL_FULL(ir->coulombtype))
1314 warning_error(wi, "AdresS simulation does not support long-range electrostatics");
1319 /* count the number of text elemets separated by whitespace in a string.
1320 str = the input string
1321 maxptr = the maximum number of allowed elements
1322 ptr = the output array of pointers to the first character of each element
1323 returns: the number of elements. */
1324 int str_nelem(const char *str, int maxptr, char *ptr[])
1329 copy0 = strdup(str);
1332 while (*copy != '\0')
1336 gmx_fatal(FARGS, "Too many groups on line: '%s' (max is %d)",
1344 while ((*copy != '\0') && !isspace(*copy))
1363 /* interpret a number of doubles from a string and put them in an array,
1364 after allocating space for them.
1365 str = the input string
1366 n = the (pre-allocated) number of doubles read
1367 r = the output array of doubles. */
1368 static void parse_n_real(char *str, int *n, real **r)
1373 *n = str_nelem(str, MAXPTR, ptr);
1376 for (i = 0; i < *n; i++)
1378 (*r)[i] = strtod(ptr[i], NULL);
1382 static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN], char weights[STRLEN])
1385 int i, j, max_n_lambda, nweights, nfep[efptNR];
1386 t_lambda *fep = ir->fepvals;
1387 t_expanded *expand = ir->expandedvals;
1388 real **count_fep_lambdas;
1389 gmx_bool bOneLambda = TRUE;
1391 snew(count_fep_lambdas, efptNR);
1393 /* FEP input processing */
1394 /* first, identify the number of lambda values for each type.
1395 All that are nonzero must have the same number */
1397 for (i = 0; i < efptNR; i++)
1399 parse_n_real(fep_lambda[i], &(nfep[i]), &(count_fep_lambdas[i]));
1402 /* now, determine the number of components. All must be either zero, or equal. */
1405 for (i = 0; i < efptNR; i++)
1407 if (nfep[i] > max_n_lambda)
1409 max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
1410 must have the same number if its not zero.*/
1415 for (i = 0; i < efptNR; i++)
1419 ir->fepvals->separate_dvdl[i] = FALSE;
1421 else if (nfep[i] == max_n_lambda)
1423 if (i != efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
1424 respect to the temperature currently */
1426 ir->fepvals->separate_dvdl[i] = TRUE;
1431 gmx_fatal(FARGS, "Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
1432 nfep[i], efpt_names[i], max_n_lambda);
1435 /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
1436 ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
1438 /* the number of lambdas is the number we've read in, which is either zero
1439 or the same for all */
1440 fep->n_lambda = max_n_lambda;
1442 /* allocate space for the array of lambda values */
1443 snew(fep->all_lambda, efptNR);
1444 /* if init_lambda is defined, we need to set lambda */
1445 if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
1447 ir->fepvals->separate_dvdl[efptFEP] = TRUE;
1449 /* otherwise allocate the space for all of the lambdas, and transfer the data */
1450 for (i = 0; i < efptNR; i++)
1452 snew(fep->all_lambda[i], fep->n_lambda);
1453 if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
1456 for (j = 0; j < fep->n_lambda; j++)
1458 fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
1460 sfree(count_fep_lambdas[i]);
1463 sfree(count_fep_lambdas);
1465 /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
1466 bookkeeping -- for now, init_lambda */
1468 if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0))
1470 for (i = 0; i < fep->n_lambda; i++)
1472 fep->all_lambda[efptFEP][i] = fep->init_lambda;
1476 /* check to see if only a single component lambda is defined, and soft core is defined.
1477 In this case, turn on coulomb soft core */
1479 if (max_n_lambda == 0)
1485 for (i = 0; i < efptNR; i++)
1487 if ((nfep[i] != 0) && (i != efptFEP))
1493 if ((bOneLambda) && (fep->sc_alpha > 0))
1495 fep->bScCoul = TRUE;
1498 /* Fill in the others with the efptFEP if they are not explicitly
1499 specified (i.e. nfep[i] == 0). This means if fep is not defined,
1500 they are all zero. */
1502 for (i = 0; i < efptNR; i++)
1504 if ((nfep[i] == 0) && (i != efptFEP))
1506 for (j = 0; j < fep->n_lambda; j++)
1508 fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
1514 /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
1515 if (fep->sc_r_power == 48)
1517 if (fep->sc_alpha > 0.1)
1519 gmx_fatal(FARGS, "sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
1523 expand = ir->expandedvals;
1524 /* now read in the weights */
1525 parse_n_real(weights, &nweights, &(expand->init_lambda_weights));
1528 snew(expand->init_lambda_weights, fep->n_lambda); /* initialize to zero */
1530 else if (nweights != fep->n_lambda)
1532 gmx_fatal(FARGS, "Number of weights (%d) is not equal to number of lambda values (%d)",
1533 nweights, fep->n_lambda);
1535 if ((expand->nstexpanded < 0) && (ir->efep != efepNO))
1537 expand->nstexpanded = fep->nstdhdl;
1538 /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
1540 if ((expand->nstexpanded < 0) && ir->bSimTemp)
1542 expand->nstexpanded = 2*(int)(ir->opts.tau_t[0]/ir->delta_t);
1543 /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to
1544 2*tau_t just to be careful so it's not to frequent */
1549 static void do_simtemp_params(t_inputrec *ir)
1552 snew(ir->simtempvals->temperatures, ir->fepvals->n_lambda);
1553 GetSimTemps(ir->fepvals->n_lambda, ir->simtempvals, ir->fepvals->all_lambda[efptTEMPERATURE]);
1558 static void do_wall_params(t_inputrec *ir,
1559 char *wall_atomtype, char *wall_density,
1563 char *names[MAXPTR];
1566 opts->wall_atomtype[0] = NULL;
1567 opts->wall_atomtype[1] = NULL;
1569 ir->wall_atomtype[0] = -1;
1570 ir->wall_atomtype[1] = -1;
1571 ir->wall_density[0] = 0;
1572 ir->wall_density[1] = 0;
1576 nstr = str_nelem(wall_atomtype, MAXPTR, names);
1577 if (nstr != ir->nwall)
1579 gmx_fatal(FARGS, "Expected %d elements for wall_atomtype, found %d",
1582 for (i = 0; i < ir->nwall; i++)
1584 opts->wall_atomtype[i] = strdup(names[i]);
1587 if (ir->wall_type == ewt93 || ir->wall_type == ewt104)
1589 nstr = str_nelem(wall_density, MAXPTR, names);
1590 if (nstr != ir->nwall)
1592 gmx_fatal(FARGS, "Expected %d elements for wall-density, found %d", ir->nwall, nstr);
1594 for (i = 0; i < ir->nwall; i++)
1596 sscanf(names[i], "%lf", &dbl);
1599 gmx_fatal(FARGS, "wall-density[%d] = %f\n", i, dbl);
1601 ir->wall_density[i] = dbl;
1607 static void add_wall_energrps(gmx_groups_t *groups, int nwall, t_symtab *symtab)
1615 srenew(groups->grpname, groups->ngrpname+nwall);
1616 grps = &(groups->grps[egcENER]);
1617 srenew(grps->nm_ind, grps->nr+nwall);
1618 for (i = 0; i < nwall; i++)
1620 sprintf(str, "wall%d", i);
1621 groups->grpname[groups->ngrpname] = put_symtab(symtab, str);
1622 grps->nm_ind[grps->nr++] = groups->ngrpname++;
1627 void read_expandedparams(int *ninp_p, t_inpfile **inp_p,
1628 t_expanded *expand, warninp_t wi)
1630 int ninp, nerror = 0;
1636 /* read expanded ensemble parameters */
1637 CCTYPE ("expanded ensemble variables");
1638 ITYPE ("nstexpanded", expand->nstexpanded, -1);
1639 EETYPE("lmc-stats", expand->elamstats, elamstats_names);
1640 EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
1641 EETYPE("lmc-weights-equil", expand->elmceq, elmceq_names);
1642 ITYPE ("weight-equil-number-all-lambda", expand->equil_n_at_lam, -1);
1643 ITYPE ("weight-equil-number-samples", expand->equil_samples, -1);
1644 ITYPE ("weight-equil-number-steps", expand->equil_steps, -1);
1645 RTYPE ("weight-equil-wl-delta", expand->equil_wl_delta, -1);
1646 RTYPE ("weight-equil-count-ratio", expand->equil_ratio, -1);
1647 CCTYPE("Seed for Monte Carlo in lambda space");
1648 ITYPE ("lmc-seed", expand->lmc_seed, -1);
1649 RTYPE ("mc-temperature", expand->mc_temp, -1);
1650 ITYPE ("lmc-repeats", expand->lmc_repeats, 1);
1651 ITYPE ("lmc-gibbsdelta", expand->gibbsdeltalam, -1);
1652 ITYPE ("lmc-forced-nstart", expand->lmc_forced_nstart, 0);
1653 EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
1654 ITYPE("nst-transition-matrix", expand->nstTij, -1);
1655 ITYPE ("mininum-var-min", expand->minvarmin, 100); /*default is reasonable */
1656 ITYPE ("weight-c-range", expand->c_range, 0); /* default is just C=0 */
1657 RTYPE ("wl-scale", expand->wl_scale, 0.8);
1658 RTYPE ("wl-ratio", expand->wl_ratio, 0.8);
1659 RTYPE ("init-wl-delta", expand->init_wl_delta, 1.0);
1660 EETYPE("wl-oneovert", expand->bWLoneovert, yesno_names);
1668 void get_ir(const char *mdparin, const char *mdparout,
1669 t_inputrec *ir, t_gromppopts *opts,
1673 double dumdub[2][6];
1677 char warn_buf[STRLEN];
1678 t_lambda *fep = ir->fepvals;
1679 t_expanded *expand = ir->expandedvals;
1681 inp = read_inpfile(mdparin, &ninp, NULL, wi);
1683 snew(dumstr[0], STRLEN);
1684 snew(dumstr[1], STRLEN);
1686 /* remove the following deprecated commands */
1689 REM_TYPE("domain-decomposition");
1690 REM_TYPE("andersen-seed");
1692 REM_TYPE("dihre-fc");
1693 REM_TYPE("dihre-tau");
1694 REM_TYPE("nstdihreout");
1695 REM_TYPE("nstcheckpoint");
1697 /* replace the following commands with the clearer new versions*/
1698 REPL_TYPE("unconstrained-start", "continuation");
1699 REPL_TYPE("foreign-lambda", "fep-lambdas");
1701 CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
1702 CTYPE ("Preprocessor information: use cpp syntax.");
1703 CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
1704 STYPE ("include", opts->include, NULL);
1705 CTYPE ("e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)");
1706 STYPE ("define", opts->define, NULL);
1708 CCTYPE ("RUN CONTROL PARAMETERS");
1709 EETYPE("integrator", ir->eI, ei_names);
1710 CTYPE ("Start time and timestep in ps");
1711 RTYPE ("tinit", ir->init_t, 0.0);
1712 RTYPE ("dt", ir->delta_t, 0.001);
1713 STEPTYPE ("nsteps", ir->nsteps, 0);
1714 CTYPE ("For exact run continuation or redoing part of a run");
1715 STEPTYPE ("init-step", ir->init_step, 0);
1716 CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
1717 ITYPE ("simulation-part", ir->simulation_part, 1);
1718 CTYPE ("mode for center of mass motion removal");
1719 EETYPE("comm-mode", ir->comm_mode, ecm_names);
1720 CTYPE ("number of steps for center of mass motion removal");
1721 ITYPE ("nstcomm", ir->nstcomm, 100);
1722 CTYPE ("group(s) for center of mass motion removal");
1723 STYPE ("comm-grps", vcm, NULL);
1725 CCTYPE ("LANGEVIN DYNAMICS OPTIONS");
1726 CTYPE ("Friction coefficient (amu/ps) and random seed");
1727 RTYPE ("bd-fric", ir->bd_fric, 0.0);
1728 ITYPE ("ld-seed", ir->ld_seed, 1993);
1731 CCTYPE ("ENERGY MINIMIZATION OPTIONS");
1732 CTYPE ("Force tolerance and initial step-size");
1733 RTYPE ("emtol", ir->em_tol, 10.0);
1734 RTYPE ("emstep", ir->em_stepsize, 0.01);
1735 CTYPE ("Max number of iterations in relax-shells");
1736 ITYPE ("niter", ir->niter, 20);
1737 CTYPE ("Step size (ps^2) for minimization of flexible constraints");
1738 RTYPE ("fcstep", ir->fc_stepsize, 0);
1739 CTYPE ("Frequency of steepest descents steps when doing CG");
1740 ITYPE ("nstcgsteep", ir->nstcgsteep, 1000);
1741 ITYPE ("nbfgscorr", ir->nbfgscorr, 10);
1743 CCTYPE ("TEST PARTICLE INSERTION OPTIONS");
1744 RTYPE ("rtpi", ir->rtpi, 0.05);
1746 /* Output options */
1747 CCTYPE ("OUTPUT CONTROL OPTIONS");
1748 CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
1749 ITYPE ("nstxout", ir->nstxout, 0);
1750 ITYPE ("nstvout", ir->nstvout, 0);
1751 ITYPE ("nstfout", ir->nstfout, 0);
1752 ir->nstcheckpoint = 1000;
1753 CTYPE ("Output frequency for energies to log file and energy file");
1754 ITYPE ("nstlog", ir->nstlog, 1000);
1755 ITYPE ("nstcalcenergy", ir->nstcalcenergy, 100);
1756 ITYPE ("nstenergy", ir->nstenergy, 1000);
1757 CTYPE ("Output frequency and precision for .xtc file");
1758 ITYPE ("nstxtcout", ir->nstxtcout, 0);
1759 RTYPE ("xtc-precision", ir->xtcprec, 1000.0);
1760 CTYPE ("This selects the subset of atoms for the .xtc file. You can");
1761 CTYPE ("select multiple groups. By default all atoms will be written.");
1762 STYPE ("xtc-grps", xtc_grps, NULL);
1763 CTYPE ("Selection of energy groups");
1764 STYPE ("energygrps", energy, NULL);
1766 /* Neighbor searching */
1767 CCTYPE ("NEIGHBORSEARCHING PARAMETERS");
1768 CTYPE ("cut-off scheme (group: using charge groups, Verlet: particle based cut-offs)");
1769 EETYPE("cutoff-scheme", ir->cutoff_scheme, ecutscheme_names);
1770 CTYPE ("nblist update frequency");
1771 ITYPE ("nstlist", ir->nstlist, 10);
1772 CTYPE ("ns algorithm (simple or grid)");
1773 EETYPE("ns-type", ir->ns_type, ens_names);
1774 /* set ndelta to the optimal value of 2 */
1776 CTYPE ("Periodic boundary conditions: xyz, no, xy");
1777 EETYPE("pbc", ir->ePBC, epbc_names);
1778 EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
1779 CTYPE ("Allowed energy drift due to the Verlet buffer in kJ/mol/ps per atom,");
1780 CTYPE ("a value of -1 means: use rlist");
1781 RTYPE("verlet-buffer-drift", ir->verletbuf_drift, 0.005);
1782 CTYPE ("nblist cut-off");
1783 RTYPE ("rlist", ir->rlist, 1.0);
1784 CTYPE ("long-range cut-off for switched potentials");
1785 RTYPE ("rlistlong", ir->rlistlong, -1);
1786 ITYPE ("nstcalclr", ir->nstcalclr, -1);
1788 /* Electrostatics */
1789 CCTYPE ("OPTIONS FOR ELECTROSTATICS AND VDW");
1790 CTYPE ("Method for doing electrostatics");
1791 EETYPE("coulombtype", ir->coulombtype, eel_names);
1792 EETYPE("coulomb-modifier", ir->coulomb_modifier, eintmod_names);
1793 CTYPE ("cut-off lengths");
1794 RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
1795 RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
1796 CTYPE ("Relative dielectric constant for the medium and the reaction field");
1797 RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
1798 RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
1799 CTYPE ("Method for doing Van der Waals");
1800 EETYPE("vdw-type", ir->vdwtype, evdw_names);
1801 EETYPE("vdw-modifier", ir->vdw_modifier, eintmod_names);
1802 CTYPE ("cut-off lengths");
1803 RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
1804 RTYPE ("rvdw", ir->rvdw, 1.0);
1805 CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
1806 EETYPE("DispCorr", ir->eDispCorr, edispc_names);
1807 CTYPE ("Extension of the potential lookup tables beyond the cut-off");
1808 RTYPE ("table-extension", ir->tabext, 1.0);
1809 CTYPE ("Separate tables between energy group pairs");
1810 STYPE ("energygrp-table", egptable, NULL);
1811 CTYPE ("Spacing for the PME/PPPM FFT grid");
1812 RTYPE ("fourierspacing", ir->fourier_spacing, 0.12);
1813 CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
1814 ITYPE ("fourier-nx", ir->nkx, 0);
1815 ITYPE ("fourier-ny", ir->nky, 0);
1816 ITYPE ("fourier-nz", ir->nkz, 0);
1817 CTYPE ("EWALD/PME/PPPM parameters");
1818 ITYPE ("pme-order", ir->pme_order, 4);
1819 RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
1820 EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
1821 RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
1822 EETYPE("optimize-fft", ir->bOptFFT, yesno_names);
1824 CCTYPE("IMPLICIT SOLVENT ALGORITHM");
1825 EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
1827 CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
1828 CTYPE ("Algorithm for calculating Born radii");
1829 EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
1830 CTYPE ("Frequency of calculating the Born radii inside rlist");
1831 ITYPE ("nstgbradii", ir->nstgbradii, 1);
1832 CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
1833 CTYPE ("between rlist and rgbradii is updated every nstlist steps");
1834 RTYPE ("rgbradii", ir->rgbradii, 1.0);
1835 CTYPE ("Dielectric coefficient of the implicit solvent");
1836 RTYPE ("gb-epsilon-solvent", ir->gb_epsilon_solvent, 80.0);
1837 CTYPE ("Salt concentration in M for Generalized Born models");
1838 RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
1839 CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
1840 RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
1841 RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
1842 RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
1843 RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
1844 EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
1845 CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
1846 CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
1847 RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
1849 /* Coupling stuff */
1850 CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
1851 CTYPE ("Temperature coupling");
1852 EETYPE("tcoupl", ir->etc, etcoupl_names);
1853 ITYPE ("nsttcouple", ir->nsttcouple, -1);
1854 ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
1855 EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
1856 CTYPE ("Groups to couple separately");
1857 STYPE ("tc-grps", tcgrps, NULL);
1858 CTYPE ("Time constant (ps) and reference temperature (K)");
1859 STYPE ("tau-t", tau_t, NULL);
1860 STYPE ("ref-t", ref_t, NULL);
1861 CTYPE ("pressure coupling");
1862 EETYPE("pcoupl", ir->epc, epcoupl_names);
1863 EETYPE("pcoupltype", ir->epct, epcoupltype_names);
1864 ITYPE ("nstpcouple", ir->nstpcouple, -1);
1865 CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
1866 RTYPE ("tau-p", ir->tau_p, 1.0);
1867 STYPE ("compressibility", dumstr[0], NULL);
1868 STYPE ("ref-p", dumstr[1], NULL);
1869 CTYPE ("Scaling of reference coordinates, No, All or COM");
1870 EETYPE ("refcoord-scaling", ir->refcoord_scaling, erefscaling_names);
1873 CCTYPE ("OPTIONS FOR QMMM calculations");
1874 EETYPE("QMMM", ir->bQMMM, yesno_names);
1875 CTYPE ("Groups treated Quantum Mechanically");
1876 STYPE ("QMMM-grps", QMMM, NULL);
1877 CTYPE ("QM method");
1878 STYPE("QMmethod", QMmethod, NULL);
1879 CTYPE ("QMMM scheme");
1880 EETYPE("QMMMscheme", ir->QMMMscheme, eQMMMscheme_names);
1881 CTYPE ("QM basisset");
1882 STYPE("QMbasis", QMbasis, NULL);
1883 CTYPE ("QM charge");
1884 STYPE ("QMcharge", QMcharge, NULL);
1885 CTYPE ("QM multiplicity");
1886 STYPE ("QMmult", QMmult, NULL);
1887 CTYPE ("Surface Hopping");
1888 STYPE ("SH", bSH, NULL);
1889 CTYPE ("CAS space options");
1890 STYPE ("CASorbitals", CASorbitals, NULL);
1891 STYPE ("CASelectrons", CASelectrons, NULL);
1892 STYPE ("SAon", SAon, NULL);
1893 STYPE ("SAoff", SAoff, NULL);
1894 STYPE ("SAsteps", SAsteps, NULL);
1895 CTYPE ("Scale factor for MM charges");
1896 RTYPE ("MMChargeScaleFactor", ir->scalefactor, 1.0);
1897 CTYPE ("Optimization of QM subsystem");
1898 STYPE ("bOPT", bOPT, NULL);
1899 STYPE ("bTS", bTS, NULL);
1901 /* Simulated annealing */
1902 CCTYPE("SIMULATED ANNEALING");
1903 CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
1904 STYPE ("annealing", anneal, NULL);
1905 CTYPE ("Number of time points to use for specifying annealing in each group");
1906 STYPE ("annealing-npoints", anneal_npoints, NULL);
1907 CTYPE ("List of times at the annealing points for each group");
1908 STYPE ("annealing-time", anneal_time, NULL);
1909 CTYPE ("Temp. at each annealing point, for each group.");
1910 STYPE ("annealing-temp", anneal_temp, NULL);
1913 CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
1914 EETYPE("gen-vel", opts->bGenVel, yesno_names);
1915 RTYPE ("gen-temp", opts->tempi, 300.0);
1916 ITYPE ("gen-seed", opts->seed, 173529);
1919 CCTYPE ("OPTIONS FOR BONDS");
1920 EETYPE("constraints", opts->nshake, constraints);
1921 CTYPE ("Type of constraint algorithm");
1922 EETYPE("constraint-algorithm", ir->eConstrAlg, econstr_names);
1923 CTYPE ("Do not constrain the start configuration");
1924 EETYPE("continuation", ir->bContinuation, yesno_names);
1925 CTYPE ("Use successive overrelaxation to reduce the number of shake iterations");
1926 EETYPE("Shake-SOR", ir->bShakeSOR, yesno_names);
1927 CTYPE ("Relative tolerance of shake");
1928 RTYPE ("shake-tol", ir->shake_tol, 0.0001);
1929 CTYPE ("Highest order in the expansion of the constraint coupling matrix");
1930 ITYPE ("lincs-order", ir->nProjOrder, 4);
1931 CTYPE ("Number of iterations in the final step of LINCS. 1 is fine for");
1932 CTYPE ("normal simulations, but use 2 to conserve energy in NVE runs.");
1933 CTYPE ("For energy minimization with constraints it should be 4 to 8.");
1934 ITYPE ("lincs-iter", ir->nLincsIter, 1);
1935 CTYPE ("Lincs will write a warning to the stderr if in one step a bond");
1936 CTYPE ("rotates over more degrees than");
1937 RTYPE ("lincs-warnangle", ir->LincsWarnAngle, 30.0);
1938 CTYPE ("Convert harmonic bonds to morse potentials");
1939 EETYPE("morse", opts->bMorse, yesno_names);
1941 /* Energy group exclusions */
1942 CCTYPE ("ENERGY GROUP EXCLUSIONS");
1943 CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
1944 STYPE ("energygrp-excl", egpexcl, NULL);
1948 CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
1949 ITYPE ("nwall", ir->nwall, 0);
1950 EETYPE("wall-type", ir->wall_type, ewt_names);
1951 RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
1952 STYPE ("wall-atomtype", wall_atomtype, NULL);
1953 STYPE ("wall-density", wall_density, NULL);
1954 RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
1957 CCTYPE("COM PULLING");
1958 CTYPE("Pull type: no, umbrella, constraint or constant-force");
1959 EETYPE("pull", ir->ePull, epull_names);
1960 if (ir->ePull != epullNO)
1963 pull_grp = read_pullparams(&ninp, &inp, ir->pull, &opts->pull_start, wi);
1966 /* Enforced rotation */
1967 CCTYPE("ENFORCED ROTATION");
1968 CTYPE("Enforced rotation: No or Yes");
1969 EETYPE("rotation", ir->bRot, yesno_names);
1973 rot_grp = read_rotparams(&ninp, &inp, ir->rot, wi);
1977 CCTYPE("NMR refinement stuff");
1978 CTYPE ("Distance restraints type: No, Simple or Ensemble");
1979 EETYPE("disre", ir->eDisre, edisre_names);
1980 CTYPE ("Force weighting of pairs in one distance restraint: Conservative or Equal");
1981 EETYPE("disre-weighting", ir->eDisreWeighting, edisreweighting_names);
1982 CTYPE ("Use sqrt of the time averaged times the instantaneous violation");
1983 EETYPE("disre-mixed", ir->bDisreMixed, yesno_names);
1984 RTYPE ("disre-fc", ir->dr_fc, 1000.0);
1985 RTYPE ("disre-tau", ir->dr_tau, 0.0);
1986 CTYPE ("Output frequency for pair distances to energy file");
1987 ITYPE ("nstdisreout", ir->nstdisreout, 100);
1988 CTYPE ("Orientation restraints: No or Yes");
1989 EETYPE("orire", opts->bOrire, yesno_names);
1990 CTYPE ("Orientation restraints force constant and tau for time averaging");
1991 RTYPE ("orire-fc", ir->orires_fc, 0.0);
1992 RTYPE ("orire-tau", ir->orires_tau, 0.0);
1993 STYPE ("orire-fitgrp", orirefitgrp, NULL);
1994 CTYPE ("Output frequency for trace(SD) and S to energy file");
1995 ITYPE ("nstorireout", ir->nstorireout, 100);
1997 /* free energy variables */
1998 CCTYPE ("Free energy variables");
1999 EETYPE("free-energy", ir->efep, efep_names);
2000 STYPE ("couple-moltype", couple_moltype, NULL);
2001 EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
2002 EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
2003 EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
2005 RTYPE ("init-lambda", fep->init_lambda, -1); /* start with -1 so
2007 it was not entered */
2008 ITYPE ("init-lambda-state", fep->init_fep_state, -1);
2009 RTYPE ("delta-lambda", fep->delta_lambda, 0.0);
2010 ITYPE ("nstdhdl", fep->nstdhdl, 50);
2011 STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
2012 STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
2013 STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
2014 STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
2015 STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
2016 STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
2017 STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
2018 ITYPE ("calc-lambda-neighbors", fep->lambda_neighbors, 1);
2019 STYPE ("init-lambda-weights", lambda_weights, NULL);
2020 EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
2021 RTYPE ("sc-alpha", fep->sc_alpha, 0.0);
2022 ITYPE ("sc-power", fep->sc_power, 1);
2023 RTYPE ("sc-r-power", fep->sc_r_power, 6.0);
2024 RTYPE ("sc-sigma", fep->sc_sigma, 0.3);
2025 EETYPE("sc-coul", fep->bScCoul, yesno_names);
2026 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2027 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2028 EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
2029 separate_dhdl_file_names);
2030 EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
2031 ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
2032 RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
2034 /* Non-equilibrium MD stuff */
2035 CCTYPE("Non-equilibrium MD stuff");
2036 STYPE ("acc-grps", accgrps, NULL);
2037 STYPE ("accelerate", acc, NULL);
2038 STYPE ("freezegrps", freeze, NULL);
2039 STYPE ("freezedim", frdim, NULL);
2040 RTYPE ("cos-acceleration", ir->cos_accel, 0);
2041 STYPE ("deform", deform, NULL);
2043 /* simulated tempering variables */
2044 CCTYPE("simulated tempering variables");
2045 EETYPE("simulated-tempering", ir->bSimTemp, yesno_names);
2046 EETYPE("simulated-tempering-scaling", ir->simtempvals->eSimTempScale, esimtemp_names);
2047 RTYPE("sim-temp-low", ir->simtempvals->simtemp_low, 300.0);
2048 RTYPE("sim-temp-high", ir->simtempvals->simtemp_high, 300.0);
2050 /* expanded ensemble variables */
2051 if (ir->efep == efepEXPANDED || ir->bSimTemp)
2053 read_expandedparams(&ninp, &inp, expand, wi);
2056 /* Electric fields */
2057 CCTYPE("Electric fields");
2058 CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
2059 CTYPE ("and a phase angle (real)");
2060 STYPE ("E-x", efield_x, NULL);
2061 STYPE ("E-xt", efield_xt, NULL);
2062 STYPE ("E-y", efield_y, NULL);
2063 STYPE ("E-yt", efield_yt, NULL);
2064 STYPE ("E-z", efield_z, NULL);
2065 STYPE ("E-zt", efield_zt, NULL);
2067 /* AdResS defined thingies */
2068 CCTYPE ("AdResS parameters");
2069 EETYPE("adress", ir->bAdress, yesno_names);
2072 snew(ir->adress, 1);
2073 read_adressparams(&ninp, &inp, ir->adress, wi);
2076 /* User defined thingies */
2077 CCTYPE ("User defined thingies");
2078 STYPE ("user1-grps", user1, NULL);
2079 STYPE ("user2-grps", user2, NULL);
2080 ITYPE ("userint1", ir->userint1, 0);
2081 ITYPE ("userint2", ir->userint2, 0);
2082 ITYPE ("userint3", ir->userint3, 0);
2083 ITYPE ("userint4", ir->userint4, 0);
2084 RTYPE ("userreal1", ir->userreal1, 0);
2085 RTYPE ("userreal2", ir->userreal2, 0);
2086 RTYPE ("userreal3", ir->userreal3, 0);
2087 RTYPE ("userreal4", ir->userreal4, 0);
2090 write_inpfile(mdparout, ninp, inp, FALSE, wi);
2091 for (i = 0; (i < ninp); i++)
2094 sfree(inp[i].value);
2098 /* Process options if necessary */
2099 for (m = 0; m < 2; m++)
2101 for (i = 0; i < 2*DIM; i++)
2110 if (sscanf(dumstr[m], "%lf", &(dumdub[m][XX])) != 1)
2112 warning_error(wi, "Pressure coupling not enough values (I need 1)");
2114 dumdub[m][YY] = dumdub[m][ZZ] = dumdub[m][XX];
2116 case epctSEMIISOTROPIC:
2117 case epctSURFACETENSION:
2118 if (sscanf(dumstr[m], "%lf%lf",
2119 &(dumdub[m][XX]), &(dumdub[m][ZZ])) != 2)
2121 warning_error(wi, "Pressure coupling not enough values (I need 2)");
2123 dumdub[m][YY] = dumdub[m][XX];
2125 case epctANISOTROPIC:
2126 if (sscanf(dumstr[m], "%lf%lf%lf%lf%lf%lf",
2127 &(dumdub[m][XX]), &(dumdub[m][YY]), &(dumdub[m][ZZ]),
2128 &(dumdub[m][3]), &(dumdub[m][4]), &(dumdub[m][5])) != 6)
2130 warning_error(wi, "Pressure coupling not enough values (I need 6)");
2134 gmx_fatal(FARGS, "Pressure coupling type %s not implemented yet",
2135 epcoupltype_names[ir->epct]);
2139 clear_mat(ir->ref_p);
2140 clear_mat(ir->compress);
2141 for (i = 0; i < DIM; i++)
2143 ir->ref_p[i][i] = dumdub[1][i];
2144 ir->compress[i][i] = dumdub[0][i];
2146 if (ir->epct == epctANISOTROPIC)
2148 ir->ref_p[XX][YY] = dumdub[1][3];
2149 ir->ref_p[XX][ZZ] = dumdub[1][4];
2150 ir->ref_p[YY][ZZ] = dumdub[1][5];
2151 if (ir->ref_p[XX][YY] != 0 && ir->ref_p[XX][ZZ] != 0 && ir->ref_p[YY][ZZ] != 0)
2153 warning(wi, "All off-diagonal reference pressures are non-zero. Are you sure you want to apply a threefold shear stress?\n");
2155 ir->compress[XX][YY] = dumdub[0][3];
2156 ir->compress[XX][ZZ] = dumdub[0][4];
2157 ir->compress[YY][ZZ] = dumdub[0][5];
2158 for (i = 0; i < DIM; i++)
2160 for (m = 0; m < i; m++)
2162 ir->ref_p[i][m] = ir->ref_p[m][i];
2163 ir->compress[i][m] = ir->compress[m][i];
2168 if (ir->comm_mode == ecmNO)
2173 opts->couple_moltype = NULL;
2174 if (strlen(couple_moltype) > 0)
2176 if (ir->efep != efepNO)
2178 opts->couple_moltype = strdup(couple_moltype);
2179 if (opts->couple_lam0 == opts->couple_lam1)
2181 warning(wi, "The lambda=0 and lambda=1 states for coupling are identical");
2183 if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
2184 opts->couple_lam1 == ecouplamNONE))
2186 warning(wi, "For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
2191 warning(wi, "Can not couple a molecule with free_energy = no");
2194 /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
2195 if (ir->efep != efepNO)
2197 if (fep->delta_lambda > 0)
2199 ir->efep = efepSLOWGROWTH;
2205 fep->bPrintEnergy = TRUE;
2206 /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
2207 if the temperature is changing. */
2210 if ((ir->efep != efepNO) || ir->bSimTemp)
2212 ir->bExpanded = FALSE;
2213 if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
2215 ir->bExpanded = TRUE;
2217 do_fep_params(ir, fep_lambda, lambda_weights);
2218 if (ir->bSimTemp) /* done after fep params */
2220 do_simtemp_params(ir);
2225 ir->fepvals->n_lambda = 0;
2228 /* WALL PARAMETERS */
2230 do_wall_params(ir, wall_atomtype, wall_density, opts);
2232 /* ORIENTATION RESTRAINT PARAMETERS */
2234 if (opts->bOrire && str_nelem(orirefitgrp, MAXPTR, NULL) != 1)
2236 warning_error(wi, "ERROR: Need one orientation restraint fit group\n");
2239 /* DEFORMATION PARAMETERS */
2241 clear_mat(ir->deform);
2242 for (i = 0; i < 6; i++)
2246 m = sscanf(deform, "%lf %lf %lf %lf %lf %lf",
2247 &(dumdub[0][0]), &(dumdub[0][1]), &(dumdub[0][2]),
2248 &(dumdub[0][3]), &(dumdub[0][4]), &(dumdub[0][5]));
2249 for (i = 0; i < 3; i++)
2251 ir->deform[i][i] = dumdub[0][i];
2253 ir->deform[YY][XX] = dumdub[0][3];
2254 ir->deform[ZZ][XX] = dumdub[0][4];
2255 ir->deform[ZZ][YY] = dumdub[0][5];
2256 if (ir->epc != epcNO)
2258 for (i = 0; i < 3; i++)
2260 for (j = 0; j <= i; j++)
2262 if (ir->deform[i][j] != 0 && ir->compress[i][j] != 0)
2264 warning_error(wi, "A box element has deform set and compressibility > 0");
2268 for (i = 0; i < 3; i++)
2270 for (j = 0; j < i; j++)
2272 if (ir->deform[i][j] != 0)
2274 for (m = j; m < DIM; m++)
2276 if (ir->compress[m][j] != 0)
2278 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.");
2279 warning(wi, warn_buf);
2291 static int search_QMstring(char *s, int ng, const char *gn[])
2293 /* same as normal search_string, but this one searches QM strings */
2296 for (i = 0; (i < ng); i++)
2298 if (gmx_strcasecmp(s, gn[i]) == 0)
2304 gmx_fatal(FARGS, "this QM method or basisset (%s) is not implemented\n!", s);
2308 } /* search_QMstring */
2311 int search_string(char *s, int ng, char *gn[])
2315 for (i = 0; (i < ng); i++)
2317 if (gmx_strcasecmp(s, gn[i]) == 0)
2324 "Group %s referenced in the .mdp file was not found in the index file.\n"
2325 "Group names must match either [moleculetype] names or custom index group\n"
2326 "names, in which case you must supply an index file to the '-n' option\n"
2333 static gmx_bool do_numbering(int natoms, gmx_groups_t *groups, int ng, char *ptrs[],
2334 t_blocka *block, char *gnames[],
2335 int gtype, int restnm,
2336 int grptp, gmx_bool bVerbose,
2339 unsigned short *cbuf;
2340 t_grps *grps = &(groups->grps[gtype]);
2341 int i, j, gid, aj, ognr, ntot = 0;
2344 char warn_buf[STRLEN];
2348 fprintf(debug, "Starting numbering %d groups of type %d\n", ng, gtype);
2351 title = gtypes[gtype];
2354 /* Mark all id's as not set */
2355 for (i = 0; (i < natoms); i++)
2360 snew(grps->nm_ind, ng+1); /* +1 for possible rest group */
2361 for (i = 0; (i < ng); i++)
2363 /* Lookup the group name in the block structure */
2364 gid = search_string(ptrs[i], block->nr, gnames);
2365 if ((grptp != egrptpONE) || (i == 0))
2367 grps->nm_ind[grps->nr++] = gid;
2371 fprintf(debug, "Found gid %d for group %s\n", gid, ptrs[i]);
2374 /* Now go over the atoms in the group */
2375 for (j = block->index[gid]; (j < block->index[gid+1]); j++)
2380 /* Range checking */
2381 if ((aj < 0) || (aj >= natoms))
2383 gmx_fatal(FARGS, "Invalid atom number %d in indexfile", aj);
2385 /* Lookup up the old group number */
2389 gmx_fatal(FARGS, "Atom %d in multiple %s groups (%d and %d)",
2390 aj+1, title, ognr+1, i+1);
2394 /* Store the group number in buffer */
2395 if (grptp == egrptpONE)
2408 /* Now check whether we have done all atoms */
2412 if (grptp == egrptpALL)
2414 gmx_fatal(FARGS, "%d atoms are not part of any of the %s groups",
2415 natoms-ntot, title);
2417 else if (grptp == egrptpPART)
2419 sprintf(warn_buf, "%d atoms are not part of any of the %s groups",
2420 natoms-ntot, title);
2421 warning_note(wi, warn_buf);
2423 /* Assign all atoms currently unassigned to a rest group */
2424 for (j = 0; (j < natoms); j++)
2426 if (cbuf[j] == NOGID)
2432 if (grptp != egrptpPART)
2437 "Making dummy/rest group for %s containing %d elements\n",
2438 title, natoms-ntot);
2440 /* Add group name "rest" */
2441 grps->nm_ind[grps->nr] = restnm;
2443 /* Assign the rest name to all atoms not currently assigned to a group */
2444 for (j = 0; (j < natoms); j++)
2446 if (cbuf[j] == NOGID)
2455 if (grps->nr == 1 && (ntot == 0 || ntot == natoms))
2457 /* All atoms are part of one (or no) group, no index required */
2458 groups->ngrpnr[gtype] = 0;
2459 groups->grpnr[gtype] = NULL;
2463 groups->ngrpnr[gtype] = natoms;
2464 snew(groups->grpnr[gtype], natoms);
2465 for (j = 0; (j < natoms); j++)
2467 groups->grpnr[gtype][j] = cbuf[j];
2473 return (bRest && grptp == egrptpPART);
2476 static void calc_nrdf(gmx_mtop_t *mtop, t_inputrec *ir, char **gnames)
2479 gmx_groups_t *groups;
2481 int natoms, ai, aj, i, j, d, g, imin, jmin, nc;
2483 int *nrdf2, *na_vcm, na_tot;
2484 double *nrdf_tc, *nrdf_vcm, nrdf_uc, n_sub = 0;
2485 gmx_mtop_atomloop_all_t aloop;
2487 int mb, mol, ftype, as;
2488 gmx_molblock_t *molb;
2489 gmx_moltype_t *molt;
2492 * First calc 3xnr-atoms for each group
2493 * then subtract half a degree of freedom for each constraint
2495 * Only atoms and nuclei contribute to the degrees of freedom...
2500 groups = &mtop->groups;
2501 natoms = mtop->natoms;
2503 /* Allocate one more for a possible rest group */
2504 /* We need to sum degrees of freedom into doubles,
2505 * since floats give too low nrdf's above 3 million atoms.
2507 snew(nrdf_tc, groups->grps[egcTC].nr+1);
2508 snew(nrdf_vcm, groups->grps[egcVCM].nr+1);
2509 snew(na_vcm, groups->grps[egcVCM].nr+1);
2511 for (i = 0; i < groups->grps[egcTC].nr; i++)
2515 for (i = 0; i < groups->grps[egcVCM].nr+1; i++)
2520 snew(nrdf2, natoms);
2521 aloop = gmx_mtop_atomloop_all_init(mtop);
2522 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
2525 if (atom->ptype == eptAtom || atom->ptype == eptNucleus)
2527 g = ggrpnr(groups, egcFREEZE, i);
2528 /* Double count nrdf for particle i */
2529 for (d = 0; d < DIM; d++)
2531 if (opts->nFreeze[g][d] == 0)
2536 nrdf_tc [ggrpnr(groups, egcTC, i)] += 0.5*nrdf2[i];
2537 nrdf_vcm[ggrpnr(groups, egcVCM, i)] += 0.5*nrdf2[i];
2542 for (mb = 0; mb < mtop->nmolblock; mb++)
2544 molb = &mtop->molblock[mb];
2545 molt = &mtop->moltype[molb->type];
2546 atom = molt->atoms.atom;
2547 for (mol = 0; mol < molb->nmol; mol++)
2549 for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
2551 ia = molt->ilist[ftype].iatoms;
2552 for (i = 0; i < molt->ilist[ftype].nr; )
2554 /* Subtract degrees of freedom for the constraints,
2555 * if the particles still have degrees of freedom left.
2556 * If one of the particles is a vsite or a shell, then all
2557 * constraint motion will go there, but since they do not
2558 * contribute to the constraints the degrees of freedom do not
2563 if (((atom[ia[1]].ptype == eptNucleus) ||
2564 (atom[ia[1]].ptype == eptAtom)) &&
2565 ((atom[ia[2]].ptype == eptNucleus) ||
2566 (atom[ia[2]].ptype == eptAtom)))
2584 imin = min(imin, nrdf2[ai]);
2585 jmin = min(jmin, nrdf2[aj]);
2588 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2589 nrdf_tc [ggrpnr(groups, egcTC, aj)] -= 0.5*jmin;
2590 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2591 nrdf_vcm[ggrpnr(groups, egcVCM, aj)] -= 0.5*jmin;
2593 ia += interaction_function[ftype].nratoms+1;
2594 i += interaction_function[ftype].nratoms+1;
2597 ia = molt->ilist[F_SETTLE].iatoms;
2598 for (i = 0; i < molt->ilist[F_SETTLE].nr; )
2600 /* Subtract 1 dof from every atom in the SETTLE */
2601 for (j = 0; j < 3; j++)
2604 imin = min(2, nrdf2[ai]);
2606 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2607 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2612 as += molt->atoms.nr;
2616 if (ir->ePull == epullCONSTRAINT)
2618 /* Correct nrdf for the COM constraints.
2619 * We correct using the TC and VCM group of the first atom
2620 * in the reference and pull group. If atoms in one pull group
2621 * belong to different TC or VCM groups it is anyhow difficult
2622 * to determine the optimal nrdf assignment.
2625 if (pull->eGeom == epullgPOS)
2628 for (i = 0; i < DIM; i++)
2640 for (i = 0; i < pull->ngrp; i++)
2643 if (pull->grp[0].nat > 0)
2645 /* Subtract 1/2 dof from the reference group */
2646 ai = pull->grp[0].ind[0];
2647 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] > 1)
2649 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5;
2650 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5;
2654 /* Subtract 1/2 dof from the pulled group */
2655 ai = pull->grp[1+i].ind[0];
2656 nrdf_tc [ggrpnr(groups, egcTC, ai)] -= 0.5*imin;
2657 nrdf_vcm[ggrpnr(groups, egcVCM, ai)] -= 0.5*imin;
2658 if (nrdf_tc[ggrpnr(groups, egcTC, ai)] < 0)
2660 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)]]);
2665 if (ir->nstcomm != 0)
2667 /* Subtract 3 from the number of degrees of freedom in each vcm group
2668 * when com translation is removed and 6 when rotation is removed
2671 switch (ir->comm_mode)
2674 n_sub = ndof_com(ir);
2681 gmx_incons("Checking comm_mode");
2684 for (i = 0; i < groups->grps[egcTC].nr; i++)
2686 /* Count the number of atoms of TC group i for every VCM group */
2687 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2692 for (ai = 0; ai < natoms; ai++)
2694 if (ggrpnr(groups, egcTC, ai) == i)
2696 na_vcm[ggrpnr(groups, egcVCM, ai)]++;
2700 /* Correct for VCM removal according to the fraction of each VCM
2701 * group present in this TC group.
2703 nrdf_uc = nrdf_tc[i];
2706 fprintf(debug, "T-group[%d] nrdf_uc = %g, n_sub = %g\n",
2710 for (j = 0; j < groups->grps[egcVCM].nr+1; j++)
2712 if (nrdf_vcm[j] > n_sub)
2714 nrdf_tc[i] += nrdf_uc*((double)na_vcm[j]/(double)na_tot)*
2715 (nrdf_vcm[j] - n_sub)/nrdf_vcm[j];
2719 fprintf(debug, " nrdf_vcm[%d] = %g, nrdf = %g\n",
2720 j, nrdf_vcm[j], nrdf_tc[i]);
2725 for (i = 0; (i < groups->grps[egcTC].nr); i++)
2727 opts->nrdf[i] = nrdf_tc[i];
2728 if (opts->nrdf[i] < 0)
2733 "Number of degrees of freedom in T-Coupling group %s is %.2f\n",
2734 gnames[groups->grps[egcTC].nm_ind[i]], opts->nrdf[i]);
2743 static void decode_cos(char *s, t_cosines *cosine, gmx_bool bTime)
2746 char format[STRLEN], f1[STRLEN];
2758 sscanf(t, "%d", &(cosine->n));
2765 snew(cosine->a, cosine->n);
2766 snew(cosine->phi, cosine->n);
2768 sprintf(format, "%%*d");
2769 for (i = 0; (i < cosine->n); i++)
2772 strcat(f1, "%lf%lf");
2773 if (sscanf(t, f1, &a, &phi) < 2)
2775 gmx_fatal(FARGS, "Invalid input for electric field shift: '%s'", t);
2778 cosine->phi[i] = phi;
2779 strcat(format, "%*lf%*lf");
2786 static gmx_bool do_egp_flag(t_inputrec *ir, gmx_groups_t *groups,
2787 const char *option, const char *val, int flag)
2789 /* The maximum number of energy group pairs would be MAXPTR*(MAXPTR+1)/2.
2790 * But since this is much larger than STRLEN, such a line can not be parsed.
2791 * The real maximum is the number of names that fit in a string: STRLEN/2.
2793 #define EGP_MAX (STRLEN/2)
2794 int nelem, i, j, k, nr;
2795 char *names[EGP_MAX];
2799 gnames = groups->grpname;
2801 nelem = str_nelem(val, EGP_MAX, names);
2804 gmx_fatal(FARGS, "The number of groups for %s is odd", option);
2806 nr = groups->grps[egcENER].nr;
2808 for (i = 0; i < nelem/2; i++)
2812 gmx_strcasecmp(names[2*i], *(gnames[groups->grps[egcENER].nm_ind[j]])))
2818 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2819 names[2*i], option);
2823 gmx_strcasecmp(names[2*i+1], *(gnames[groups->grps[egcENER].nm_ind[k]])))
2829 gmx_fatal(FARGS, "%s in %s is not an energy group\n",
2830 names[2*i+1], option);
2832 if ((j < nr) && (k < nr))
2834 ir->opts.egp_flags[nr*j+k] |= flag;
2835 ir->opts.egp_flags[nr*k+j] |= flag;
2843 void do_index(const char* mdparin, const char *ndx,
2846 t_inputrec *ir, rvec *v,
2850 gmx_groups_t *groups;
2854 char warnbuf[STRLEN], **gnames;
2855 int nr, ntcg, ntau_t, nref_t, nacc, nofg, nSA, nSA_points, nSA_time, nSA_temp;
2858 int nacg, nfreeze, nfrdim, nenergy, nvcm, nuser;
2859 char *ptr1[MAXPTR], *ptr2[MAXPTR], *ptr3[MAXPTR];
2860 int i, j, k, restnm;
2862 gmx_bool bExcl, bTable, bSetTCpar, bAnneal, bRest;
2863 int nQMmethod, nQMbasis, nQMcharge, nQMmult, nbSH, nCASorb, nCASelec,
2864 nSAon, nSAoff, nSAsteps, nQMg, nbOPT, nbTS;
2865 char warn_buf[STRLEN];
2869 fprintf(stderr, "processing index file...\n");
2875 snew(grps->index, 1);
2877 atoms_all = gmx_mtop_global_atoms(mtop);
2878 analyse(&atoms_all, grps, &gnames, FALSE, TRUE);
2879 free_t_atoms(&atoms_all, FALSE);
2883 grps = init_index(ndx, &gnames);
2886 groups = &mtop->groups;
2887 natoms = mtop->natoms;
2888 symtab = &mtop->symtab;
2890 snew(groups->grpname, grps->nr+1);
2892 for (i = 0; (i < grps->nr); i++)
2894 groups->grpname[i] = put_symtab(symtab, gnames[i]);
2896 groups->grpname[i] = put_symtab(symtab, "rest");
2898 srenew(gnames, grps->nr+1);
2899 gnames[restnm] = *(groups->grpname[i]);
2900 groups->ngrpname = grps->nr+1;
2902 set_warning_line(wi, mdparin, -1);
2904 ntau_t = str_nelem(tau_t, MAXPTR, ptr1);
2905 nref_t = str_nelem(ref_t, MAXPTR, ptr2);
2906 ntcg = str_nelem(tcgrps, MAXPTR, ptr3);
2907 if ((ntau_t != ntcg) || (nref_t != ntcg))
2909 gmx_fatal(FARGS, "Invalid T coupling input: %d groups, %d ref-t values and "
2910 "%d tau-t values", ntcg, nref_t, ntau_t);
2913 bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI == eiBD || EI_TPI(ir->eI));
2914 do_numbering(natoms, groups, ntcg, ptr3, grps, gnames, egcTC,
2915 restnm, bSetTCpar ? egrptpALL : egrptpALL_GENREST, bVerbose, wi);
2916 nr = groups->grps[egcTC].nr;
2918 snew(ir->opts.nrdf, nr);
2919 snew(ir->opts.tau_t, nr);
2920 snew(ir->opts.ref_t, nr);
2921 if (ir->eI == eiBD && ir->bd_fric == 0)
2923 fprintf(stderr, "bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
2930 gmx_fatal(FARGS, "Not enough ref-t and tau-t values!");
2934 for (i = 0; (i < nr); i++)
2936 ir->opts.tau_t[i] = strtod(ptr1[i], NULL);
2937 if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
2939 sprintf(warn_buf, "With integrator %s tau-t should be larger than 0", ei_names[ir->eI]);
2940 warning_error(wi, warn_buf);
2943 if (ir->etc != etcVRESCALE && ir->opts.tau_t[i] == 0)
2945 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.");
2948 if (ir->opts.tau_t[i] >= 0)
2950 tau_min = min(tau_min, ir->opts.tau_t[i]);
2953 if (ir->etc != etcNO && ir->nsttcouple == -1)
2955 ir->nsttcouple = ir_optimal_nsttcouple(ir);
2960 if ((ir->etc == etcNOSEHOOVER) && (ir->epc == epcBERENDSEN))
2962 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");
2964 if ((ir->epc == epcMTTK) && (ir->etc > etcNO))
2966 if (ir->nstpcouple != ir->nsttcouple)
2968 int mincouple = min(ir->nstpcouple, ir->nsttcouple);
2969 ir->nstpcouple = ir->nsttcouple = mincouple;
2970 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);
2971 warning_note(wi, warn_buf);
2975 /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
2976 primarily for testing purposes, and does not work with temperature coupling other than 1 */
2978 if (ETC_ANDERSEN(ir->etc))
2980 if (ir->nsttcouple != 1)
2983 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");
2984 warning_note(wi, warn_buf);
2987 nstcmin = tcouple_min_integration_steps(ir->etc);
2990 if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
2992 sprintf(warn_buf, "For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
2993 ETCOUPLTYPE(ir->etc),
2995 ir->nsttcouple*ir->delta_t);
2996 warning(wi, warn_buf);
2999 for (i = 0; (i < nr); i++)
3001 ir->opts.ref_t[i] = strtod(ptr2[i], NULL);
3002 if (ir->opts.ref_t[i] < 0)
3004 gmx_fatal(FARGS, "ref-t for group %d negative", i);
3007 /* set the lambda mc temperature to the md integrator temperature (which should be defined
3008 if we are in this conditional) if mc_temp is negative */
3009 if (ir->expandedvals->mc_temp < 0)
3011 ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
3015 /* Simulated annealing for each group. There are nr groups */
3016 nSA = str_nelem(anneal, MAXPTR, ptr1);
3017 if (nSA == 1 && (ptr1[0][0] == 'n' || ptr1[0][0] == 'N'))
3021 if (nSA > 0 && nSA != nr)
3023 gmx_fatal(FARGS, "Not enough annealing values: %d (for %d groups)\n", nSA, nr);
3027 snew(ir->opts.annealing, nr);
3028 snew(ir->opts.anneal_npoints, nr);
3029 snew(ir->opts.anneal_time, nr);
3030 snew(ir->opts.anneal_temp, nr);
3031 for (i = 0; i < nr; i++)
3033 ir->opts.annealing[i] = eannNO;
3034 ir->opts.anneal_npoints[i] = 0;
3035 ir->opts.anneal_time[i] = NULL;
3036 ir->opts.anneal_temp[i] = NULL;
3041 for (i = 0; i < nr; i++)
3043 if (ptr1[i][0] == 'n' || ptr1[i][0] == 'N')
3045 ir->opts.annealing[i] = eannNO;
3047 else if (ptr1[i][0] == 's' || ptr1[i][0] == 'S')
3049 ir->opts.annealing[i] = eannSINGLE;
3052 else if (ptr1[i][0] == 'p' || ptr1[i][0] == 'P')
3054 ir->opts.annealing[i] = eannPERIODIC;
3060 /* Read the other fields too */
3061 nSA_points = str_nelem(anneal_npoints, MAXPTR, ptr1);
3062 if (nSA_points != nSA)
3064 gmx_fatal(FARGS, "Found %d annealing-npoints values for %d groups\n", nSA_points, nSA);
3066 for (k = 0, i = 0; i < nr; i++)
3068 ir->opts.anneal_npoints[i] = strtol(ptr1[i], NULL, 10);
3069 if (ir->opts.anneal_npoints[i] == 1)
3071 gmx_fatal(FARGS, "Please specify at least a start and an end point for annealing\n");
3073 snew(ir->opts.anneal_time[i], ir->opts.anneal_npoints[i]);
3074 snew(ir->opts.anneal_temp[i], ir->opts.anneal_npoints[i]);
3075 k += ir->opts.anneal_npoints[i];
3078 nSA_time = str_nelem(anneal_time, MAXPTR, ptr1);
3081 gmx_fatal(FARGS, "Found %d annealing-time values, wanter %d\n", nSA_time, k);
3083 nSA_temp = str_nelem(anneal_temp, MAXPTR, ptr2);
3086 gmx_fatal(FARGS, "Found %d annealing-temp values, wanted %d\n", nSA_temp, k);
3089 for (i = 0, k = 0; i < nr; i++)
3092 for (j = 0; j < ir->opts.anneal_npoints[i]; j++)
3094 ir->opts.anneal_time[i][j] = strtod(ptr1[k], NULL);
3095 ir->opts.anneal_temp[i][j] = strtod(ptr2[k], NULL);
3098 if (ir->opts.anneal_time[i][0] > (ir->init_t+GMX_REAL_EPS))
3100 gmx_fatal(FARGS, "First time point for annealing > init_t.\n");
3106 if (ir->opts.anneal_time[i][j] < ir->opts.anneal_time[i][j-1])
3108 gmx_fatal(FARGS, "Annealing timepoints out of order: t=%f comes after t=%f\n",
3109 ir->opts.anneal_time[i][j], ir->opts.anneal_time[i][j-1]);
3112 if (ir->opts.anneal_temp[i][j] < 0)
3114 gmx_fatal(FARGS, "Found negative temperature in annealing: %f\n", ir->opts.anneal_temp[i][j]);
3119 /* Print out some summary information, to make sure we got it right */
3120 for (i = 0, k = 0; i < nr; i++)
3122 if (ir->opts.annealing[i] != eannNO)
3124 j = groups->grps[egcTC].nm_ind[i];
3125 fprintf(stderr, "Simulated annealing for group %s: %s, %d timepoints\n",
3126 *(groups->grpname[j]), eann_names[ir->opts.annealing[i]],
3127 ir->opts.anneal_npoints[i]);
3128 fprintf(stderr, "Time (ps) Temperature (K)\n");
3129 /* All terms except the last one */
3130 for (j = 0; j < (ir->opts.anneal_npoints[i]-1); j++)
3132 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3135 /* Finally the last one */
3136 j = ir->opts.anneal_npoints[i]-1;
3137 if (ir->opts.annealing[i] == eannSINGLE)
3139 fprintf(stderr, "%9.1f- %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3143 fprintf(stderr, "%9.1f %5.1f\n", ir->opts.anneal_time[i][j], ir->opts.anneal_temp[i][j]);
3144 if (fabs(ir->opts.anneal_temp[i][j]-ir->opts.anneal_temp[i][0]) > GMX_REAL_EPS)
3146 warning_note(wi, "There is a temperature jump when your annealing loops back.\n");
3155 if (ir->ePull != epullNO)
3157 make_pull_groups(ir->pull, pull_grp, grps, gnames);
3162 make_rotation_groups(ir->rot, rot_grp, grps, gnames);
3165 nacc = str_nelem(acc, MAXPTR, ptr1);
3166 nacg = str_nelem(accgrps, MAXPTR, ptr2);
3167 if (nacg*DIM != nacc)
3169 gmx_fatal(FARGS, "Invalid Acceleration input: %d groups and %d acc. values",
3172 do_numbering(natoms, groups, nacg, ptr2, grps, gnames, egcACC,
3173 restnm, egrptpALL_GENREST, bVerbose, wi);
3174 nr = groups->grps[egcACC].nr;
3175 snew(ir->opts.acc, nr);
3176 ir->opts.ngacc = nr;
3178 for (i = k = 0; (i < nacg); i++)
3180 for (j = 0; (j < DIM); j++, k++)
3182 ir->opts.acc[i][j] = strtod(ptr1[k], NULL);
3185 for (; (i < nr); i++)
3187 for (j = 0; (j < DIM); j++)
3189 ir->opts.acc[i][j] = 0;
3193 nfrdim = str_nelem(frdim, MAXPTR, ptr1);
3194 nfreeze = str_nelem(freeze, MAXPTR, ptr2);
3195 if (nfrdim != DIM*nfreeze)
3197 gmx_fatal(FARGS, "Invalid Freezing input: %d groups and %d freeze values",
3200 do_numbering(natoms, groups, nfreeze, ptr2, grps, gnames, egcFREEZE,
3201 restnm, egrptpALL_GENREST, bVerbose, wi);
3202 nr = groups->grps[egcFREEZE].nr;
3203 ir->opts.ngfrz = nr;
3204 snew(ir->opts.nFreeze, nr);
3205 for (i = k = 0; (i < nfreeze); i++)
3207 for (j = 0; (j < DIM); j++, k++)
3209 ir->opts.nFreeze[i][j] = (gmx_strncasecmp(ptr1[k], "Y", 1) == 0);
3210 if (!ir->opts.nFreeze[i][j])
3212 if (gmx_strncasecmp(ptr1[k], "N", 1) != 0)
3214 sprintf(warnbuf, "Please use Y(ES) or N(O) for freezedim only "
3215 "(not %s)", ptr1[k]);
3216 warning(wi, warn_buf);
3221 for (; (i < nr); i++)
3223 for (j = 0; (j < DIM); j++)
3225 ir->opts.nFreeze[i][j] = 0;
3229 nenergy = str_nelem(energy, MAXPTR, ptr1);
3230 do_numbering(natoms, groups, nenergy, ptr1, grps, gnames, egcENER,
3231 restnm, egrptpALL_GENREST, bVerbose, wi);
3232 add_wall_energrps(groups, ir->nwall, symtab);
3233 ir->opts.ngener = groups->grps[egcENER].nr;
3234 nvcm = str_nelem(vcm, MAXPTR, ptr1);
3236 do_numbering(natoms, groups, nvcm, ptr1, grps, gnames, egcVCM,
3237 restnm, nvcm == 0 ? egrptpALL_GENREST : egrptpPART, bVerbose, wi);
3240 warning(wi, "Some atoms are not part of any center of mass motion removal group.\n"
3241 "This may lead to artifacts.\n"
3242 "In most cases one should use one group for the whole system.");
3245 /* Now we have filled the freeze struct, so we can calculate NRDF */
3246 calc_nrdf(mtop, ir, gnames);
3252 /* Must check per group! */
3253 for (i = 0; (i < ir->opts.ngtc); i++)
3255 ntot += ir->opts.nrdf[i];
3257 if (ntot != (DIM*natoms))
3259 fac = sqrt(ntot/(DIM*natoms));
3262 fprintf(stderr, "Scaling velocities by a factor of %.3f to account for constraints\n"
3263 "and removal of center of mass motion\n", fac);
3265 for (i = 0; (i < natoms); i++)
3267 svmul(fac, v[i], v[i]);
3272 nuser = str_nelem(user1, MAXPTR, ptr1);
3273 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser1,
3274 restnm, egrptpALL_GENREST, bVerbose, wi);
3275 nuser = str_nelem(user2, MAXPTR, ptr1);
3276 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcUser2,
3277 restnm, egrptpALL_GENREST, bVerbose, wi);
3278 nuser = str_nelem(xtc_grps, MAXPTR, ptr1);
3279 do_numbering(natoms, groups, nuser, ptr1, grps, gnames, egcXTC,
3280 restnm, egrptpONE, bVerbose, wi);
3281 nofg = str_nelem(orirefitgrp, MAXPTR, ptr1);
3282 do_numbering(natoms, groups, nofg, ptr1, grps, gnames, egcORFIT,
3283 restnm, egrptpALL_GENREST, bVerbose, wi);
3285 /* QMMM input processing */
3286 nQMg = str_nelem(QMMM, MAXPTR, ptr1);
3287 nQMmethod = str_nelem(QMmethod, MAXPTR, ptr2);
3288 nQMbasis = str_nelem(QMbasis, MAXPTR, ptr3);
3289 if ((nQMmethod != nQMg) || (nQMbasis != nQMg))
3291 gmx_fatal(FARGS, "Invalid QMMM input: %d groups %d basissets"
3292 " and %d methods\n", nQMg, nQMbasis, nQMmethod);
3294 /* group rest, if any, is always MM! */
3295 do_numbering(natoms, groups, nQMg, ptr1, grps, gnames, egcQMMM,
3296 restnm, egrptpALL_GENREST, bVerbose, wi);
3297 nr = nQMg; /*atoms->grps[egcQMMM].nr;*/
3298 ir->opts.ngQM = nQMg;
3299 snew(ir->opts.QMmethod, nr);
3300 snew(ir->opts.QMbasis, nr);
3301 for (i = 0; i < nr; i++)
3303 /* input consists of strings: RHF CASSCF PM3 .. These need to be
3304 * converted to the corresponding enum in names.c
3306 ir->opts.QMmethod[i] = search_QMstring(ptr2[i], eQMmethodNR,
3308 ir->opts.QMbasis[i] = search_QMstring(ptr3[i], eQMbasisNR,
3312 nQMmult = str_nelem(QMmult, MAXPTR, ptr1);
3313 nQMcharge = str_nelem(QMcharge, MAXPTR, ptr2);
3314 nbSH = str_nelem(bSH, MAXPTR, ptr3);
3315 snew(ir->opts.QMmult, nr);
3316 snew(ir->opts.QMcharge, nr);
3317 snew(ir->opts.bSH, nr);
3319 for (i = 0; i < nr; i++)
3321 ir->opts.QMmult[i] = strtol(ptr1[i], NULL, 10);
3322 ir->opts.QMcharge[i] = strtol(ptr2[i], NULL, 10);
3323 ir->opts.bSH[i] = (gmx_strncasecmp(ptr3[i], "Y", 1) == 0);
3326 nCASelec = str_nelem(CASelectrons, MAXPTR, ptr1);
3327 nCASorb = str_nelem(CASorbitals, MAXPTR, ptr2);
3328 snew(ir->opts.CASelectrons, nr);
3329 snew(ir->opts.CASorbitals, nr);
3330 for (i = 0; i < nr; i++)
3332 ir->opts.CASelectrons[i] = strtol(ptr1[i], NULL, 10);
3333 ir->opts.CASorbitals[i] = strtol(ptr2[i], NULL, 10);
3335 /* special optimization options */
3337 nbOPT = str_nelem(bOPT, MAXPTR, ptr1);
3338 nbTS = str_nelem(bTS, MAXPTR, ptr2);
3339 snew(ir->opts.bOPT, nr);
3340 snew(ir->opts.bTS, nr);
3341 for (i = 0; i < nr; i++)
3343 ir->opts.bOPT[i] = (gmx_strncasecmp(ptr1[i], "Y", 1) == 0);
3344 ir->opts.bTS[i] = (gmx_strncasecmp(ptr2[i], "Y", 1) == 0);
3346 nSAon = str_nelem(SAon, MAXPTR, ptr1);
3347 nSAoff = str_nelem(SAoff, MAXPTR, ptr2);
3348 nSAsteps = str_nelem(SAsteps, MAXPTR, ptr3);
3349 snew(ir->opts.SAon, nr);
3350 snew(ir->opts.SAoff, nr);
3351 snew(ir->opts.SAsteps, nr);
3353 for (i = 0; i < nr; i++)
3355 ir->opts.SAon[i] = strtod(ptr1[i], NULL);
3356 ir->opts.SAoff[i] = strtod(ptr2[i], NULL);
3357 ir->opts.SAsteps[i] = strtol(ptr3[i], NULL, 10);
3359 /* end of QMMM input */
3363 for (i = 0; (i < egcNR); i++)
3365 fprintf(stderr, "%-16s has %d element(s):", gtypes[i], groups->grps[i].nr);
3366 for (j = 0; (j < groups->grps[i].nr); j++)
3368 fprintf(stderr, " %s", *(groups->grpname[groups->grps[i].nm_ind[j]]));
3370 fprintf(stderr, "\n");
3374 nr = groups->grps[egcENER].nr;
3375 snew(ir->opts.egp_flags, nr*nr);
3377 bExcl = do_egp_flag(ir, groups, "energygrp-excl", egpexcl, EGP_EXCL);
3378 if (bExcl && ir->cutoff_scheme == ecutsVERLET)
3380 warning_error(wi, "Energy group exclusions are not (yet) implemented for the Verlet scheme");
3382 if (bExcl && EEL_FULL(ir->coulombtype))
3384 warning(wi, "Can not exclude the lattice Coulomb energy between energy groups");
3387 bTable = do_egp_flag(ir, groups, "energygrp-table", egptable, EGP_TABLE);
3388 if (bTable && !(ir->vdwtype == evdwUSER) &&
3389 !(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
3390 !(ir->coulombtype == eelPMEUSERSWITCH))
3392 gmx_fatal(FARGS, "Can only have energy group pair tables in combination with user tables for VdW and/or Coulomb");
3395 decode_cos(efield_x, &(ir->ex[XX]), FALSE);
3396 decode_cos(efield_xt, &(ir->et[XX]), TRUE);
3397 decode_cos(efield_y, &(ir->ex[YY]), FALSE);
3398 decode_cos(efield_yt, &(ir->et[YY]), TRUE);
3399 decode_cos(efield_z, &(ir->ex[ZZ]), FALSE);
3400 decode_cos(efield_zt, &(ir->et[ZZ]), TRUE);
3404 do_adress_index(ir->adress, groups, gnames, &(ir->opts), wi);
3407 for (i = 0; (i < grps->nr); i++)
3419 static void check_disre(gmx_mtop_t *mtop)
3421 gmx_ffparams_t *ffparams;
3422 t_functype *functype;
3424 int i, ndouble, ftype;
3425 int label, old_label;
3427 if (gmx_mtop_ftype_count(mtop, F_DISRES) > 0)
3429 ffparams = &mtop->ffparams;
3430 functype = ffparams->functype;
3431 ip = ffparams->iparams;
3434 for (i = 0; i < ffparams->ntypes; i++)
3436 ftype = functype[i];
3437 if (ftype == F_DISRES)
3439 label = ip[i].disres.label;
3440 if (label == old_label)
3442 fprintf(stderr, "Distance restraint index %d occurs twice\n", label);
3450 gmx_fatal(FARGS, "Found %d double distance restraint indices,\n"
3451 "probably the parameters for multiple pairs in one restraint "
3452 "are not identical\n", ndouble);
3457 static gmx_bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys,
3458 gmx_bool posres_only,
3462 gmx_mtop_ilistloop_t iloop;
3472 for (d = 0; d < DIM; d++)
3474 AbsRef[d] = (d < ndof_com(ir) ? 0 : 1);
3476 /* Check for freeze groups */
3477 for (g = 0; g < ir->opts.ngfrz; g++)
3479 for (d = 0; d < DIM; d++)
3481 if (ir->opts.nFreeze[g][d] != 0)
3489 /* Check for position restraints */
3490 iloop = gmx_mtop_ilistloop_init(sys);
3491 while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
3494 (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
3496 for (i = 0; i < ilist[F_POSRES].nr; i += 2)
3498 pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
3499 for (d = 0; d < DIM; d++)
3501 if (pr->posres.fcA[d] != 0)
3510 return (AbsRef[XX] != 0 && AbsRef[YY] != 0 && AbsRef[ZZ] != 0);
3513 void triple_check(const char *mdparin, t_inputrec *ir, gmx_mtop_t *sys,
3517 int i, m, g, nmol, npct;
3518 gmx_bool bCharge, bAcc;
3519 real gdt_max, *mgrp, mt;
3521 gmx_mtop_atomloop_block_t aloopb;
3522 gmx_mtop_atomloop_all_t aloop;
3525 char warn_buf[STRLEN];
3527 set_warning_line(wi, mdparin, -1);
3529 if (EI_DYNAMICS(ir->eI) && !EI_SD(ir->eI) && ir->eI != eiBD &&
3530 ir->comm_mode == ecmNO &&
3531 !(absolute_reference(ir, sys, FALSE, AbsRef) || ir->nsteps <= 10) &&
3532 !ETC_ANDERSEN(ir->etc))
3534 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");
3537 /* Check for pressure coupling with absolute position restraints */
3538 if (ir->epc != epcNO && ir->refcoord_scaling == erscNO)
3540 absolute_reference(ir, sys, TRUE, AbsRef);
3542 for (m = 0; m < DIM; m++)
3544 if (AbsRef[m] && norm2(ir->compress[m]) > 0)
3546 warning(wi, "You are using pressure coupling with absolute position restraints, this will give artifacts. Use the refcoord_scaling option.");
3554 aloopb = gmx_mtop_atomloop_block_init(sys);
3555 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
3557 if (atom->q != 0 || atom->qB != 0)
3565 if (EEL_FULL(ir->coulombtype))
3568 "You are using full electrostatics treatment %s for a system without charges.\n"
3569 "This costs a lot of performance for just processing zeros, consider using %s instead.\n",
3570 EELTYPE(ir->coulombtype), EELTYPE(eelCUT));
3571 warning(wi, err_buf);
3576 if (ir->coulombtype == eelCUT && ir->rcoulomb > 0 && !ir->implicit_solvent)
3579 "You are using a plain Coulomb cut-off, which might produce artifacts.\n"
3580 "You might want to consider using %s electrostatics.\n",
3582 warning_note(wi, err_buf);
3586 /* Generalized reaction field */
3587 if (ir->opts.ngtc == 0)
3589 sprintf(err_buf, "No temperature coupling while using coulombtype %s",
3591 CHECK(ir->coulombtype == eelGRF);
3595 sprintf(err_buf, "When using coulombtype = %s"
3596 " ref-t for temperature coupling should be > 0",
3598 CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
3601 if (ir->eI == eiSD1 &&
3602 (gmx_mtop_ftype_count(sys, F_CONSTR) > 0 ||
3603 gmx_mtop_ftype_count(sys, F_SETTLE) > 0))
3605 sprintf(warn_buf, "With constraints integrator %s is less accurate, consider using %s instead", ei_names[ir->eI], ei_names[eiSD2]);
3606 warning_note(wi, warn_buf);
3610 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3612 for (m = 0; (m < DIM); m++)
3614 if (fabs(ir->opts.acc[i][m]) > 1e-6)
3623 snew(mgrp, sys->groups.grps[egcACC].nr);
3624 aloop = gmx_mtop_atomloop_all_init(sys);
3625 while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
3627 mgrp[ggrpnr(&sys->groups, egcACC, i)] += atom->m;
3630 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3632 for (m = 0; (m < DIM); m++)
3634 acc[m] += ir->opts.acc[i][m]*mgrp[i];
3638 for (m = 0; (m < DIM); m++)
3640 if (fabs(acc[m]) > 1e-6)
3642 const char *dim[DIM] = { "X", "Y", "Z" };
3644 "Net Acceleration in %s direction, will %s be corrected\n",
3645 dim[m], ir->nstcomm != 0 ? "" : "not");
3646 if (ir->nstcomm != 0 && m < ndof_com(ir))
3649 for (i = 0; (i < sys->groups.grps[egcACC].nr); i++)
3651 ir->opts.acc[i][m] -= acc[m];
3659 if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
3660 !gmx_within_tol(sys->ffparams.reppow, 12.0, 10*GMX_DOUBLE_EPS))
3662 gmx_fatal(FARGS, "Soft-core interactions are only supported with VdW repulsion power 12");
3665 if (ir->ePull != epullNO)
3667 if (ir->pull->grp[0].nat == 0)
3669 absolute_reference(ir, sys, FALSE, AbsRef);
3670 for (m = 0; m < DIM; m++)
3672 if (ir->pull->dim[m] && !AbsRef[m])
3674 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.");
3680 if (ir->pull->eGeom == epullgDIRPBC)
3682 for (i = 0; i < 3; i++)
3684 for (m = 0; m <= i; m++)
3686 if ((ir->epc != epcNO && ir->compress[i][m] != 0) ||
3687 ir->deform[i][m] != 0)
3689 for (g = 1; g < ir->pull->ngrp; g++)
3691 if (ir->pull->grp[g].vec[m] != 0)
3693 gmx_fatal(FARGS, "Can not have dynamic box while using pull geometry '%s' (dim %c)", EPULLGEOM(ir->pull->eGeom), 'x'+m);
3705 void double_check(t_inputrec *ir, matrix box, gmx_bool bConstr, warninp_t wi)
3709 char warn_buf[STRLEN];
3712 ptr = check_box(ir->ePBC, box);
3715 warning_error(wi, ptr);
3718 if (bConstr && ir->eConstrAlg == econtSHAKE)
3720 if (ir->shake_tol <= 0.0)
3722 sprintf(warn_buf, "ERROR: shake-tol must be > 0 instead of %g\n",
3724 warning_error(wi, warn_buf);
3727 if (IR_TWINRANGE(*ir) && ir->nstlist > 1)
3729 sprintf(warn_buf, "With twin-range cut-off's and SHAKE the virial and the pressure are incorrect.");
3730 if (ir->epc == epcNO)
3732 warning(wi, warn_buf);
3736 warning_error(wi, warn_buf);
3741 if ( (ir->eConstrAlg == econtLINCS) && bConstr)
3743 /* If we have Lincs constraints: */
3744 if (ir->eI == eiMD && ir->etc == etcNO &&
3745 ir->eConstrAlg == econtLINCS && ir->nLincsIter == 1)
3747 sprintf(warn_buf, "For energy conservation with LINCS, lincs_iter should be 2 or larger.\n");
3748 warning_note(wi, warn_buf);
3751 if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder < 8))
3753 sprintf(warn_buf, "For accurate %s with LINCS constraints, lincs-order should be 8 or more.", ei_names[ir->eI]);
3754 warning_note(wi, warn_buf);
3756 if (ir->epc == epcMTTK)
3758 warning_error(wi, "MTTK not compatible with lincs -- use shake instead.");
3762 if (bConstr && ir->epc == epcMTTK)
3764 warning_note(wi, "MTTK with constraints is deprecated, and will be removed in GROMACS 5.1");
3767 if (ir->LincsWarnAngle > 90.0)
3769 sprintf(warn_buf, "lincs-warnangle can not be larger than 90 degrees, setting it to 90.\n");
3770 warning(wi, warn_buf);
3771 ir->LincsWarnAngle = 90.0;
3774 if (ir->ePBC != epbcNONE)
3776 if (ir->nstlist == 0)
3778 warning(wi, "With nstlist=0 atoms are only put into the box at step 0, therefore drifting atoms might cause the simulation to crash.");
3780 bTWIN = (ir->rlistlong > ir->rlist);
3781 if (ir->ns_type == ensGRID)
3783 if (sqr(ir->rlistlong) >= max_cutoff2(ir->ePBC, box))
3785 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",
3786 bTWIN ? (ir->rcoulomb == ir->rlistlong ? "rcoulomb" : "rvdw") : "rlist");
3787 warning_error(wi, warn_buf);
3792 min_size = min(box[XX][XX], min(box[YY][YY], box[ZZ][ZZ]));
3793 if (2*ir->rlistlong >= min_size)
3795 sprintf(warn_buf, "ERROR: One of the box lengths is smaller than twice the cut-off length. Increase the box size or decrease rlist.");
3796 warning_error(wi, warn_buf);
3799 fprintf(stderr, "Grid search might allow larger cut-off's than simple search with triclinic boxes.");
3806 void check_chargegroup_radii(const gmx_mtop_t *mtop, const t_inputrec *ir,
3810 real rvdw1, rvdw2, rcoul1, rcoul2;
3811 char warn_buf[STRLEN];
3813 calc_chargegroup_radii(mtop, x, &rvdw1, &rvdw2, &rcoul1, &rcoul2);
3817 printf("Largest charge group radii for Van der Waals: %5.3f, %5.3f nm\n",
3822 printf("Largest charge group radii for Coulomb: %5.3f, %5.3f nm\n",
3828 if (rvdw1 + rvdw2 > ir->rlist ||
3829 rcoul1 + rcoul2 > ir->rlist)
3832 "The sum of the two largest charge group radii (%f) "
3833 "is larger than rlist (%f)\n",
3834 max(rvdw1+rvdw2, rcoul1+rcoul2), ir->rlist);
3835 warning(wi, warn_buf);
3839 /* Here we do not use the zero at cut-off macro,
3840 * since user defined interactions might purposely
3841 * not be zero at the cut-off.
3843 if ((EVDW_IS_ZERO_AT_CUTOFF(ir->vdwtype) ||
3844 ir->vdw_modifier != eintmodNONE) &&
3845 rvdw1 + rvdw2 > ir->rlistlong - ir->rvdw)
3847 sprintf(warn_buf, "The sum of the two largest charge group "
3848 "radii (%f) is larger than %s (%f) - rvdw (%f).\n"
3849 "With exact cut-offs, better performance can be "
3850 "obtained with cutoff-scheme = %s, because it "
3851 "does not use charge groups at all.",
3853 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3854 ir->rlistlong, ir->rvdw,
3855 ecutscheme_names[ecutsVERLET]);
3858 warning(wi, warn_buf);
3862 warning_note(wi, warn_buf);
3865 if ((EEL_IS_ZERO_AT_CUTOFF(ir->coulombtype) ||
3866 ir->coulomb_modifier != eintmodNONE) &&
3867 rcoul1 + rcoul2 > ir->rlistlong - ir->rcoulomb)
3869 sprintf(warn_buf, "The sum of the two largest charge group radii (%f) is larger than %s (%f) - rcoulomb (%f).\n"
3870 "With exact cut-offs, better performance can be obtained with cutoff-scheme = %s, because it does not use charge groups at all.",
3872 ir->rlistlong > ir->rlist ? "rlistlong" : "rlist",
3873 ir->rlistlong, ir->rcoulomb,
3874 ecutscheme_names[ecutsVERLET]);
3877 warning(wi, warn_buf);
3881 warning_note(wi, warn_buf);